gflib.c

Proakis《contemporarycommunication systems using matlab》matlab源代码

	    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
   *      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