gflib.c

Proakis《contemporarycommunication systems using matlab》matlab源代码

   *                c(i+j-1) = gfmul(a(i), b(j), tp);
   *            end;  
   *        end;
   *    end;
   *end;
   */
  One = 1;
  for (i=0; i < len_a; i++){
    for (j=0; j < len_b; j++){
      if (pc[i+j] >= 0){
	gfpmul(&pa[i],1,&pb[j],1,pp,np,mp,tmp,&One,Iwork);
	gfpadd(&pc[i+j],1,tmp,1,pp,np,mp,&pc[i+j],&One,Iwork+2+mp+np);
      }else{
	gfpmul(&pa[i],1,&pb[j],1,pp,np,mp,&pc[i+j],&One,Iwork+3+2*mp+2*np);
      }
    }
  }
}
/*--end of GFpCONV()--*/
/*
 * GFdeconv() when P is a scalar.
 *  Iwork --- 14*len_b-2*len_a+3
 *  Iwork = pfb
 *      + len_b = pfa
 *          + len_a = px
 *              + len_b-len_a+1 = Iwork for gffilter()
 *                  + 3*len_b-len_a = Iwork for fliplr()
 *                      + len_b-len_a+1 = tmp
 *                          + len_b = Iwork for gfconv()
 *                              + 6*len_b+2 = Iwork for gfadd()
 *                                  + len_b = bottom of Iwork
 */
static void gfdeconv(pb, len_b, pa, len_a, pp, pq, len_q, pr, len_r, Iwork)
     int *pb, len_b, *pa, len_a, pp, *pq, len_q, *pr, *len_r, *Iwork;
{
  int  i, t, len_x;
  int  *pfb, *pfa,*px, *tmp;

  /* handle the case: len_a > len_b  */
  if (len_a > len_b){
    *pq = 0;
    for(i=0; i < len_b; i++)
      pr[i] = pb[i];
  }else{
    /* make a(len_a) to be one. */
    if (pa[len_a-1] != 1){
      if ( isprime(pp) ){/* by Fermat's theory (pp 96 in MacWilliams) 
			  * for any integer t less than p, t^(p-1) = 1. 
			  */
	t = 1;
	for(i=0; i < pp-2; i++)
	  t = t*pa[len_a-1];
#ifdef MATLAB_MEX_FILE
	if (((t * pa[len_a-1]) % pp) != 1){
	  printf("Warning: The computation may be not accurate because a large integer\n");
	  printf("         is introduced in to the calculation.\n");
	}
      } else {
	printf("P must be a prime in GFDECONV.\n");
      }
#endif
      /* made pa[len_a-1] = 1. */
      for( i=0; i < len_a; i++ )
	pa[i] = (pa[i]*t) % pp;
    }
    pfb = Iwork;
    pfa = pfb + len_b;
    for(i=0; i<len_b; i++)
      pfb[i] = pb[len_b-1-i]; 
    for(i=0; i<len_a; i++)
      pfa[i] = pa[len_a-1-i]; 
    len_x = len_b-len_a+1;
    px = pfa+len_a;
    px[0] = 1;
    
    gffilter(pfb, len_b, pfa, len_a, px, len_x, pp, pq, px+len_x);
    fliplr(pq, len_q, px+3*len_x+len_b+len_a-2);
    
    for (i=0; i < len_b; i++)
      pb[i] = (pp - pb[i]) % pp;
    /* r= gfadd(rem(p-b,p), gfconv(q,a,p),p,[]);*/            
    /* just take advantage of the memory of pb */
    tmp = px+3*len_x+2*len_b-1;
    gfconv(pq, len_q, pa, len_a, pp, tmp, tmp+len_b);
    gfadd(pb, len_b, tmp, len_b, pp, 0, pr, len_r, tmp+7*len_b+2);
  }
}
/*--end of GEDECONV()--*/    
/*
 * GFPLUS  Galois Field addition (vectorized)
 */
