gflib.c

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

/* ============================================================================
 * gflib.c is a library file including the following functions:
 *  P.S.: All the parameters in these functions have the type of integer. 
 *  static void gftrunc(*pa, *len_a, len_p, *Iwork)
 *          M_FILE: C = GFTRUNC(A); Iwork --- len_a
 *          truncates  the results of every GF's computation .
 *  static void gfpadd(*pa, len_a, *pb, len_b, *pp, np, mp, *pc, *nc, *Iwork)
 *          M_FILE: C = GFADD(A,B,P); Iwork --- *nc+*nc*mp+np
 *          computes A adds B in GF when P is a matrix.
 *  static void gfadd(*pa, len_a, *pb, len_b, pp, len, *pc, *nc, *Iwork)
 *          M_FILE: C = GFADD(A,B,P,LEN), Iwork --- *nc
 *          computes A adds B in GF when P is a scalar.
 *  static void gfpmul(*pa, len_a, *pb, len_b, *pp, np, mp, *pc, *nc, *Iwork)
 *          M_FILE: C = GFMUL(A,B,P); Iwork --- 2*nc+nc*mp+np
 *          computes A multiply B in GF(P) when P is a matrix.
 *  static void gfmul(*pa, len_a, *pb, len_b, pp, *pc)
 *          M_FILE: C = GFMUL(A,B,P); no Iwork 
 *          computes A multiply B in GF(P) when P is a scalar.
 *  static void gfconv(*pa, len_a, *pb, len_b, pp, *pc, *Iwork)
 *          M_FILE: C = GFCONV(A,B,P); Iwork --- 6*(len_b+len_a)-5
 *          computes the convolution between two GF(P) polynomials when P is a scalar.
 *  static void gfpconv(*pa, len_a, *pb, len_b, *pp, np, mp, *pc, *Iwork)
 *          M_FILE: C = GFCONV(A,B,P); Iwork --- 5+3*(np+mp)
 *          computes the convolution between two GF(P) when P is a Matrix. 
 *  static void gfdeconv(*pb, len_b, *pa, len_a, pp, *pq, len_q, *pr, *len_r, *Iwork)
 *          M_FILE: [Q, R] = GFDECONV(B,A,P); Iwork --- 14*len_b-2*len_a+2
 *          computes the deconvolution between two GF(P) polynomials when P is a scalar.
 *  static void gfpdeconv(*pb, len_b, *pa, len_a, *pp, np, mp, *pq, *pr, *len_r, *Iwork)
 *          M_FILE: [Q, R] = GFDECONV(B,A,P); Iwork --- (2+mp)*np+len_a-1
 *          computes the deconvolution between two GF(P) when P is a Matrix. 
 *  static void gfplus(*pi, mi, ni, *pj, mj, nj, *alpha, len_alpha, *beta, len_beta, *pk)
 *          M_FILE: C = GFPLUS(A,B,ALPHA,BETA) 
 *          Galois Field addition (vectorized).
 *  static void gffilter(*pb, len_b, *pa, len_a, *px, len_x, p, *pOut, *Iwork)
 *          M_FILE: Y = GFFILTER(B,A,X,P); Iwork --- 2*len_x+len_b+len_a-2
 *          filters the data in GF(P).
 *  static void flxor(*px, mx, nx, *py, my, ny, *pz, *mz, *nz)  
 *          M_FILE: Z = FLXOR(X,Y) 
 *          exclusive OR computation.
 *  static void fliplr(*pa, len_a, *Iwork)
 *          M_FILE: Y = FLIPLR(X)
 *          Iwork --- len_a
 *          flips matrix in the left/right direction.
 *  static int isprime(p) 
 *          M_FILE: ISP = ISPRIME(X) 
 *          return TRUE for prime numbers.
 *  static void errlocp1(*syndr,t,*pNum,*pInv,pow_dim,err_In,*Iwork,*sigma_Out,*len_Out)
 *          M_FILE: [SIGMA, ERR] = ERRLOCP(SYNDROME,T,TP,POW_DIM,ERR,1);
 *          Iwork --- 5*(2+t)+(t+4)*(t+1)   PS: t=2*T (T is in M_file)
 *  static void errlocp0(*syndr,t,*pNum,*pInv,pow_dim,err_In,*Iwork,*sigma_Out,*len_Out)
 *          M_FILE: [SIGMA, ERR] = ERRLOCP(SYNDROME,T,TP,POW_DIM,ERR,0);
 *          Iwork --- 5*(2+t)+(t+4)*(t+1)   PS: t=T (T is in M_file)
 *  static void bchcore(*code,pow_dim,dim,k,t,*pp,*Iwork,*msg,*err,*ccode)
 *          M_FILE: [SIGMA, ERR, CCODE] = BCHCORE(CODE,POW_DIM,DIM,K,T,P);
 *          Iwork --- (2+dim)*(2*pow_dim+1)+t*t+15*t+16
 *  static void rscore(*code,k,*pp,dim,pow_dim,*Iwork,*msg,*err,*ccode)
 *          M_FILE: [SIGMA, ERR, CCODE] = RSCORE(CODE,K,P,DIM,POW_DIM);
 *          Iwork --- (pow_dim-2*k+5*dim+18)*pow_dim+3*dim+k*(k-13)+27
 *  static void pNumInv(*pp, np, mp, prim, *pNum, *pInv, *Iwork)
 *          M_CODE: tp_num = tp * prim.^[0 : tp_m-1]'; -- vector pNum
 *                  tp_inv(tp_num+1) = 0:pow_dim; -- vector pInv
 *          Iwork --- np*mp
 *  static void bi2de(*pbi, np, mp, prim, *pde)
 * ============================================================================
 *
 *     Original designed by Wes Wang,
 *     Jun Wu,     The Mathworks, Inc.
 *     Jan-25, 1996
 *
 *     Copyright (c) 1996 by The MAthWorks, Inc.
 *     All Rights Reserved
 *
 * ===========================================================================
 */
