gfconv.c

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

/*=============================================================================
 *       Syntax:     c = gfconv(a, b, p)
 * GFCONV GF(P) polynomial convolution or GF(P^M) element multiplication. 
 *      C = GFCONV(A, B) computes the convolution between two GF(2) 
 *      polynomials A and B. The polynomial degree of the resulted GF(2) 
 *      polynomial C equals degree(A) + degree(B).
 *
 *      C = GFADD(A, B, P) computes the convolution between two GF(P)
 *          polynomials 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 polynomial computation, A, B, and C are ascending ordered, 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 GFADD, GFDIV, GFTUPLE
 *
 *=============================================================================
 *     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:08 $
 *===========================================================================
 */
#include <math.h>
#include "mex.h"
#include "gflib.c"
void mexFunction(int nlhs, Matrix *plhs[], int nrhs, Matrix *prhs[])
{
    int     ma, na, mb, nb, nc, np, mp, len_a, len_b, len_p, i;
    int     *paa, *pbb, *pp, *pcc, *Iwork1, *Iwork2;
    double  *pa, *pb, *p, *pc;

    if ( nrhs < 2 ){
        mexErrMsgTxt("Not enough input for GFCONV!");
    }else if ( nrhs == 2 ){
        len_p = 1;
    }else if ( nrhs > 2 ){
        p = mxGetPr(prhs[2]);
        np= mxGetM(prhs[2]);
        mp= mxGetN(prhs[2]);
        len_p = np*mp;
    }
    /* 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]);
    len_a = ma*na;
    len_b = mb*nb;

    /* 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(len_p, sizeof(int));
    for (i=0; i < len_a; i++)
        paa[i] = (int) pa[i];        
    for (i=0; i < len_b; i++)
        pbb[i] = (int) pb[i];
    if( nrhs == 2 ){
        pp[0] = 2;
	}else{
        for (i=0; i < len_p; i++)
            pp[i] = (int) p[i];
    }

	/* truncate input */
	Iwork1 = (int *)mxCalloc(len_a+len_b, sizeof(int));
    gftrunc(paa,&len_a,len_p,Iwork1); 
    gftrunc(pbb,&len_b,len_p,Iwork1+len_a); 

    /* computation */    
    if (len_p <= 1){  /* case of polynomial calculation */
        /* input check up */
        for (i=0; i < len_a; i++){
            if (pa[i] < 0 || pa[i] != floor(pa[i]) || pa[i] >= pp[0] )
                mexErrMsgTxt("The polynomial coeficients must be in GF(P)");
        }
        for (i=0; i < len_b; i++){
            if (pb[i] < 0 || pb[i] != floor(pb[i]) || pb[i] >= pp[0])
                mexErrMsgTxt("The polynomial coeficients must be in GF(P)");
	    }
        /* call gfconv() in gflib.c */
        nc = len_a+len_b-1;
        pcc = (int *)mxCalloc(nc, sizeof(int));
        Iwork2 = (int *)mxCalloc(6*nc+2, sizeof(int));
        gfconv(paa, len_a, pbb, len_b, *pp, pcc, Iwork2);
   } else { /* computation in GF(p^m) field */   
        /* call gfpconv() in gflib.c */
        nc = len_a+len_b+1;
        pcc = (int *)mxCalloc(nc, sizeof(int));
        Iwork2 = (int *)mxCalloc(5+3*(np+mp), sizeof(int));        
        gfpconv(paa, len_a, pbb, len_b, pp, np, mp, pcc, Iwork2);
    }
    pc = mxGetPr(plhs[0]=mxCreateFull(1, nc, 0));
    for(i=0; i < nc; i++){
        if(pcc[i] < 0){
            pc[i] = -mexGetInf();
        }else{
            pc[i] = (double)pcc[i];
        }
    }
    return;
}
/*--end of GFCONV.C--*/