static void gfplus(pi, mi, ni, pj, mj, nj, alpha, len_alpha, beta, len_beta, pk)
     int     *pi, *pj, *alpha, *beta, *pk;
     int     mi, ni, mj, nj, len_alpha, len_beta;
{
  int     i, r, inci, incj, indi, indj, len_k;
  
  len_alpha = len_alpha - 1;
  if ((mi == 1) && (ni == 1) && (mj == 1) && (nj == 1)) {
    len_k = 1;
    inci = 0;
    incj = 0;
  } else if ((mi == mj) && (ni == nj)) {
    len_k = mi*ni;
    inci = 1;
    incj = 1;
  } else if ((mj == 1) && (nj == 1)) {
    len_k = mi*ni;
    inci = 1;
    incj = 0;
  } else if ((mi == 1) && (ni == 1)) {
    len_k = mj*nj;
    inci = 0;
    incj = 1;
  } else {
#ifdef MATLAB_MEX_FILE
    printf("Matrix dimensions must agree.\n");
#endif
  }
  for (i = 0; i < len_k; i++) {
    indi = (*pi < 0) ? 0 : *pi % len_alpha + 1;
    indj = (*pj < 0) ? 0 : *pj % len_alpha + 1;
    r = alpha[indi] ^ alpha[indj];
    pk[i] = beta[r]-1;
    pi += inci;
    pj += incj;
  }
}
/*--end of GFPLUS()--*/
/*
 * GFdeconv() when P is a matrix.
 *  Iwork --- (2+mp)*np+len_a-1
 *  Iwork = pNum
 *      + np = pInv
 *          + np = Iwork for pNumInv()
 *              + np*mp = Iwork for gftrunc()
 *                  + len_a-1 = bottom of Iwork in gfpdeconv()
 */