#define Inf 32766

/*
 *bi2de()
 *  Iwork --- np*mp
 */
static void bi2de(pbi, np, mp, prim, pde)
    int *pbi, np,  mp, prim, *pde;
{
  int     i, j, k;
  
  for(i=0; i<np; i++)
    pde[i] = 0;
  
  for (i=0; i < np; i++){
    for (j=0; j < mp; j++){
      if (pbi[i+j*np] != 0){
	if(j > 0){
	  for (k=0; k < j; k++)
	    pbi[i+j*np] = pbi[i+j*np]*prim;
	}
      }
      pde[i] = pde[i] + pbi[i+j*np];
    }
  }
}
/*end of bi2de()*/
/*
 *fliplr()
 *  Iwork --- len_a
 */
static void fliplr(pa, len_a, Iwork)
     int *pa, len_a, *Iwork;
{
  int     i;
  int     *tmp;

  tmp = Iwork;
  for(i=0;i < len_a; i++)
    tmp[i]= pa[len_a-1-i];
  for(i=0;i<len_a;i++)
    pa[i]= tmp[i];
}
/*--end of FLIPLR()--*/
/*
 *pNumInv()
 *  Iwork --- np*mp
 */
static void pNumInv(p, np, mp, prim, pNum, pInv, Iwork)
     int *p, np, mp, prim, *pNum, *pInv, *Iwork;
{
  int     i, j, k;
  int     *pp;
  
  pp = Iwork;
  for(i=0; i < np*mp; i++)
    pp[i] = p[i];
  for(i=0; i<np; i++){
    pNum[i] = 0;
    pInv[i] = 0;
  }
  for (i=0; i < np; i++){
    for (j=0; j < mp; j++){
      if (pp[i+j*np] != 0){
	if(j > 0){
	  for (k=0; k < j; k++)
	    pp[i+j*np] = pp[i+j*np]*prim;
	}
      }
      pNum[i] = pNum[i] + pp[i+j*np];
    }
    pInv[pNum[i]] = i;
  }
}
/*--end of pNumInv()--*/
/*
 * Check if P is a prime number.
 *   Return 1 when p is a prime number and return 0 when it is not.
 */
