gflib.c

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

   *      b = [b zeros(1, na-nb)];
   *    elseif nb > na
   *      a = [a zeros(1, nb-na)];
   *    end;
   *    c = rem(a+b, p);
   */
  if (len_a > len_b){
    for (i=len_b; i <len_a; i++)
      pc[i] = pa[i] % pp;
    for (i=0; i < len_b; i++)
      pc[i] = (pa[i]+pb[i]) % pp;
  } else if (len_b > len_a){
    for (i=len_a; i <len_b; i++)
      pc[i] = pb[i] % pp;
    for (i=0; i < len_a; i++)
      pc[i] = (pa[i]+pb[i]) % pp;       
  } else {
    for(i=0; i < len_a;  i++)
      pc[i] = (pa[i]+pb[i]) % pp;
  }
  /*    % truncate the result as requested.
   *    if (nargin > 3)
   *        nc = max(na, nb);
   *        if isempty(len)
   *           c = gftrunc(c);
   *        elseif len < 0
   *           c = gftrunc(c);
   *        elseif len <= nc
   *            c = c(1:len);
   *        else
   *            c = [c zeros(1, len-nc)];
   *        end;
   *    end;
   */
  /* truncate the result as requested.*/
  if(len <= 0){
    gftrunc(pc,nc,1,Iwork);
  } else if(len > 0){
    nc[0] = len;
  }
}
/*--end of GFADD()--*/
/*
 * Multiply for Matrix P
 *      Iwork --- 2*nc+nc*mp+np
 *      Iwork = indx
 *          + *nc = sum_ab
 *              + *nc = tmp
 *                  +(*nc)*mp = sum
 *                      + np = bottom of Iwork
 */
static void gfpmul(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     i, j, k, prim, len_indx;
  int     *sum_ab, *indx, *tmp, *sum;
  
  indx = Iwork;
  sum_ab = indx + nc[0];
  tmp = sum_ab + nc[0];
  sum = tmp + nc[0]*mp;
  /* prim = max(max(p))+1; */
  prim = 0;
  for(i=0; i < np*mp; i++){
    if(prim < pp[i])                              
      prim = pp[i] + 1;           
  }
  /* fingure out input elements less than zero and chang them to -Inf */
  /* indx = [find(sum_ab < 0)  find(isnan(sum_ab))];*/
  for(i=0; i < len_a; i++){
    if(pa[i] < 0)
      pa[i] = -Inf;
  }
  for(i=0; i < len_b; i++){
    if(pb[i] < 0)
      pb[i] = -Inf;
  }
  if (len_a == len_b || len_a == 1 || len_b == 1){
    if ( len_a == len_b ){
      for (i = 0; i < nc[0]; i++){
	if (pa[i] == -Inf || pb[i] == -Inf)
	  sum_ab[i] = -1;
	else
	  sum_ab[i] =(pa[i] + pb[i]) % (np - 1);
      }
    }else{
      if(len_a == 1 && len_b != 1){
	for (i = 0; i < len_b; i++){
	  if (pa[0] < 0 || pb[i] < 0)
	    sum_ab[i] = -1;
	  else
	    sum_ab[i] = (pa[0] + pb[i]) % (np - 1);
	}
      }
      if (len_b == 1 && len_a != 1){
	for (i = 0; i < len_a; i++){
	  if (pb[0] < 0 || pa[i] < 0)
	    sum_ab[i] = -1;
	  else
	    sum_ab[i] = (pa[i] + pb[0]) % (np - 1);
	}
      }
    }
  } else {
    if (len_a <= len_b)
      nc[0] = len_a;
    else
      nc[0] = len_b;
    
    for (i=0; i < nc[0]; i++){
      if (pa[i] < 0 || pb[i] < 0)
	sum_ab[i] = -1;
      else
	sum_ab[i] =(pa[i] + pb[i]) % (np - 1);
    }
  }
  /*tmp = p(sum_ab +2, :); */
  for(i=0; i < nc[0]; i++){
    for(j=0; j < mp; j++)
      tmp[i+j*nc[0]] = pp[j*np+sum_ab[i]+1];
  }
  for(i=0; i < nc[0]; 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
    }
  }
  len_indx = 0;
  for(i=0; i < nc[0]; i++){
    if( pc[i] < 0 ){
      indx[len_indx] = i;
      len_indx++;
    }
  }
  if(len_indx != 0){
    for(i=0; i < len_indx; i++)
      pc[indx[i]] = -Inf;
  }
}
/*--end of GFpMUL()---*/
/*
 * multiply for scalar p
 */
