gflib.c

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

	if(mu2_l_mu[indx[j]] == max)
	  rho = indx[j];
      }
      /* M-file
       *    % by (6.28)
       *    % 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;
      }
      /*M-file
       *shifting = (mu(de_i - 1) - mu(rho))*2;
       */
      shifting = mu[de_i-1] - mu[rho];
      
      /*M-file
       *% calculate new sigma_mu
       *if ~isempty(find(sigma_mu(1:de_i-2,max(1,t-shifting+2):t+1) >=0))
       *    disp('potential error more than posible')
       *    disp(['de_i=',num2str(de_i), ' shifting=',num2str(shifting)]);
       *    disp(sigma_mu)
       *    disp(['shifted=',mat2str(shifted)]);
       *end;
       */
      
      if ( t-shifting+2 > 1 )
	max = t - shifting + 1;
      else
	max = 0;
      print_flag = 0;
      for(j=0; j < de_i-2; j++){
	for(k=max; k < t+1; k++){
	  if( sigma_mu[j+k*(t+2)] >= 0 )
	    print_flag++;
	}
      }
      if ( print_flag != 0 ){
	/*                printf("potential error more than posible\n");
			  printf("de_i=%d,shifting=%d\n",de_i,shifting);
			  for(j=0; j<(t+2)*(t+1); j++)
			  printf("sigma_mu=%d\n",sigma_mu[j]);
			  printf("shifted=");
			  for(j=0; j < t+1; j++)
			  printf(" %d\n",shifted[j]); */
      }
      /*M-file
       *if ~isempty(find(shifted(max(t-shifting+2, 1) : 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;
       *end;--this end for "if (d_mu(de_i - 1) < 0) | err"
       *l_mu(de_i) = max(find(sigma_mu(de_i,:) >= 0)) - 1;
       */
      /* calculate new sigma_mu */
      len_tmp = 0;
      for(j=max; j < t+1; j++){
	if( shifted[j] >= 0 )
	  len_tmp++;
      }
      if (len_tmp > 0){
	err_In[0] = 1;
      }else{
	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)]);
      }
    }
    for(j=0; j < t+1; j++){
      if(sigma_mu[de_i+j*(t+2)] >= 0)
	l_mu[de_i] = j;
    }
    /*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)*2 + 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++;
	}
      }
      
      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]);
      }
    }
    mu2_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+1;
  }
  for(i=0; i < len_Out[0]; i++)
    sigma_Out[i] = sigma_mu[(t+1)+i*(t+2)];
}
/*-- end of errlocp0()--*/
/*
 *bchcore()
 *  Integer Working Space list:
 *  total for function: Iwork = (2+dim)*(2*pow_dim+1)+t*t+15*t+16;
 *  Iwork = pNum
 *      + pow_dim+1 = pInv
 *          + pow_dim+1 = Iwork for pNumInv()
 *              + (pow_dim+1)*dim = non_zeros_itm
 *                  + pow_dim = syndrome
 *                      + 2*t = tmpRoom
 *                          + 2*t = sigma 
 *                              + (t+1) = len_sigma
 *                                  + 1 = Iwork for errlocp0()
 *                                      + 5*(t+2)+(t+4)*(t+1) = loc_err
 *                                          + pow_dim = pos_err
 *                                              + pow_dim*dim = bottom of Iwork
 */
