gfmul.c

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

/* ============================================================================
 *  Syntax:     c = gfmul(a, b, p)
 * GFMUL  Element-by-element multiply in GF(P) or GF(P^M).
 *      C = BFMUL(A, B) multiplies A with B in GF(2) element-by-element. A and
 *      B are scalar or vectors in GF(2).
 *
 *      C = BFMUL(A, B, P) multiplies A with B in GF(P) element-by-element
 *      when P is a prime number. In this case, A and B are scalars or vectors
 *      in GF(P). When P is a matrix that contains the list of all tuple
 *      elements in GF(p^M), this function multiplies A with B in GF(P^M). The
 *      output C is alpha^C = alpha^A * alpha^B in GF(Q^M). In this case A and
 *      B are scalars or vectors in GF(p^M).  You can generate the tuple of all
 *      elements in GF(Q^M) by P = GFTUPLE([-1:Q^M-2]', M, Q). 
 *
 *      The element in GF(p^M) is represented by power 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:31 $
 * ==========================================================================
 */
#include <math.h>
#include "mex.h"
#include "gflib.c"
void mexFunction(int nlhs, Matrix *plhs[], int nrhs, Matrix *prhs[])
{
    int     i, j, row, col, mat_flag;
    int     ma, na, mb, nb, np, mp, nc, len_a, len_b;
    int     *paa, *pbb, *pcc, *pp, *Iwork;
    double  *pa, *pb, *p, *pc;

    if ( nrhs < 2 ){
        mexErrMsgTxt("Not enough input for GFMUL!");
    }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]);
    len_a = ma*na;
    len_b = mb*nb;

    if( ma != 1 ){
        if( na != 1 ){
            if( ma != mb && na != nb)
                mexErrMsgTxt("The input matrix of GFMUL in GF() must be exactly matched!");
            else
                mat_flag = 1;
		}else{
            mat_flag = 0;
        }
    }else{
        mat_flag = 0;
    }
    if( mb != 1 ){
        if( nb != 1 ){
            if( mb != ma && nb != na)
                mexErrMsgTxt("The input matrix of GFMUL in GF() must be exactly matched!");
            else
                mat_flag = 1;
		}else{
            mat_flag = 0;
        }
    }else{
        mat_flag = 0;
    }

    /* the calculation lenth is max(len_a, len_b) */
    if ( mat_flag == 0 ){
        if (len_a >= len_b)
            nc = len_a;
        else
            nc = len_b;
    } else {
        row = ma;
        col = na;
        nc = row*col;
    }

    /* variable type conversion for calling functions in gflib.c */
    paa = (int *)mxCalloc(len_a, sizeof(int));
    pbb = (int *)mxCalloc(len_b, sizeof(int));
    if ( mat_flag == 0 ){
        for (i=0; i < len_a; i++)
            paa[i] = (int) pa[i];
        for (i=0; i < len_b; i++)
            pbb[i] = (int) pb[i];
    } else {
        for ( i=0; i < row; i++){
            for ( j=0; j < col; j++){
                paa[i+j*row] = (int)pa[i+j*row];
                pbb[i+j*row] = (int)pb[i+j*row];
            }
        }
    }
    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){
        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 gfmul() in gflib.c */
        pcc = (int *)mxCalloc(nc, sizeof(int));
        gfmul(paa, len_a, pbb, len_b, *pp, pcc);
        if ( mat_flag == 0 ){
            if ( ma == 1 ) 
                pc = mxGetPr(plhs[0]=mxCreateFull(1, nc, 0));
            else if ( na == 1 )
                pc = mxGetPr(plhs[0]=mxCreateFull(nc, 1, 0));
    
            for(i=0; i < nc; i++){
                if(pcc[i] == -Inf){
                    pc[i] = -mexGetInf();
                }else{
                    pc[i] = (double)pcc[i];
                }
            }
        } else {
            pc = mxGetPr(plhs[0]=mxCreateFull(row, col, 0));
            for ( i=0; i < row; i++){
                for ( j=0; j < col; j++){
                    if(pcc[i+j*row] == -Inf)
                        pc[i+j*row] = -mexGetInf();
                    else
                        pc[i+j*row] = (double)pcc[i+j*row];
                }
            }
        }
        return;
	} else {
        /* call gfpmul() in gflib.c */
        pcc = (int *)mxCalloc(nc, sizeof(int));
        Iwork = (int *)mxCalloc(2*nc+nc*mp+np+1, sizeof(int));
        Iwork[0] = nc;
        gfpmul(paa, len_a, pbb, len_b, pp, np, mp, pcc, Iwork, Iwork+1);
        if ( mat_flag == 0 ){
            if ( ma == 1 ) 
                pc = mxGetPr(plhs[0]=mxCreateFull(1, Iwork[0], 0));
            else if ( na == 1 )
                pc = mxGetPr(plhs[0]=mxCreateFull(Iwork[0], 1, 0));
    
            for(i=0; i < Iwork[0]; i++){
                if(pcc[i] == -Inf)
                    pc[i] = -mexGetInf();
                else
                    pc[i] = (double)pcc[i];
            }
        } else {
            pc = mxGetPr(plhs[0]=mxCreateFull(row, col, 0));
            for ( i=0; i < row; i++){
                for ( j=0; j < col; j++){
                    if ( i+j*row >= Iwork[0] ){
                        pcc[i+j*row] = 0;
                    } else {
                        if(pcc[i+j*row] == -Inf)
                            pc[i+j*row] = -mexGetInf();
                        else
                            pc[i+j*row] = (double)pcc[i+j*row];
                    }
                }
            }
        }
        return;
    }
}
/*--end of GFMUL.C--*/