gflib.c

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

	 */
	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]]';
   *   mu_l_mu = mu - l_mu;
   */
  /* allocate (t+2) for *mu */
  mu = Iwork;
  for(i=0; i < t+2; i++)
    mu[i] = i - 1;
  
  /* allocate (t+2)*(t+1) for *sigma_mu */    
  sigma_mu = mu + (t + 2);
  for(i=0; i < (t+2)*(t+1); i++){
    if( i > t+1 )
      sigma_mu[i] = -1;
  }
  
  /* allocate (t+2) for *d_mu */
  d_mu = sigma_mu + (t+1)*(t+2);
  d_mu[1] = syndr[0];
  
  /* allocate (t+2) for *l_mu */
  l_mu = d_mu + (t+2);
  for(i=0; i < t; i++)
    l_mu[i+2] = i + 1;
  
  /* allocate (t+2) for *mu_l_mu */
  mu_l_mu = l_mu + (t+2);
  for(i=0; i < t+2; i++)
    mu_l_mu[i] = mu[i] - l_mu[i];

  /* allocate (t+2) for *indx */
  /* allocate (t+1) for *shifted */
  /* allocate (t+1) for *tmpRoom */
  indx = mu_l_mu + (t+2);
  shifted = indx + (t+2);
  tmpRoom = shifted+(t+1);
  /* M-file
   *% iteratiev start with row three. The first two rows are filled.
   *for de_i = 3:t+2
   */
  for(de_i=2; de_i < t+2; de_i++){
    /*M-file
     *% no more effort to failed situation
     *if (d_mu(de_i - 1) < 0) | err
     *    sigma_mu(de_i, :) = sigma_mu(de_i-1, :);
     *    l_mu(de_i) = l_mu(de_i - 1);
     */
    if( d_mu[de_i-1] < 0 || err_In[0] != 0 ){
      for(j=0; j < t+1; j++)
	sigma_mu[de_i+j*(t+2)] = sigma_mu[de_i-1+j*(t+2)];
      l_mu[de_i] = l_mu[de_i-1];
      /*M-file
       *  else
       *      % find another row proceeding to row de_i -1
       *      % d_mu equals to zero
       *      % and mu - l_mu is the largest.
       *      indx = find(d_mu(1:de_i - 2) >= 0);
       *      rho  = find(mu_l_mu(indx) == max(mu_l_mu(indx)));
       *      rho = indx(rho(1));
       */
    } else {
      len_indx = 0;
      for(j=0; j < de_i-1; j++){
	if(d_mu[j] >= 0){
	  indx[len_indx] = j;
	  len_indx++;
	}
      }
      max = mu_l_mu[indx[0]];
      for(j=0; j < len_indx; j++){
	if(mu_l_mu[indx[j]] >= max)
	  max = mu_l_mu[indx[j]];
      }
      len_tmp = 0;
      for(j=0; j < len_indx; j++){
	if(mu_l_mu[indx[j]] != max){
	  len_tmp++;
	}else{
	  rho = indx[len_tmp];
	  j = len_indx;
	}
      }
      /* M-file
       *    % by (6.25)
       *    % shifted = gfmul(d_mu(de_i - 1), pow_dim - d_mu(rho), tp);
       *    % shifted = gfmul(shifted, sigma_mu(rho, :), tp)';
       *    % multiply inreplace the above two lines.
       *    shifted = -ones(1, t+1);
       *    if (d_mu(de_i - 1) >= 0) & (pow_dim - d_mu(rho) >= 0)
       *        tmp = rem(pow_dim - d_mu(rho) + d_mu(de_i - 1), pow_dim);
       *        indx = find(sigma_mu(rho,:) >= 0);
       *        for de_k = 1 : length(indx)
       *            shifted(indx(de_k)) = rem(tmp + sigma_mu(rho, indx(de_k)), pow_dim);
       *        end;            
       *    end;
       *    % end multiply
       */
      for(j=0; j < t+1; j++)
	shifted[j] = -1;
      
      if( d_mu[de_i-1]>=0 && (pow_dim-d_mu[rho])>=0 ){
	tmp = (pow_dim - d_mu[rho] + d_mu[de_i-1]) % pow_dim;
	/* take advantage of the memory of *indx */
	len_indx = 0;
	for(j=0; j < t+1; j++){
	  if( sigma_mu[rho+j*(t+2)] >= 0 ){
	    indx[len_indx] = j;
	    len_indx++;
	  }
	}
	for(de_k=0; de_k < len_indx; de_k++)
	  shifted[indx[de_k]] = (tmp + sigma_mu[rho+indx[de_k]*(t+2)]) % pow_dim;
      }
      shifting = mu[de_i-1] - mu[rho];
      
      /* M-file
       *  % calculate new sigma_mu
       *  if ~isempty(find(shifted(t-shifting+2 : t+1) >= 0))
       *      % more than t errors, BCH code fails.
       *        err = 1;
       *  else
       *      % calculate the new sigma
       *      shifted = [-ones(1, shifting) shifted(1:t-shifting+1)];
       *      sigma_mu(de_i, :) = gfplus(sigma_mu(de_i-1,:), shifted, tp_num, tp_inv);
       *  end;
       *  l_mu(de_i) = max(l_mu(de_i-1), l_mu(rho) + (de_i - 1) - rho);
       *end;
       */
      /* calculate new sigma_mu */
      len_tmp = 0;
      for(j=t-shifting+1; j < t+1; j++){
	if( shifted[j] >= 0 )
	  len_tmp++;
      }
      
      if (len_tmp != 0){
	err_In[0] = 1;
      }else{
	/* calculate the new sigma */
	for(j=0;j<t+1;j++)
	  tmpRoom[j]=shifted[j];
	
	for(j=0; j < shifting; j++)
	  shifted[j] = -1;
	if(shifting < t+1){
	  for(j=shifting; j < t+1; j++)
	    shifted[j] = tmpRoom[j-shifting];
                }
	for(j=0; j < t+1; j++)
	  gfplus(&sigma_mu[de_i-1+j*(t+2)],1,1, &shifted[j],1,1, pNum,pow_dim+1,pInv,pow_dim+1, &sigma_mu[de_i+j*(t+2)]);
      }
      if( l_mu[de_i-1] < l_mu[rho]+(de_i-1)-rho )
	l_mu[de_i] = l_mu[rho] + (de_i-1) - rho-1;
    }
    /*M-file
     *% calculate d_mu. It is not necessary to do so if mu(de_i) == t
     *if de_i < t+2
     *    % the constant term
     *    d_mu(de_i) = syndrome(mu(de_i) + 1);
     *    indx = find(sigma_mu(de_i, 2:t) >= 0);
     */
    if( de_i < t+1 ){
      d_mu[de_i] = syndr[mu[de_i]];
      len_indx = 0;
      for(j=1; j < t; j++){
	if( sigma_mu[de_i+j*(t+2)] >= 0 ){
	  indx[len_indx] = j-1;
	  len_indx++;
	}
      }
      /*M-file
       *for de_j = 1 : length(indx)
       *    de_j_tmp = indx(de_j);
       *        % Before the "end", it is equivalent to
       *        % d_mu(de_i) = gfadd(d_mu(de_i), ...
       *        %              gfmul(sigma_mu(de_i + 1, de_j_tmp+1), ...
       *        %              syndrome(mu(de_i) * 2 - de_j_tmp + 1), tp), tp);
       *    tmp = syndrome(mu(de_i) - de_j_tmp + 1);
       *    if (tmp < 0) | (sigma_mu(de_i, de_j_tmp + 1) < 0)
       *        tmp = -1;
       *    else
       *        tmp = rem(tmp + sigma_mu(de_i, de_j_tmp + 1), pow_dim);
       *    end;
       *    d_mu(de_i) = gfplus(d_mu(de_i), tmp, tp_num, tp_inv);
       *end;
       *end;--- this 'end' for 'if de_i < t+2
       */
      for(de_j=0; de_j < len_indx; de_j++){
	de_j_tmp = indx[de_j];
	tmp = syndr[ mu[de_i] - de_j_tmp -1];
	if( tmp < 0 || sigma_mu[de_i + (de_j_tmp+1)*(t+2)] < 0 )
	  tmp = -1;
	else
	  tmp = (tmp + sigma_mu[de_i+(de_j_tmp+1)*(t+2)] ) % pow_dim;
	
	gfplus(&d_mu[de_i],1,1,&tmp,1,1,pNum,pow_dim+1,pInv,pow_dim+1, &d_mu[de_i]);
      }
    }
    /*M-file
     *% calculate mu-l_mu
     *mu_l_mu(de_i) = mu(de_i) - l_mu(de_i);
     *end;
     */
    mu_l_mu[de_i] = mu[de_i] - l_mu[de_i];
  }
  /* truncate the reduancy */
  len_Out[0] = 0;
  for(i=0; i < t+1; i++){
    if ( sigma_mu[(t+1) + i*(t+2)] >= 0 )
      len_Out[0] = i;
  }
  len_Out[0] = len_Out[0] + 1;
  for(i=0; i < len_Out[0]; i++)
    sigma_Out[i] = sigma_mu[(t+1)+i*(t+2)];
}
/*--end of errlocp1()--*/
/*
 * errlocp0()
 *  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 = mu2_l_mu
 *                      + t+2 = indx
 *                          + t+2 = shifted
 *                              + t+1 = tmpRoom
 *                                  + t+1 = bottom of iwork in errlocp0() 
 */