static int isprime(p)
     int p;
{
  int i, Fox, lenFox;
  
  lenFox = 0;
  for ( i=2; i < p; i++){
    Fox = p % i;
    if(Fox == 0)
      lenFox++;
  }
  if(lenFox == 0)
    return(1);
  else
    return(0);
}
/*--end of ISPRIME()--*/
/*
 * GFFilter()
 *  Iwork --- 2*len_x+len_b+len_a-2
 *  Iwork = tmp_x
 *      + len_b+len_x-1 = tmp_y
 *          + len_a+len_x-1 = bottom of Iwork in gffilter()
 */
static void gffilter(pb, len_b, pa, len_a, px, len_x, p, pOut, Iwork)
     int *pb, len_b, *pa, len_a, *px, len_x, p, *pOut, *Iwork;
{
  int     i, j, q, indx;
  int     tmp_b, tmp_a, len_xx, len_yy;
  int     *tmp_x, *tmp_y;
  
  /* make a(1) to be one. */
  if (pa[0] == 0){
#ifdef MATLAB_MEX_FILE
    printf("First denominator filter coeficient must be non-zero!");
#endif
  } else if (pa[0] != 1){
    if( isprime(p) ){     /* by Fermat's theory (pp 96 in MacWilliams) 
			   * for any integer q less than p.
			   */
      q = 1;
      for(i=0; i < p-2; i++){
	    q = q*pa[0];
      }
      if (((q*pa[0]) % p) != 1){
#ifdef MATLAB_MEX_FILE
	    printf("Warning: The computation may be not accurate because a large integer\n");
	    printf("         is introduced in to the calculation.\n");
#endif
      }
    } else {
#ifdef MATLAB_MEX_FILE
      printf("P must be a prime in GFFILTER.");
#endif
    }
    /* made pa[0] = 1. */
    for(i=0; i < len_a; i++)
      pa[i] = (pa[i]*q) % p;
  }
  /* computation length adjustment */
  len_xx = len_b+len_x-1;
  len_yy = len_a+len_x-1;
  tmp_x = Iwork;
  tmp_y = Iwork + len_xx;
  for(i=len_b-1; i < len_xx; i++)
    tmp_x[i] = px[i-len_b+1];
  
  /*
   * flip the computation variable
   * difference equation iteration
   */
  len_a--; 
  for (i=0; i < len_x; i++){
    indx = i + len_a;
    tmp_a = 0;
    tmp_b = 0;
    if ( len_a > 0 ){
      /*y(indx) = fb * xx(i:len_b +i-1)' + fa * y(i:indx-1)'; */
      for (j=0; j < len_b; j++)
	tmp_b = tmp_b + pb[len_b-1-j]*tmp_x[i+j];
      for (j=0; j < len_a; j++)
	tmp_a = tmp_a + (p-pa[len_a-j])*tmp_y[i+j];
      tmp_y[indx] = tmp_b + tmp_a;    
    }else{
      /* y(indx) = fb * xx(i:len_b+i-1)'; */
      for (j=0; j < len_b; j++)
	    tmp_b = tmp_b + pb[len_b-1-j]*tmp_x[i+j];
      tmp_y[indx] = tmp_b;        
    }
    tmp_y[indx] = tmp_y[indx] % p;
  }
  /* remove the added zeros. */
  for(i=0; i < len_x; i++)
    pOut[i] = tmp_y[i+len_a];
}
/*--end of GFFILTER()--*/
/*
 * gftrunc(a)
 *      Iwork --- *len_a
 */
