gfadd.c

数字通信第四版原书的例程

/*=============================================================================
 *       Syntax:     c = gfadd(a, b, p, len)
 * GFADD  Add two GF(P) polynomials or two GF(P^M) elements.
 *       C = GFADD(A, B) adds two GF(2) polynomials A and B. The resulting
 *       GF(2) polynomial C will keep the larger length between A and B.
 *
 *       C = GFADD(A, B, P) adds two polynomials A and B in GF(P) when P is
 *       a scalar prime number.
 *       When P is a matrix that contains the tuple of all elements in 
 *       GF(Q^M), this function takes A and B as indices (power number of the
 *       exponential form) of GF(Q^M) elements, the output C is 
 *       alpha^C = alpha^A + alpha^B in GF(Q^M). The computation is 
 *       element-by-element computation. You can generate the tuple of all
 *       elements in GF(Q^M) by P = GFTUPLE([-1:Q^M-2]', M, Q).  In power
 *       representation form, [-Inf, 0, 1, 2, ...] represents 
 *       [0, 1, alpha, alpha^2, ...] in GF(p^m).
 *
 *       C = GFADD(A, B, P, LEN) adds two GF(P) polynomials A and B. The
 *       resulting GF(P) polynomial C keeps the given length LEN. When LEN is a
 *       negative number, the length of C equals degree(C) + 1. P must a scalar
 *       integer in using this format.
 *
 *       In polynomial computation, A, B, and C are in ascending order, i.e.,
 *       A = [a_0, a_1, a_2,..., a_(n-1), a_n] represents
 *       A(X) = a_0 + a_1 X + a_2 X^2 +...+ a_(n-1) X^(n-1) + a_n X^n
 *       a_i must be a element in GF(P).
 *
 *       In power representation form, [-Inf, 0, 1, 2, ...] represents 
 *       [0, 1, alpha, alpha^2, ...] in GF(p^m).
 *
 *       See also GFMUL, GFDIV, GFTUPLE, GFPLUS.
 *=============================================================================
 *     Original designed by Wes Wang,
 *     Jun Wu,     The Mathworks, Inc.
 *     Dec-12, 1995
 *
 *     Copyright (c) 1995-96 by The MAthWorks, Inc.
 *     All Rights Reserved
 *     $Revision: 1.1 $  $Date: 1996/04/01 18:14:02 $
 *===========================================================================
 */
#include <math.h>

#ifdef MATLAB_MEX_FILE
#include "mex.h"
#endif