static void gfmul(pa, len_a, pb, len_b, pp, pc)
     int *pa, len_a, *pb, len_b, pp, *pc;
{
  int     i;

  /*    len_a = length(a);
   *    len_b = length(b);
   *    if (len_a == len_b) | (len_b == 1)
   *        c = rem(a .* b, p);
   *    elseif len_a == 1
   *        c = rem(b .* a, p);   
   *    elseif len_a > len_b
   *        c = a;
   *        c(1:len_b) = rem(a(1:len_b) .* b, p);
   *    else
   *        c = b;
   *        c(1:len_a) = rem(b(1:len_a) .* a, p);
   *    end;
   */
  if (len_a == len_b) {
    for (i=0; i < len_a; i++)
      pc[i] = (pa[i]*pb[i]) % pp;
  }else if(len_b == 1){
    for (i=0; i < len_a; i++)
      pc[i] = (pa[i]*pb[0]) % pp;
  } else if (len_a == 1){
    for (i=0; i < len_b; i++)
      pc[i] = (pb[i]*pa[0]) % pp;
  } else if (len_a > len_b){ 
    for (i=len_b; i < len_a; i++)
		    pc[i] = pa[i];
	    for (i=0; i < len_b; i++)
	        pc[i] = (pa[i]*pb[i]) % pp;
    } else {
        for (i=len_a; i < len_b; i++)
		    pc[i] = pb[i] ;
	    for (i=0; i < len_a; i++)
	        pc[i] = (pb[i]*pa[i]) % pp;
   	}
}   
/*--end of GFMUL()--*/
/*
 * GFconv() when P is a scalar.
 *      Iwork is needed to be assigned for nc.
 *      Iwork --- 6*(len_b+len_a)-4
 *      Iwork = One
 *          + 1 = pfb / pfa
 *              + len_a+len_b = pfa/pfb
 *                  + len_b+len_a-1 = Iwork for gffilter()
 *                      + 3*len_b+3*len_a-3 = Iwork for fliplr(pc)
 *                          + len_b+len_a-1 = bottom of Iwork in gfconv()
 */
static void gfconv(pa, len_a, pb, len_b, pp, pc, Iwork)
    int *pa, len_a, *pb, len_b, pp, *pc, *Iwork;
{
  int     i, nc;
  int     *pfa, *pfb, *One;
  
  /*    %variable assignment
   *    a = fliplr(a);
   *    b = fliplr(b);
   */
  One = Iwork;
  One[0] = 1; 
  nc = len_a+len_b-1;
  /* convolution computation */
  /*    if na > nb
   *        if nb > 1
   *            a(na+nb-1) = 0;
   *        end
   *        c = gffilter(b, 1, a, p);
   *    else
   */
  if(len_a > len_b){
    pfb = Iwork+1;
    for(i=0; i<len_b; i++)
      pfb[i] = pb[len_b-1-i];
    pfa = pfb+len_b+len_a;
    for(i=0; i<nc; i++){
      if(i<len_a)
	pfa[i] = pa[len_a-1-i];
      else
	pfa[i] = 0;
    }
    /* c = gffilter(b,1,a,p);*/
    gffilter(pfb, len_b, One, 1, pfa, nc, pp, pc, pfa+nc);
        /* c = fliplr(c) */
    fliplr(pc,nc,pfa+3*nc+len_b-1);
  }else{
    /*        if na > 1
     *            b(na+nb-1) = 0;
     *        end
     *        c = gffilter(a, 1, b, p);
     *    end
     */
    pfa = Iwork+1;
    for(i=0; i<len_a; i++)
      pfa[i] = pa[len_a-1-i];
    pfb = pfa+len_a+len_b;
    for(i=0; i < nc; i++){
      if(i<len_b)
	pfb[i] = pb[len_b-1-i];
      else
	pfb[i] = 0;
    }
    /* c = gffilter(a,1,b,p);*/
    gffilter(pfa, len_a, One, 1, pfb, nc, pp, pc, pfb+nc);
    /* c = fliplr(c) */
    fliplr(pc,nc,pfb+3*nc+len_a-1);
  }
}
/*--end of GFCONV()--*/
/*
 * GFpconv() when this computation is in GF(p^m) field.
 *  where the Iwork need to be assigned for 5+3*(np+mp).
 *  Iwork --- 5+3*(np+mp)
 *  Iwork = Iwork for gfpmul()
 *      + 2+mp+np = Iwork for gfpadd()
 *          + 1+mp+np = Iwork for gfpmul()
 *              + 2+mp+np = bottom of Iwork in gfpconv()
 */
static void gfpconv(pa, len_a, pb, len_b, pp, np, mp, pc, Iwork)
     int *pa, len_a, *pb, len_b, *pp, np, mp, *pc, *Iwork;
{
  int     i, j, k, prim, pow_dim;
  int     *tmp, One;
  
  /* p = max(max(tp))+1 */
  prim = 2;
  for(i=0; i < np*mp; i++){
    if(prim < pp[i])                              
            prim = pp[i] + 1;           
  }
  /* check up if the given GF(P^M) M-tuple list is a valid one */
  /*    [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;
   *    c = [0 : na + nb] - Inf;
   */
  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.");
#endif
  for (i=0; i < len_a+len_b+1; i++)
    pc[i] = -Inf;
  /*    for i = 1 : na
   *        for j = 1 : nb
   *            k = i+j-1;
   *            if c(k) >= 0
   *                c(i+j-1) = gfadd(c(i+j-1), gfmul(a(i), b(j), tp), tp);
   *            else
   *                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;