static void errlocp0(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, print_flag;
  int     *mu, *sigma_mu, *d_mu, *l_mu, *mu2_l_mu, *tmpRoom, *indx, *shifted;
  
  /*M-file
   *% use simplified algorithm
   *mu = [-1/2, 0:t]';
   *sigma_mu = [zeros(t+2,1), -ones(t+2, t)];
   *d_mu = [0; syndrome(1); zeros(t, 1)];
   *l_mu = [0 0 2*(1:t)]';
   *mu2_l_mu = 2*mu - l_mu;
   */
  /* allocate (t+2) for *mu
   * this *mu is different with M-file, *mu = 2*mu
   */
  mu = Iwork ;
  mu[0] = -1;
  for(i=0; i < t+1; i++)
    mu[i+1] = 2*i;
  
  /* allocate (t+2)*(t+1) for *sigma_mu */    
  sigma_mu = mu + (t+2);
  for(i=0; i < (t+2)*(t+1); i++){
    if( i > t+1 )
      sigma_mu[i] = -1;
    else
      sigma_mu[i] = 0;     
  }
  
  /* allocate (t+2) for *d_mu */
  d_mu = sigma_mu + (t+2)*(t+1);
  d_mu[1] = syndr[0];
  
  /* allocate (t+2) for *l_mu */
  l_mu = d_mu + (t+2);
  for(i=0; i < t; i++)
    l_mu[i+2] = 2*(i + 1);
  
  /* allocate (t+2) for *mu2_l_mu */
  mu2_l_mu = l_mu + (t+2);
  for(i=0; i < t+2; i++)
    mu2_l_mu[i] = mu[i] - l_mu[i];
  
  /* allocate (t+2) for *indx */
  /* allocate (t+1) for *shifted */
  /* allocate (t+1) for *tmpRoom */
  indx = mu2_l_mu + (t+2);
  shifted = indx + (t+2);
  tmpRoom = shifted+(t+1);
  /* M-file
   *% iteratiev start with row three. The first two rows are filled.
   *for de_i = 3:t+2
   */
  for(de_i=2; de_i < t+2; de_i++){
    /*M-file
     *% no more effort to failed situation
     *if (d_mu(de_i - 1) < 0) | err
     *    sigma_mu(de_i, :) = sigma_mu(de_i-1, :);
     */
    if( d_mu[de_i-1] < 0 || err_In[0] != 0 ){
      for(j=0; j < t+1; j++)
	sigma_mu[de_i+j*(t+2)] = sigma_mu[de_i-1+j*(t+2)];
      /*M-file
       *  else
       *      % find another row proceeding to row de_i-1
       *      % d_mu equals to zero
       *      % and 2*mu - l_mu is the largest.
       *      indx = find(d_mu(1:de_i - 2) >= 0);
       *      rho  = find(mu2_l_mu(indx) == max(mu2_l_mu(indx)));
       *      rho = indx(rho(1));
       */
    } else {
      len_indx = 0;
      for(j=0; j < de_i-1; j++){
	if(d_mu[j] >= 0){
	  indx[len_indx] = j;
	  len_indx++;
	}
      }
      max = mu2_l_mu[0];
      for(j=0; j < len_indx; j++){
	if(mu2_l_mu[indx[j]] > max)
	  max = mu2_l_mu[indx[j]];
      }
      for(j=0; j < len_indx; j++){