#include "gflib.c"
void mexFunction(int nlhs, Matrix *plhs[], int nrhs, Matrix *prhs[])
{
    int     i, j, maxrow, maxcol;
    int     ma, na, mb, nb, nc, np, mp, len, len_a, len_b;
    int     *paa, *pbb, *pp, *pcc, *Iwork;
    double  *pa, *pb, *p, *pc;

    if ( nrhs < 2 ){
        mexErrMsgTxt("Not enough input for GFADD!");
    }else if ( nrhs == 2 ){
        np= 1;
        mp= 1;
    }else if ( nrhs >= 2 ){
        p = mxGetPr(prhs[2]);
        np= mxGetM(prhs[2]);
        mp= mxGetN(prhs[2]);
    }
    
    /* get input arguments */
    pa = mxGetPr(prhs[0]);
    pb = mxGetPr(prhs[1]);
    ma = mxGetM(prhs[0]);
    na = mxGetN(prhs[0]);
    mb = mxGetM(prhs[1]);
    nb = mxGetN(prhs[1]);

    if ( nrhs > 3 )
        len = (int)mxGetScalar(prhs[3]);
	else
        len = 0;

    if ( ma >= mb )
        maxrow = ma;
    else
        maxrow = mb;

    if ( na >= nb )
        maxcol = na;
    else
        maxcol = nb;

    /* the calculation lenth is max(len_a, len_b) */
    nc = maxrow*maxcol;
    len_a = nc;
    len_b = nc;

    /* % match the length of vectors. */
    if( nrhs > 3 && np*mp <= 1 )
        nc = len;
    else
        nc = nc;

    /* variable type conversion for calling functions in gflib.c */
    paa = (int *)mxCalloc(len_a, sizeof(int));
    pbb = (int *)mxCalloc(len_b, sizeof(int));
    pp = (int *)mxCalloc(np*mp, sizeof(int));
    if( nrhs > 2 ){
        for (i=0; i < np*mp; i++)
            pp[i] = (int) p[i];
    }else{
        pp[0] = 2;
    }

    /* computation */    
    if (np <= 1){  
        /* input check up */
        for (i=0; i < ma*na; i++){
            if ( pa[i] < 0 || pa[i] >= pp[0] )
                mexErrMsgTxt(" polynomial coeficients must be in GF(P)");
        }
        for (i=0; i < mb*nb; i++){
            if (pb[i] < 0 || pb[i] >= pp[0])
                mexErrMsgTxt("The polynomial coeficients must be in GF(P)");
	    }
        for (i=0; i < maxrow; i++){
            for(j=0; j < maxcol; j++){
                if( i >= ma || j >= na ){
                    paa[i+j*maxrow] = 0;
                } else {
                    paa[i+j*maxrow] = (int)pa[i+j*ma];
                }
            }
        }
        for (i=0; i < maxrow; i++){
            for(j=0; j < maxcol; j++){
                if( i >= mb || j >= nb ){
                    pbb[i+j*maxrow] = 0;
                } else {
                    pbb[i+j*maxrow] = (int)pb[i+j*mb];
                }
            }
        }
        pcc = (int *)mxCalloc(nc, sizeof(int));
        Iwork = (int *)mxCalloc(nc+1, sizeof(int));
        Iwork[0] = nc;
        gfadd(paa, len_a, pbb, len_b, *pp, len, pcc, Iwork, Iwork+1);
	} else {
        /* call gfpadd() in gflib.c */
        for (i=0; i < maxrow; i++){
            for(j=0; j < maxcol; j++){
                if( i >= ma || j >= na ){
                    paa[i+j*maxrow] = -1;
				} else {
                    if ( pa[i+j*maxrow] < 0 ){
                        paa[i+j*maxrow] = -1;
                    } else {
                        paa[i+j*maxrow] = (int)pa[i+j*ma];
                    }
    			}
    		}
		}
        for (i=0; i < maxrow; i++){
            for(j=0; j < maxcol; j++){
                if( i >= mb || j >= nb ){
                    pbb[i+j*maxrow] = -1;
				} else {
                    if ( pb[i+j*maxrow] < 0 ){
                        pbb[i+j*maxrow] = -1;
                    } else {
                        pbb[i+j*maxrow] = (int)pb[i+j*mb];
                    }
    			}
    		}
		}
        pcc = (int *)mxCalloc(nc, sizeof(int));        	    
        Iwork = (int *)mxCalloc(nc+nc*mp+np+1, sizeof(int));
        Iwork[0] = nc;
        gfpadd(paa, len_a, pbb, len_b, pp, np, mp, pcc, Iwork, Iwork+1);
    }
    if ( nrhs > 3 && np <= 1 ){ 
        pc = mxGetPr(plhs[0]=mxCreateFull(1, nc, 0));
        for(i=0; i < nc; i++){
            if( i < maxrow*maxcol ){
                pc[i] = (double)pcc[i];
            } else {
                pc[i] = 0;
            }
        }
    } else {
        pc = mxGetPr(plhs[0]=mxCreateFull(maxrow, maxcol, 0));
        for(i=0; i < maxrow; i++){
            for(j=0; j < maxcol; j++){
                if( pcc[i+j*maxrow] < 0 )
                    pc[i+j*maxrow] = -mexGetInf();
                else
                    pc[i+j*maxrow] = (double)pcc[i+j*maxrow];
            }
        }
    }
    return;
}
/*--end of GFADD.C--*/