static void gftrunc(pa, len_a, len_p, Iwork)
     int *pa, *len_a, len_p, *Iwork;
{
  int     i, len_ind, cut_0;
  int     *ind;
  
  /*if isempty(a) 
   *      c = a;
   *else
   *    cut_0 = 1;
   *    if nargin > 1
   *        if length(tp) > 1
   *            cut_0 = 0;
   *        end;  
   *    end;
   */
  ind = Iwork;
  if ( len_p > 1 )
    cut_0 = 0;
  else
    cut_0 = 1;
  /*    if cut_0
   *        ind = find(a ~= 0);
   *    else
   *        ind = find(a >= 0);
   *    end;
   */
  len_ind = 0;
  if ( cut_0 != 0 ){               
    for (i=0; i < len_a[0]; i++){
      if ( pa[i] != 0 ){
	    ind[len_ind] = i;
	    len_ind++;
      }
    }
  } else {                        
    for (i=0; i < len_a[0]; i++){  
      if ( pa[i] >= 0 ){
	    ind[len_ind] = i;
	    len_ind++;
      }
    }
  }
  /*  if ~isempty(ind)
   *        c = a(1 : ind(length(ind)));
   *    else
   *        if cut_0
   *                    c = 0;
   *        else
   *            c = -Inf;
   *        end;
   *     end;
   *end;
   */
  if ( len_ind > 0 ){
    len_a[0] = ind[len_ind-1]+1;
  } else {
    if(cut_0 != 0){
      pa[0] = 0;
      len_a[0] = 1;
    }else{
      pa[0] = -Inf;
      len_a[0] = 1;
    }
  }
}
/*--end of GFTRUNC()--*/
/*
 * gfpadd(a, b, p) adds two GF(P) polynomials A and B.
 *      When P is a matrix.
 *      Iwork --- nc+nc*mp+np
 *      Iwork = indx
 *          + *nc = tmp
 *              + (*nc)*mp = sum
 *                  + np = bottom of Iwork
 */