static void gfpdeconv(pb, len_b, pa, len_a, pp, np, mp, pq, pr, len_r, Iwork)
     int *pb, len_b, *pa, len_a, *pp, np, mp, *pq, *pr, *len_r, *Iwork;  
{
  int     i, j, k, len_q, prim, pow_dim, tmp, flag, factor;
  int     *pNum, *pInv;
  
  /* handle the case: len_a > len_b  */
  /*if len_a > len_b
   *    q = 0;
   *    if len_p > 1
   *        q = -Inf;
   *    end;
   *    r = b;
   *    return
   *end;
   */
  if (len_a > len_b){
    pq[0] = -Inf;
    for(i=0; i < len_b; i++)
      pr[i] = pb[i];
  }else{
    /*tp = p;
     *p = max(max(tp)) + 1; % the prime number.
     *[n_tp, m_tp] = size(tp);
     *pow_dim = p^m_tp - 1;
     *if pow_dim+1 ~= n_tp
     *    error('The given GF(P^M) M-tuple list is not a valid one');
     *end;
     */
    prim = 2;
    for(i=0; i < np*mp; i++){
      if(prim < pp[i])                              
	prim = pp[i] + 1;           
    }
    pow_dim = 1;
    for(i=0; i < mp; i++)
      pow_dim = pow_dim * prim ;
    
    if ( pow_dim != np ){
#ifdef MATLAB_MEX_FILE
      printf("The given GF(P^M) M-tuple list is not a valid one.\n");
#endif
    } else {
      pow_dim--;
      /*
       *tp_num = tp * 2.^[0 : m_tp-1]';
       *tp_inv(tp_num+1) = 0:pow_dim;
       */
      pNum = Iwork;
      pInv = Iwork+np;
      pNumInv(pp, np, mp, prim, pNum, pInv, pInv+np);
      /*
       *a = fliplr(a);
       *b = fliplr(b);
       *flag = 0;
       */
      flag = 0;
      if( pa[len_a-1] != 0 ){
	/* have to make a(1) = 0 *******/
	factor = pow_dim-pa[len_a-1];
	/*  indxa = find( a < 0 );
	 *  indxb = find( b < 0 );
	 *  a = rem(a + factor, pow_dim);
	 *  b = rem(b + factor, pow_dim);
	 *  flag_a1 = 1;
	 *  if ~isempty(indxa)
	 *      a(indxa) = indxa - Inf;
	 *  end;
	 *  if ~isempty(indxb)
	 *     b(indxb) = indxb - Inf;
	 *  end;
	 * end;
	 */
	for (i=0; i < len_a; i++){
	  if(pa[i] < 0){
	    pa[i] = -Inf;
	  }else{
	    pa[i] =(pa[i] + factor) % pow_dim;
	  }
	}
	for (i=0; i < len_b; i++){
	  if(pb[i] < 0){
	    pb[i] = -Inf;
	  }else{
	    pb[i] =(pb[i] + factor) % pow_dim;
	  }
	}
	flag = 1;
      }
      /*      len_q = len_b - len_a + 1;
       *      for i = 1 : len_q
       *          for k = 1 : len_a - 1
       *              if (b(i) < 0) | (a(k+1) < 0)
       *                  tmp = -1;
       *              else
       *                  tmp = rem(b(i) + a(k+1), pow_dim);
       *              end;
       *              b(i+k) = gfplus(b(i+k), tmp, tp_num, tp_inv);
       *          end;
       *      end;
       */
      len_q = len_b - len_a + 1; 
      for(i=0; i < len_q; i++){
	for(j=0; j < len_a-1; j++){
	  if(pb[len_b-1-i] < 0 || pa[len_a-2-j] < 0){
	    tmp = -1;
	  }else{
	    tmp = (pb[len_b-1-i] + pa[len_a-2-j]) % pow_dim;
	  }
	  gfplus(&pb[len_b-2-i-j], 1, 1, &tmp, 1, 1, pNum, np, pInv, np, &pb[len_b-2-i-j]);
	}
      }
      /*
       *    if len_a > 1
       *        r = gftrunc(fliplr(b(len_b-len_a+2 : len_b)), tp);
       *        if flag_a1
       *            r = rem(r + pow_dim-factor, pow_dim);
       *        end;
       *    else
       *        r = -Inf;
       *    end;
       *    q = fliplr(b(1:len_q));
       *end;
       */
      /* computes the quotient Q */ 
      for (i=len_a-1; i < len_b ; i++)
	pq[i-len_a+1] = pb[i];
      /* computes the remainder R */ 
      if (len_a > 1){
	len_r[0] = len_a-1;
	for(i=0; i <len_r[0]; i++)
	  pr[i] = pb[i];
	gftrunc(pr, len_r, np*mp, pInv+(mp+1)*np);
	
	if ( flag != 0 ){
	  for(i=0; i < len_r[0]; i++)
	    pr[i] = (pr[i] + pow_dim - factor) % pow_dim;
	}
      }else{
	pr[0] = -Inf;
	len_r[0] = 1;
      }
    }
  }
}
/*--end of gfpdeconv()--*/
static void gfdeconvp(pb, len_b, pa, len_a, pp, np, mp, pr, len_r, Iwork)
     int *pb, len_b, *pa, len_a, *pp, np, mp, *pr, *len_r, *Iwork;
{
  int     i, j, k, prim, pow_dim, len_q, tmp, flag, factor;
  int     *pNum, *pInv;
  
  prim = 2;
  for(i=0; i < np*mp; i++){
    if(prim < pp[i])                              
      prim = pp[i] + 1;           
  }
  pow_dim = 1;
  for(i=0; i < mp; i++)
    pow_dim = pow_dim * prim ;
  
  if ( pow_dim != np ){
#ifdef MATLAB_MEX_FILE
    printf("The given GF(P^M) M-tuple list is not a valid one.\n");
#endif
  } else {
    pow_dim--;
    pNum = Iwork;
    pInv = Iwork+np;
    pNumInv(pp, np, mp, prim, pNum, pInv, pInv+np);
    flag = 0;
    if( pa[len_a-1] != 0 ){
      factor = pow_dim-pa[len_a-1];
      for (i=0; i < len_a; i++){
	if(pa[i] < 0)
	  pa[i] = -Inf;
	else
	  pa[i] =(pa[i] + factor) % pow_dim;
      }
      for (i=0; i < len_b; i++){
	if(pb[i] < 0)
	  pb[i] = -Inf;
	else
	  pb[i] =(pb[i] + factor) % pow_dim;
      }
      flag = 1;
    }
    len_q = len_b - len_a + 1; 
    for(i=0; i < len_q; i++){
      for(j=0; j < len_a-1; j++){
	if(pb[len_b-1-i] < 0 || pa[len_a-2-j] < 0){
	  tmp = -1;
	}else{
	  tmp = (pb[len_b-1-i] + pa[len_a-2-j]) % pow_dim;
	}
	gfplus(&pb[len_b-2-i-j], 1, 1, &tmp, 1, 1, pNum, np, pInv, np, &pb[len_b-2-i-j]);
      }
    }
    /* computes the remainder R */ 
    if (len_a > 1){
      len_r[0] = len_a-1;
      for(i=0; i <len_r[0]; i++)
	pr[i] = pb[i];
      gftrunc(pr, len_r, np*mp, pInv+(mp+1)*np);
      
      if ( flag != 0 ){
	for(i=0; i < len_r[0]; i++)
	  pr[i] = (pr[i] + pow_dim - factor) % pow_dim;
      }
    }else{
      pr[0] = -Inf;
      len_r[0] = 1;
    }
  }
}
/*--end of gfdeconvp()--*/
/*
 * FLXOR  Exclusive OR computation.
 */