static void bchcore(code, pow_dim, dim, k, t, pp, Iwork, err, ccode)
     int *code, pow_dim, dim, k, t, *pp, *Iwork, *err, *ccode;
{
  int     i, j, prim, np, mp, nk, tmp, er_j;
  int     *loc_err, num_err, cnt_err, len_pos_err, er_i, test_flag;
  int     *non_zeros_itm, len_non_z_itm, *tmpRoom, *len_sigma, *Iwork1;
  int     *pNum, *pInv, *sigma, *pos_err, *syndrome;
  
  np = pow_dim + 1;
  mp = dim;
  /*  code = code(:)';
   *  err = 0;
   *  tp_num = tp * 2.^[0:dim-1]';
   *  tp_inv(tp_num+1) = 0:pow_dim;
   */
  err[0] = 0; 
  pNum = Iwork;
  pInv = Iwork + np;
  prim = 2;
  pNumInv(pp, np, mp, prim, pNum, pInv, pInv+np);
  
  /*  % **(1)** syndrome computation.
   *  % initialization, find all non-zeros to do the calculation
   *  non_zero_itm = find(code > 0) - 1;
   *  len_non_z_itm = length(non_zero_itm);
   *  syndrome = -ones(1, 2*t);
   */
  /* allocate pow_dim for *non_zeros_itm */
  non_zeros_itm = pInv + (dim+1)*np;    
  for(i=0; i < pow_dim; i++)
    non_zeros_itm[i] = 0;
  len_non_z_itm = 0;    
  for(i=0; i < pow_dim; i++){
    if(code[i] > 0 ){
      non_zeros_itm[len_non_z_itm] = i;
      len_non_z_itm++;
    }
  }
  /* allocate 2*t for *syndrome */
  syndrome = non_zeros_itm + pow_dim;
  for(i=0; i < 2*t; i++)
    syndrome[i] = -1;
  
  /*  % syndrome number is 2*t where t is error correction capability
   *  if len_non_z_itm > 0
   *      tmp = 1:2*t;
   *      syndrome(tmp) = non_zero_itm(1) * tmp;
   *      if len_non_z_itm > 1
   *          for n_k = 2 : len_non_z_itm
   *              syndrome(tmp) = gfplus(syndrome(tmp), non_zero_itm(n_k) * tmp, tp_num, tp_inv);
   *          end;
   *      end;
   *  end;
   *  % complete syndrome computation
   */ 
  tmpRoom = syndrome + 2*t; /* size is 2*t */
  if( len_non_z_itm > 0 ){
    for(i=0; i < 2*t; i++)
      syndrome[i] = non_zeros_itm[0]*(i+1); 
    if( len_non_z_itm > 1 ){
      for( nk=1; nk < len_non_z_itm; nk++ ){
	for(i=0; i < 2*t; i++)
	  tmpRoom[i] = non_zeros_itm[nk]*(i+1);
	gfplus(syndrome,2*t,1,tmpRoom,2*t,1,pNum,np,pInv,np,syndrome);
      }
    }
  }
  /* complete syndrome computation */  
  /*  % **(2)** determine the error-location polynomial.
   *  % This step is the most complicated part in the BCH decode.
   *  % reference to p158 of Shu Lin's Error Control Coding,
   *  % the simplified algorithm for finding sigma_x.
   *  % the maximum degree of the error-location polynomial is t
   *  % if the degree is larger than t, there are more than t errors and
   *  % the algorithm cannot find a solution to correct the errors.
   *
   *  % for BCH code, use simplified method. Note that when you call
   *  % errlocp, the parameter should be exact.
   *  [sigma, err] = errlocp(syndrome, t, tp, pow_dim, err, 0);
   */
  /* allocate (t+1) for sigma*/
  sigma = tmpRoom+2*t;
  for(i=0; i<t+1;i++)
    sigma[i] = 0;
  len_sigma = sigma+(t+1);
  errlocp0(syndrome, t, pNum, pInv, pow_dim, err, len_sigma+1, sigma, len_sigma);
  /* need 5*(t+2)+(t+4)*(t+1) for errlocp0 */
  
  /*  % ***(3)*** computation of error-location numbers.
   *  loc_err = zeros(1, pow_dim);
   */
  /* allocate pow_dim for loc_err
   * it is automatically initialized as zero.
   */ 
  loc_err = len_sigma+1+5*(t+2)+(t+4)*(t+1);;
  for(i=0; i<pow_dim; i++)
    loc_err[i] = 0;
  /*  % in case of failed or no error, skip.
   *  num_err = length(sigma) - 1;
   */
  num_err = len_sigma[0] - 1;
  /*  if (~err) & (num_err > 0)
   *      cnt_err = 0;
   *      pos_err = [];
   *      er_i = 0;
   */
  /* allocating pow_dim*dim for pos_err */
  pos_err = loc_err + pow_dim;
  if( err[0] == 0 && num_err > 0 ){
    cnt_err = 0;
    er_i = 0;
    /*      while (cnt_err < num_err) & (er_i < pow_dim * dim)
     *          test_flag = sigma(1);
     *          for er_j = 1 : num_err
     */
    while( cnt_err < num_err && er_i < pow_dim* dim ){
      test_flag = sigma[0];
      for( er_j=0; er_j < num_err; er_j++){
	/*              if sigma(er_j + 1) >= 0
	 *                  % The following 6 lines is equivelent to
	 *                  % tmp = gfmul(er_i * er_j, sigma(er_j+1), tp);
	 *                  tmp = er_i * er_j;
	 *                  if (tmp < 0) | (sigma(er_j+1) < 0)
	 *                      tmp = -1;
	 *                  else
	 *                      tmp = rem(tmp + sigma(er_j + 1), pow_dim);
	 *                  end;
	 *                  test_flag = gfplus(test_flag, tmp, tp_num, tp_inv);
	 *              end;
	 *          end;
	 */
	if(sigma[er_j + 1] >= 0){
	  tmp = er_i * (er_j+1);
	  if( tmp < 0 || sigma[er_j+1] < 0 )
	    tmp = -1;
	  else
	    tmp = (tmp + sigma[er_j+1]) % pow_dim;
	  
	  gfplus(&test_flag,1,1,&tmp,1,1,pNum,np,pInv,np,&test_flag);
	}
      }
      /*          if test_flag < 0
       *              cnt_err = cnt_err + 1;
       *              pos_err = [pos_err, rem(pow_dim-er_i, pow_dim)];
       *          end;
       *          er_i = er_i + 1;
       *      end;
       */
      if ( test_flag < 0 ){
	pos_err[cnt_err] = (pow_dim - er_i) % pow_dim;
	cnt_err++;
      }
      er_i++;
    }
    /*      pos_err = rem(pow_dim+pos_err, pow_dim);
     *      pos_err = pos_err + 1; % shift one location because power zero is one.
     *      loc_err(pos_err) = ones(1, cnt_err);
     *      err = num_err;
     */
    for(i=0; i < cnt_err; i++)
      pos_err[i] = (pow_dim + pos_err[i]) % pow_dim;
    for(i=0; i < cnt_err; i++)
      loc_err[pos_err[i]] = 1;
    err[0] = num_err;
  }else{
    /*  else
     *      if err
     *          err = -1;
     *      end;
     *  end;
     *  % completed error location detection
     */
    if( err[0] != 0 )
      err[0] = -1;
  }
  /* completed error location detection */
  /*  % correct the error
   *  ccode = rem(code + loc_err, 2);
   *  msg = ccode(pow_dim-k+1 : pow_dim);
   *  %-- end of bchcore---
   */
  for(i=0; i < pow_dim; i++)
    ccode[i] = (code[i] + loc_err[i]) % 2;
}
/*--end of bchcore()--*/
/*
 * rscore()
 * Integer Working Space list:
 *  total for function: Iwork = (pow_dim-2*k+5*dim+18)*pow_dim+3*dim+k*(k-13)+27;
 *  Iwork = pNum
 *      + pow_dim+1 = pInv
 *          + pow_dim+1 = Iwork for pNumInv()
 *              + dim*(pow_dim+1)= tmpRoom
 *                  + pow_dim = syndrome
 *                      + pow_dim-k = sy
 *                          + 2 = Iworkgfdec
 *                              + (2+dim)*pow_dim+3+dim = sigma
 *                                  + pow_dim-k+1 = len_sigma
 *                                      + 1 = Iworkerr
 *                                          + (pow_dim-k)^2+10*(pow_dim-k)+14 = loc_err
 *                                              + pow_dim = pos_err
 *                                                  + pow_dim*dim = Z
 *                                                      + pow_dim-k+1 = IworkMul
 *                                                          + 3+pow_dim+dim = er_loc
 *                                                              + pow_dim*dim = bottom of Iwork
 */
