www.pudn.com > communicationmatlab.rar > GFLIB.C
/* ============================================================================
* gflib.c is a library file including the following functions:
* P.S.: All the parameters in these functions have the type of integer.
* static void gftrunc(*pa, *len_a, len_p, *Iwork)
* M_FILE: C = GFTRUNC(A); Iwork --- len_a
* truncates the results of every GF's computation .
* static void gfpadd(*pa, len_a, *pb, len_b, *pp, np, mp, *pc, *nc, *Iwork)
* M_FILE: C = GFADD(A,B,P); Iwork --- *nc+*nc*mp+np
* computes A adds B in GF when P is a matrix.
* static void gfadd(*pa, len_a, *pb, len_b, pp, len, *pc, *nc, *Iwork)
* M_FILE: C = GFADD(A,B,P,LEN), Iwork --- *nc
* computes A adds B in GF when P is a scalar.
* static void gfpmul(*pa, len_a, *pb, len_b, *pp, np, mp, *pc, *nc, *Iwork)
* M_FILE: C = GFMUL(A,B,P); Iwork --- 2*nc+nc*mp+np
* computes A multiply B in GF(P) when P is a matrix.
* static void gfmul(*pa, len_a, *pb, len_b, pp, *pc)
* M_FILE: C = GFMUL(A,B,P); no Iwork
* computes A multiply B in GF(P) when P is a scalar.
* static void gfconv(*pa, len_a, *pb, len_b, pp, *pc, *Iwork)
* M_FILE: C = GFCONV(A,B,P); Iwork --- 6*(len_b+len_a)-5
* computes the convolution between two GF(P) polynomials when P is a scalar.
* static void gfpconv(*pa, len_a, *pb, len_b, *pp, np, mp, *pc, *Iwork)
* M_FILE: C = GFCONV(A,B,P); Iwork --- 5+3*(np+mp)
* computes the convolution between two GF(P) when P is a Matrix.
* static void gfdeconv(*pb, len_b, *pa, len_a, pp, *pq, len_q, *pr, *len_r, *Iwork)
* M_FILE: [Q, R] = GFDECONV(B,A,P); Iwork --- 14*len_b-2*len_a+2
* computes the deconvolution between two GF(P) polynomials when P is a scalar.
* static void gfpdeconv(*pb, len_b, *pa, len_a, *pp, np, mp, *pq, *pr, *len_r, *Iwork)
* M_FILE: [Q, R] = GFDECONV(B,A,P); Iwork --- (2+mp)*np+len_a-1
* computes the deconvolution between two GF(P) when P is a Matrix.
* static void gfplus(*pi, mi, ni, *pj, mj, nj, *alpha, len_alpha, *beta, len_beta, *pk)
* M_FILE: C = GFPLUS(A,B,ALPHA,BETA)
* Galois Field addition (vectorized).
* static void gffilter(*pb, len_b, *pa, len_a, *px, len_x, p, *pOut, *Iwork)
* M_FILE: Y = GFFILTER(B,A,X,P); Iwork --- 2*len_x+len_b+len_a-2
* filters the data in GF(P).
* static void flxor(*px, mx, nx, *py, my, ny, *pz, *mz, *nz)
* M_FILE: Z = FLXOR(X,Y)
* exclusive OR computation.
* static void fliplr(*pa, len_a, *Iwork)
* M_FILE: Y = FLIPLR(X)
* Iwork --- len_a
* flips matrix in the left/right direction.
* static int isprime(p)
* M_FILE: ISP = ISPRIME(X)
* return TRUE for prime numbers.
* static void errlocp1(*syndr,t,*pNum,*pInv,pow_dim,err_In,*Iwork,*sigma_Out,*len_Out)
* M_FILE: [SIGMA, ERR] = ERRLOCP(SYNDROME,T,TP,POW_DIM,ERR,1);
* Iwork --- 5*(2+t)+(t+4)*(t+1) PS: t=2*T (T is in M_file)
* static void errlocp0(*syndr,t,*pNum,*pInv,pow_dim,err_In,*Iwork,*sigma_Out,*len_Out)
* M_FILE: [SIGMA, ERR] = ERRLOCP(SYNDROME,T,TP,POW_DIM,ERR,0);
* Iwork --- 5*(2+t)+(t+4)*(t+1) PS: t=T (T is in M_file)
* static void bchcore(*code,pow_dim,dim,k,t,*pp,*Iwork,*msg,*err,*ccode)
* M_FILE: [SIGMA, ERR, CCODE] = BCHCORE(CODE,POW_DIM,DIM,K,T,P);
* Iwork --- (2+dim)*(2*pow_dim+1)+t*t+15*t+16
* static void rscore(*code,k,*pp,dim,pow_dim,*Iwork,*msg,*err,*ccode)
* M_FILE: [SIGMA, ERR, CCODE] = RSCORE(CODE,K,P,DIM,POW_DIM);
* Iwork --- (pow_dim-2*k+5*dim+18)*pow_dim+3*dim+k*(k-13)+27
* static void pNumInv(*pp, np, mp, prim, *pNum, *pInv, *Iwork)
* M_CODE: tp_num = tp * prim.^[0 : tp_m-1]'; -- vector pNum
* tp_inv(tp_num+1) = 0:pow_dim; -- vector pInv
* Iwork --- np*mp
* static void bi2de(*pbi, np, mp, prim, *pde)
* ============================================================================
*
* Original designed by Wes Wang,
* Jun Wu, The Mathworks, Inc.
* Jan-25, 1996
*
* Copyright (c) 1996 by The MAthWorks, Inc.
* All Rights Reserved
*
* ===========================================================================
*/
#define Inf 32766
/*
*bi2de()
* Iwork --- np*mp
*/
static void bi2de(pbi, np, mp, prim, pde)
int *pbi, np, mp, prim, *pde;
{
int i, j, k;
for(i=0; i 0){
for (k=0; k < j; k++)
pbi[i+j*np] = pbi[i+j*np]*prim;
}
}
pde[i] = pde[i] + pbi[i+j*np];
}
}
}
/*end of bi2de()*/
/*
*fliplr()
* Iwork --- len_a
*/
static void fliplr(pa, len_a, Iwork)
int *pa, len_a, *Iwork;
{
int i;
int *tmp;
tmp = Iwork;
for(i=0;i < len_a; i++)
tmp[i]= pa[len_a-1-i];
for(i=0;i 0){
for (k=0; k < j; k++)
pp[i+j*np] = pp[i+j*np]*prim;
}
}
pNum[i] = pNum[i] + pp[i+j*np];
}
pInv[pNum[i]] = i;
}
}
/*--end of pNumInv()--*/
/*
* Check if P is a prime number.
* Return 1 when p is a prime number and return 0 when it is not.
*/
static int isprime(p)
int p;
{
int i, Fox, lenFox;
lenFox = 0;
for ( i=2; i < p; i++){
Fox = p % i;
if(Fox == 0)
lenFox++;
}
if(lenFox == 0)
return(1);
else
return(0);
}
/*--end of ISPRIME()--*/
/*
* GFFilter()
* Iwork --- 2*len_x+len_b+len_a-2
* Iwork = tmp_x
* + len_b+len_x-1 = tmp_y
* + len_a+len_x-1 = bottom of Iwork in gffilter()
*/
static void gffilter(pb, len_b, pa, len_a, px, len_x, p, pOut, Iwork)
int *pb, len_b, *pa, len_a, *px, len_x, p, *pOut, *Iwork;
{
int i, j, q, indx;
int tmp_b, tmp_a, len_xx, len_yy;
int *tmp_x, *tmp_y;
/* make a(1) to be one. */
if (pa[0] == 0){
#ifdef MATLAB_MEX_FILE
printf("First denominator filter coeficient must be non-zero!");
#endif
} else if (pa[0] != 1){
if( isprime(p) ){ /* by Fermat's theory (pp 96 in MacWilliams)
* for any integer q less than p.
*/
q = 1;
for(i=0; i < p-2; i++){
q = q*pa[0];
}
if (((q*pa[0]) % p) != 1){
#ifdef MATLAB_MEX_FILE
printf("Warning: The computation may be not accurate because a large integer\n");
printf(" is introduced in to the calculation.\n");
#endif
}
} else {
#ifdef MATLAB_MEX_FILE
printf("P must be a prime in GFFILTER.");
#endif
}
/* made pa[0] = 1. */
for(i=0; i < len_a; i++)
pa[i] = (pa[i]*q) % p;
}
/* computation length adjustment */
len_xx = len_b+len_x-1;
len_yy = len_a+len_x-1;
tmp_x = Iwork;
tmp_y = Iwork + len_xx;
for(i=len_b-1; i < len_xx; i++)
tmp_x[i] = px[i-len_b+1];
/*
* flip the computation variable
* difference equation iteration
*/
len_a--;
for (i=0; i < len_x; i++){
indx = i + len_a;
tmp_a = 0;
tmp_b = 0;
if ( len_a > 0 ){
/*y(indx) = fb * xx(i:len_b +i-1)' + fa * y(i:indx-1)'; */
for (j=0; j < len_b; j++)
tmp_b = tmp_b + pb[len_b-1-j]*tmp_x[i+j];
for (j=0; j < len_a; j++)
tmp_a = tmp_a + (p-pa[len_a-j])*tmp_y[i+j];
tmp_y[indx] = tmp_b + tmp_a;
}else{
/* y(indx) = fb * xx(i:len_b+i-1)'; */
for (j=0; j < len_b; j++)
tmp_b = tmp_b + pb[len_b-1-j]*tmp_x[i+j];
tmp_y[indx] = tmp_b;
}
tmp_y[indx] = tmp_y[indx] % p;
}
/* remove the added zeros. */
for(i=0; i < len_x; i++)
pOut[i] = tmp_y[i+len_a];
}
/*--end of GFFILTER()--*/
/*
* gftrunc(a)
* Iwork --- *len_a
*/
static void gftrunc(pa, len_a, len_p, Iwork)
int *pa, *len_a, len_p, *Iwork;
{
int i, len_ind, cut_0;
int *ind;
/*if isempty(a)
* c = a;
*else
* cut_0 = 1;
* if nargin > 1
* if length(tp) > 1
* cut_0 = 0;
* end;
* end;
*/
ind = Iwork;
if ( len_p > 1 )
cut_0 = 0;
else
cut_0 = 1;
/* if cut_0
* ind = find(a ~= 0);
* else
* ind = find(a >= 0);
* end;
*/
len_ind = 0;
if ( cut_0 != 0 ){
for (i=0; i < len_a[0]; i++){
if ( pa[i] != 0 ){
ind[len_ind] = i;
len_ind++;
}
}
} else {
for (i=0; i < len_a[0]; i++){
if ( pa[i] >= 0 ){
ind[len_ind] = i;
len_ind++;
}
}
}
/* if ~isempty(ind)
* c = a(1 : ind(length(ind)));
* else
* if cut_0
* c = 0;
* else
* c = -Inf;
* end;
* end;
*end;
*/
if ( len_ind > 0 ){
len_a[0] = ind[len_ind-1]+1;
} else {
if(cut_0 != 0){
pa[0] = 0;
len_a[0] = 1;
}else{
pa[0] = -Inf;
len_a[0] = 1;
}
}
}
/*--end of GFTRUNC()--*/
/*
* gfpadd(a, b, p) adds two GF(P) polynomials A and B.
* When P is a matrix.
* Iwork --- nc+nc*mp+np
* Iwork = indx
* + *nc = tmp
* + (*nc)*mp = sum
* + np = bottom of Iwork
*/
static void gfpadd(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 prim, cal_len;
int i, j, k, minus_a, minus_b, len_indx, len_indx_a, len_indx_b;
int *indx, *tmp, *sum;
indx = Iwork;
tmp = indx + nc[0];
sum = tmp + nc[0]*mp;
/* GF(Q^M) field calculation.
* The first row of P is -Inf, second is 1 and so on.
*/
/* handle vector case
* indx = find(a > n_p - 2);
* if ~isempty(indx)
* a(indx) = rem(a(indx), n_p - 1);
* end;
* indx = find(b > n_p - 2);
* if ~isempty(indx)
* b(indx) = rem(b(indx), n_p - 1);
* end;
*/
for(i=0; i < len_a; i++){
if(pa[i] > np-2)
pa[i] = pa[i] % (np - 1);
}
for(i=0; i < len_b; i++){
if(pb[i] > np-2)
pb[i] = pb[i] % (np - 1);
}
/* if (a < 0)
* c = b;
* elseif (b < 0)
* c = a;
* else
*/
minus_a = 0;
minus_b = 0;
for(i=0; i < len_a; i++){
if(pa[i] < 0)
minus_a++;
}
for(i=0; i < len_b; i++){
if(pb[i] < 0)
minus_b++;
}
if( minus_a == len_a ){
nc[0] = len_b;
for(i=0; i < len_b; i++)
pc[i] = pb[i];
} else if(minus_b == len_b){
nc[0] = len_a;
for(i=0; i < len_a; i++)
pc[i] = pa[i];
} else {
/* prim = max(max(p)) + 1;
* indx = find(~(a >= 0));
* a(indx) = - ones(1, length(indx));
* indx = find(~(b >= 0));
* b(indx) = - ones(1, length(indx));
* len_a = length(a);
* len_b = length(b);
*/
prim = 0;
for(i=0; i < np*mp; i++){
if(prim < pp[i])
prim = pp[i] + 1;
}
for(i=0; i < len_a; i++){
if(pa[i] < 0)
pa[i] = -1;
}
for(i=0; i < len_b; i++){
if(pb[i] < 0)
pb[i] = -1;
}
/* if (len_a == len_b)
* tmp = rem(p(a + 2, :) + p(b + 2, :), prim);
* cal_len = len_a;
* elseif (len_a == 1)
* tmp = rem(ones(len_b,1)*p(a + 2, :) + p(b + 2, :), prim);
* cal_len = len_b;
* elseif (len_b == 1)
* tmp = rem(p(a + 2, :) + ones(len_a, 1)*p(b + 2, :), prim);
* cal_len = len_a;
* else
* cal_len = min(len_a, len_b);
* tmp = rem(p(a(1:cal_len) + 2, :) + p(b(1:cal_len) + 2, :), prim);
* end;
*/
if( len_a == len_b && 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[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 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){
for (i=len_a; i 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 1
* b(na+nb-1) = 0;
* end
* c = gffilter(a, 1, b, p);
* end
*/
pfa = Iwork+1;
for(i=0; i= 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 */
/*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 1){
len_r[0] = len_a-1;
for(i=0; 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= 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++){
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= 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 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 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);
}