static void flxor(px, mx, nx, py, my, ny, pz, mz, nz)  
     int *px, mx, nx, *py, my, ny, *pz, *mz, *nz; 
{
  int     i, incx, incy;
  
  if ((mx == 1) && (nx == 1) && (my == 1) && (ny == 1)){
    *mz = 1;
    *nz = 1;        
    incx = 0;
    incy = 0;
  }else if ((mx == my) && (nx == ny)) {
    *mz = mx;
    *nz = nx;
    incx = 1;
    incy = 1;
  } else if ((my == 1) && (ny == 1)) {
    *mz = mx;
    *nz = nx;
    incx = 1;
    incy = 0;
  } else if ((mx == 1) && (nx == 1)) {
    *mz = my;
    *nz = ny;
    incx = 0;
    incy = 1;
  } else {
#ifdef MATLAB_MEX_FILE
    printf("Matrix dimensions must agree.\n");
#endif
  }
  for (i = 0; i < (*mz)*(*nz); i++) {
    pz[i] =  *px ^ *py;
    px += incx;
    py += incy;
  }
}
/*--end of FLXOR()--*/
/*
 * errlocp1()
 *  Iwork --- 5*(t+2)+(t+4)*(t+1)
 *  Iwork = mu
 *      + t+2 = sigma_mu
 *          + (t+2)*(t+1) = d_mu
 *              + t+2 = l_mu
 *                  + t+2 = mu_l_mu
 *                      + t+2 = indx
 *                          + t+2 = shifted
 *                              + t+1 = tmpRoom
 *                                  + t+1 = bottom of iwork in errlocp1() 
 */
static void errlocp1(syndr, t, pNum, pInv, pow_dim, err_In, Iwork, sigma_Out, len_Out)
     int *syndr, t, *pNum, *pInv, pow_dim, *err_In, *Iwork, *sigma_Out, *len_Out;
{
  int     i, j, k, de_i, de_j, de_k, de_j_tmp;
  int     len_indx, len_tmp, max, rho, shifting, tmp;
  int     *mu, *sigma_mu, *d_mu, *l_mu, *mu_l_mu, *tmpRoom, *indx, *shifted;
  
  /* M-file
   *if type_flag
   * % use the berlekamp's algorithm
   *  t = 2 * t;
   *   mu = [-1, 0:t]';
   *   sigma_mu = [zeros(t+2,1), -ones(t+2, t)];
   *   d_mu = [0; syndrome(1); zeros(t, 1)];
   *   l_mu = [0 0 [1:t]]';