static void rscore(code, k, pp, dim, pow_dim, Iwork, err, ccode)
     int *code, k, *pp, dim, pow_dim, *Iwork, *err,*ccode;
{
  int     i, j, np, mp, sy_i, er_i, er_j, am_i, am_j, syn;
  int     t, t2, prim, tmp, test_flag, len_sy, len, zero;
  int     num, den, num_err, cnt_err, pos_err_inv;
  int     *pNum, *pInv, *tmpRoom, *syndrome, *sy, *sigma, *loc_err, *pos_err, *Z, *er_loc;
  int     *Iworkgfdec, *Iworkerr, *IworkMul, *len_sigma, len_syndrome[1];
  
  np = pow_dim + 1;
  mp = dim;
  t2 = pow_dim - k;
  t =(int)(t2/2);
  
  prim = 2;
  /* allocate 2*np for both of *pNum and *pInv */
  pNum = Iwork;
  pInv = Iwork + np;
  pNumInv(pp, np, mp, prim, pNum, pInv, pInv+np);
  /*
   * allocate pow_dim for *tmpRoom, t2 for *syndrome and 2 for *sy; 
   */
  tmpRoom = pInv + np +np*dim;
  syndrome = tmpRoom + pow_dim ;
  sy = syndrome + pow_dim - k;
  Iworkgfdec = sy + 2;
  len_syndrome[0] = 1;
  for( sy_i=0; sy_i < t2; sy_i++){
    sy[0] = sy_i+1;
    sy[1]=0;
    len_sy = 2;
    for(i=0; i < pow_dim; i++)
      tmpRoom[i] = code[i];
    gfdeconvp(tmpRoom, pow_dim, sy, len_sy, pp, np, mp, &syndrome[sy_i], len_syndrome, Iworkgfdec);
  }
  /*% (2) find error location polynomial
   * allocate pow_dim-k+1 for *sigma, and (pow_dim-k)*(pow_dim-k)+10*(pow_dim-k)+14 for errlocp1(); 
   */
  sigma = Iworkgfdec+1+(dim+2)*(pow_dim+1);
  for(i=0; i<pow_dim-k+1; i++)
    sigma[i] = 0;
  len_sigma = sigma+pow_dim-k+1;
  Iworkerr = len_sigma + 1;
  err[0] = 0;
  errlocp1(syndrome, 2*t, pNum, pInv, pow_dim, err, Iworkerr, sigma, len_sigma);
  
  /*% (3) find solution for the polynomial
   * allocate pow_dim for *loc_err
   */
  loc_err = Iworkerr+(pow_dim-k)*(pow_dim-k)+10*(pow_dim-k)+14;
  for(i=0; i<pow_dim; i++)
    loc_err[i] = -Inf;
  num_err = len_sigma[0] - 1;
  
  if( num_err > t)
    err[0] = 1;
  
  pos_err = loc_err + pow_dim;
  if( err[0] == 0 && num_err > 0 ){
    cnt_err = 0;
    /* allocating pow_dim*dim for pos_err */
    er_i = 0;
    while( cnt_err < num_err && er_i < pow_dim*dim ){
      test_flag = sigma[0];
      for( er_j=0; er_j < num_err; er_j++){
	if(sigma[er_j + 1] >= 0){
	  tmp = er_i * (er_j+1);
	  if( tmp < 0 || sigma[er_j+1]<0 )
	    tmp = -1;
	  else
	    tmp = (tmp + sigma[er_j+1]) % pow_dim;
	  
	  gfplus(&test_flag,1,1,&tmp,1,1,pNum,np,pInv,np,&test_flag);
	}
      }
      if ( test_flag < 0 ){
	pos_err[cnt_err] = (pow_dim - er_i) % pow_dim;
	cnt_err++;
      }
      er_i++;
    }
    for(i=0; i < cnt_err; i++)
      pos_err[i] = (pow_dim + pos_err[i]) % pow_dim;
    err[0] = num_err;
  }else{
    if( err[0] != 0 )
      err[0] = -1;
  }
  
  /* allocate 2*t+1 */
  /* allocate 2+np+mp for IworkMul of gfpmul() */
  /* allocate Pow_dim*dim for *er_loc */
  Z = pos_err+ dim*pow_dim;
  IworkMul = Z + pow_dim - k + 1;
  er_loc = IworkMul + 2+np+mp;
  if( err[0] > 0 ){
    Z[0] = 1;
    for( am_i=1; am_i<=num_err; am_i++){
      gfplus(&syndrome[am_i-1],1,1,&sigma[am_i],1,1,pNum,np,pInv,np,&Z[am_i]);
      if( am_i > 1 ){
	len = 1;
	for(am_j=1; am_j <= am_i-1; am_j++){
	  gfpmul(&sigma[am_j],1,&syndrome[am_i-am_j-1],1,pp,np,mp,&tmpRoom[0],&len,IworkMul);
	  gfplus(&Z[am_i],1,1,&tmpRoom[0],1,1,pNum,np,pInv,np,&Z[am_i]);
	}
      }
    }
    for(am_i=0; am_i < cnt_err; am_i++){
      num = 0;
      den = 0;
      pos_err_inv = (pow_dim - pos_err[am_i]) % pow_dim;
      for(am_j=1; am_j <= num_err; am_j++){
	tmp = pos_err_inv * am_j;
	
	if( tmp < 0 || Z[am_j] < 0 )
	  tmp = -1;
	else
	  tmp = (tmp + Z[am_j]) % pow_dim;    
	
	gfplus(&num,1,1,&tmp,1,1,pNum,np,pInv,np,&num);
	
	if ( am_i != am_j-1 ){
	  if ( den < 0 ){
	    den = -1;
	  }else{
	    if( (pos_err[am_j-1]) < 0 || pos_err_inv < 0 )
	      tmp = -1;
	    else
	      tmp = ( pos_err[am_j-1]+pos_err_inv ) % pow_dim;
	    
	    zero = 0;
	    gfplus(&zero,1,1,&tmp,1,1,pNum,np,pInv,np, &tmp);
	    
	    if ( tmp < 0 )
	      den = -1;
	    else
	      den = ( den + tmp ) % pow_dim;
	  }
	}
      }
      tmp = pow_dim - den;
      if( tmp < 0 || num < 0 )
	er_loc[am_i] = -1;
      else
	er_loc[am_i] = ( tmp + num ) % pow_dim;
    } 
    for(i=0; i < cnt_err; i++)
      loc_err[ pos_err[i] ] = er_loc[i];
  }
  /*calculate *ccode
   *and don't have to calculate *msg since it just is the part of *ccode.
   */
  gfplus(loc_err,pow_dim,1,code,pow_dim,1,pNum,np,pInv,np,ccode);
}