static void gfpadd(pa, len_a, pb, len_b, pp, np, mp, pc, nc, Iwork)
     int *pa, len_a, *pb, len_b, *pp, np, mp, *pc, *nc, *Iwork;
{
  int     prim, cal_len;
  int     i, j, k, minus_a, minus_b, len_indx, len_indx_a, len_indx_b;
  int     *indx, *tmp, *sum;
  
  indx = Iwork;
  tmp = indx + nc[0];
  sum = tmp + nc[0]*mp;
  /* GF(Q^M) field calculation.
   * The first row of P is -Inf, second is 1 and so on.
   */
  /*  handle vector case
   *    indx = find(a > n_p - 2);
   *    if ~isempty(indx)
   *        a(indx) = rem(a(indx), n_p - 1);
   *    end;
   *    indx = find(b > n_p - 2);
   *    if ~isempty(indx)
   *        b(indx) = rem(b(indx), n_p - 1);
   *    end;
   */
  for(i=0; i < len_a; i++){
    if(pa[i] > np-2)
      pa[i] = pa[i] % (np - 1);    
  }
  for(i=0; i < len_b; i++){
    if(pb[i] > np-2)
      pb[i] = pb[i] % (np - 1);
  }
  /*    if (a < 0)
   *        c = b;
   *    elseif (b < 0)
   *        c = a;
   *    else
   */
  minus_a = 0;
  minus_b = 0;
  for(i=0; i < len_a; i++){                 
    if(pa[i] < 0)                         
      minus_a++;                                          
  }
  for(i=0; i < len_b; i++){             
    if(pb[i] < 0)
      minus_b++;
  }
  if( minus_a == len_a ){
    nc[0] = len_b;
    for(i=0; i < len_b; i++)
      pc[i] = pb[i];
  } else if(minus_b == len_b){
    nc[0] = len_a;
    for(i=0; i < len_a; i++)
      pc[i] = pa[i];
  } else {
    /*        prim = max(max(p)) + 1;
     *        indx = find(~(a >= 0));
     *        a(indx) = - ones(1, length(indx));
     *        indx = find(~(b >= 0));
     *        b(indx) = - ones(1, length(indx));
     *        len_a = length(a);
     *        len_b = length(b);
     */
    prim = 0;
    for(i=0; i < np*mp; i++){
      if(prim < pp[i])                              
	    prim = pp[i] + 1;           
    }
    for(i=0; i < len_a; i++){
      if(pa[i] < 0)
	    pa[i] = -1;
    }
    for(i=0; i < len_b; i++){
      if(pb[i] < 0)
	    pb[i] = -1;
    }
    /*        if (len_a == len_b) 
     *            tmp = rem(p(a + 2, :) + p(b + 2, :), prim);
     *            cal_len = len_a;
     *        elseif  (len_a == 1)
     *            tmp = rem(ones(len_b,1)*p(a + 2, :) + p(b + 2, :), prim);
     *            cal_len = len_b;
     *        elseif  (len_b == 1)
     *            tmp = rem(p(a + 2, :) + ones(len_a, 1)*p(b + 2, :), prim);
     *            cal_len = len_a;
     *        else
     *            cal_len = min(len_a, len_b);
     *            tmp = rem(p(a(1:cal_len) + 2, :) + p(b(1:cal_len) + 2, :), prim);
     *        end;
     */
    if( len_a == len_b && len_a != 1){
      for(i=0; i < len_a; i++){
	    for(j=0; j < mp; j++){
	      tmp[i+j*nc[0]] = (pp[j*np+pa[i]+1] + pp[j*np+pb[i]+1]) % prim;
	      cal_len = len_a;
	    }
      }
    }else if( len_a ==1 && len_b != 1){
      for(i=0; i < len_b; i++){   
	    for(j=0; j < mp; j++){ 
	      tmp[i+j*nc[0]] = (pp[j*np+pa[0]+1]+pp[j*np+pb[i]+1]) % prim;
	      cal_len = len_b;
	    }
      }
    }else if(len_b == 1 && len_a != 1){
      for(i=0; i < len_a; i++){
	for(j=0; j < mp; j++){
	  tmp[i+j*nc[0]] = (pp[j*np+pa[i]+1]+pp[j*np+pb[0]+1]) % prim;
	  cal_len = len_a;             
	}
      }
    }else if( len_b == len_a && len_a == 1){
      for(j=0; j < mp; j++){
	tmp[j*nc[0]] = (pp[j*np+pa[0]+1]+pp[j*np+pb[0]+1]) % prim;
	cal_len = 1;
      }
    }else{
      if (len_a >= len_b)
	cal_len = len_b;                
      else
	cal_len = len_a;
      for(i=0; i < cal_len; i++){
	for(j=0; j < mp; j++)                   
	  tmp[i+j*nc[0]] =(pp[j*np+pa[i] + 1] + pp[j*np + pb[i] + 1]) % prim;
      }
    }
    for(i=0; i < cal_len; i++){
      len_indx = 0;
      for(j=0; j < np; j++){
	sum[j] = 0;
	for(k=0; k < mp; k++)
	  sum [j] = sum[j] + (pp[j+k*np]+prim-tmp[i+k*nc[0]]) % prim;
	
	if( sum[j]==0 ){
	  indx[len_indx] = j;
	  len_indx++;
	}
      }
      if (len_indx == 1){
	pc[i] = indx[0] - 1;
      } else {
#ifdef MATLAB_MEX_FILE
	if( len_indx == 0 )
	  printf("The list of Galois field is not a complete one.\n");
	else
	  printf("The list of Galois field has repeat element.\n");
#endif
      }
    }
    if( cal_len < len_a ){
      for(i=0; i<cal_len; i++)                
	pc[i] = pc[i];
      for(i=cal_len; i < len_a; i++)
	pc[i] = pa[i];
    }else if(cal_len < len_b){
      for(i=0; i<cal_len; i++)
	pc[i] = pc[i];                      
      for(i=cal_len; i < len_b; i++)
	pc[i] = pb[i];                      
    }
  }
  /*    indx = find(c < 0);
   *    if ~isempty(indx)
   *        c(indx) = indx - Inf;
   *    end;
   *    c = c(:)';
   *    if cal_len < len_a
   *        c = [c a(cal_len + 1 : len_a)];
   *    elseif cal_len < len_b
   *        c = [c b(cal_len + 1 : len_b)];
   *    end;
   *end;
   */
  for(i=0; i < nc[0]; i++){
    if( pc[i] < 0 )
      pc[i] = -Inf;
  }
}
/*--end of GFpADD()--*/
/*
 * 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.
 *  gfadd(a, b, p) adds two polynomials A and B in GF(P) when P is
 *          a scalar prime number.
 *    Iwork --- *nc
 */
static void gfadd(pa, len_a, pb, len_b, pp, len, pc, nc, Iwork)
     int *pa, len_a, *pb, len_b, pp, len, *pc, *nc, *Iwork;
{
  int i, indx;

  /*    % match the length of vectors.
   *    na = length(a);
   *    nb = length(b);
   *    if na > nb