libfriidump/rs.c

   1 #include <stdlib.h>
   2 #include <math.h>
   3 #include <stdio.h>
   4 #include <string.h>
   5
   6 #define mm  8           /* RS code over GF(2**mm) - change to suit */
   7 #define n   256             /* n = size of the field */
   8 #define nn  182         /* nn=2**mm -1   length of codeword */
   9 #define kk  172         /* kk = nn-2*tt  */ /* Degree of g(x) = 2*tt */
  10
  11 //#define       NN              n-1
  12 //#define       FCR             0
  13 //#define       PRIM    1
  14 #define _NROOTS nn-kk
  15 //#define       PAD             NN-nn
  16 //#define       A0              NN
  17 //#define       IPRIM   1
  18
  19 const int       NN     = n-1;
  20 const int       FCR    = 0;
  21 const int       PRIM   = 1;
  22 const int       NROOTS = nn-kk;
  23 const int       PAD    = (n-1)-nn;
  24 const int       A0     = n-1;
  25 const int       IPRIM  = 1;
  26
  27
  28 #ifndef min
  29 #define min(a,b)        ((a) < (b) ? (a) : (b))
  30 #endif
  31
  32 /**** Primitive polynomial ****/
  33 int pp [mm+1] = { 1, 0, 1, 1, 1, 0, 0, 0, 1}; /* 1+x^2+x^3+x^4+x^8 */
  34
  35 /* generator polynomial, tables for Galois field */
  36 int alpha_to[n], index_of[n], gg[nn-kk+1];
  37
  38 int b0 = 1;
  39
  40 /* data[] is the info vector, bb[] is the parity vector, recd[] is the
  41   noise corrupted received vector  */
  42 int recd[nn], data[kk], bb[nn-kk];
  43
  44 int modnn(int x){
  45   while (x >= 0xff) {
  46     x -= 0xff;
  47     x = (x >> 0xff) + (x & 0xff);
  48   }
  49   return x;
  50 }
  51
  52
  53 void generate_gf()
  54  {
  55         register int i, mask ;
  56
  57   mask = 1 ;
  58   alpha_to[mm] = 0 ;
  59   for (i=0; i<mm; i++)
  60    { alpha_to[i] = mask ;
  61      index_of[alpha_to[i]] = i ;
  62      if (pp[i]!=0) /* If pp[i] == 1 then, term @^i occurs in poly-repr of @^mm */
  63        alpha_to[mm] ^= mask ;  /* Bit-wise EXOR operation */
  64      mask <<= 1 ; /* single left-shift */
  65    }
  66   index_of[alpha_to[mm]] = mm ;
  67   /* Have obtained poly-repr of @^mm. Poly-repr of @^(i+1) is given by
  68      poly-repr of @^i shifted left one-bit and accounting for any @^mm
  69      term that may occur when poly-repr of @^i is shifted. */
  70   mask >>= 1 ;
  71   for (i=mm+1; i<255; i++)
  72    { if (alpha_to[i-1] >= mask)
  73         alpha_to[i] = alpha_to[mm] ^ ((alpha_to[i-1]^mask)<<1) ;
  74      else alpha_to[i] = alpha_to[i-1]<<1 ;
  75      index_of[alpha_to[i]] = i ;
  76    }
  77   index_of[0] = A0 ;//-1
  78  }
  79
  80
  81 void gen_poly()
  82 /* Obtain the generator polynomial of the tt-error correcting, length */
  83  {
  84         register int i, j, root;
  85
  86         gg[0] = 1;
  87
  88         for (i = 0,root=0*1; i < nn-kk; i++,root += 1) {
  89                 gg[i+1] = 1;
  90
  91                 for (j = i; j > 0; j--){
  92                         if (gg[j] != 0)
  93                                 gg[j] = gg[j-1] ^ alpha_to[modnn(index_of[gg[j]] + root)];
  94                         else
  95                                 gg[j] = gg[j-1];
  96                 }
  97
  98                 gg[0] = alpha_to[modnn(index_of[gg[0]] + root)];
  99         }
 100         for (i=0; i <= nn-kk; i++) {
 101                 gg[i] = index_of[gg[i]];
 102         }
 103  }
 104
 105
 106 void rs_encode(unsigned char *data, unsigned char *bb)
 107  {
 108         register int i,j ;
 109         int feedback;
 110
 111         for (i=0; i<NROOTS; i++)   bb[i] = 0; //nullify result
 112
 113         for(i=0;i<NN-NROOTS-PAD;i++){
 114                 feedback = index_of[data[i] ^ bb[0]];
 115
 116                 if(feedback != A0){      /* feedback term is non-zero */
 117                         for(j=1;j<NROOTS;j++) {
 118                                 bb[j] ^= alpha_to[modnn(feedback + gg[NROOTS-j])];
 119                         }
 120                 }
 121                 /* Shift */
 122                 memmove(&bb[0],&bb[1], NROOTS-1);
 123                 //for (j=0; j<NROOTS-1; j++)   bb[j] = bb[j+1];
 124
 125                 if(feedback != A0)
 126                         bb[NROOTS-1] = alpha_to[modnn(feedback + gg[0])];
 127                 else
 128                         bb[NROOTS-1] = 0;
 129         }
 130  }
 131 ///*
 132 int rs_decode(unsigned char *data, int *eras_pos, int no_eras){
 133   int deg_lambda, el, deg_omega;
 134   int i, j, r,k;
 135   unsigned char u,q,tmp,num1,num2,den,discr_r;
 136   unsigned char lambda[_NROOTS+1], s[_NROOTS];
 137   unsigned char b[_NROOTS+1], t[_NROOTS+1], omega[_NROOTS+1];
 138   unsigned char root[_NROOTS], reg[_NROOTS+1], loc[_NROOTS];
 139   int syn_error, count;
 140
 141
 142   // form the syndromes; i.e., evaluate data(x) at roots of g(x)
 143   for(i=0;i<NROOTS;i++)
 144     s[i] = data[0];
 145
 146   for(j=1;j<NN-PAD;j++){
 147     for(i=0;i<NROOTS;i++){
 148       if(s[i] == 0){
 149         s[i] = data[j];
 150       } else {
 151         s[i] = data[j] ^ alpha_to[modnn(index_of[s[i]] + (FCR+i)*PRIM)];
 152       }
 153     }
 154   }
 155
 156   // Convert syndromes to index form, checking for nonzero condition
 157   syn_error = 0;
 158   for(i=0;i<NROOTS;i++){
 159     syn_error |= s[i];
 160     s[i] = index_of[s[i]];
 161   }
 162
 163   if (!syn_error) {
 164     // if syndrome is zero, data[] is a codeword and there are no
 165     // errors to correct. So return data[] unmodified
 166     //
 167     count = 0;
 168     goto finish;
 169   }
 170   memset(&lambda[1],0,NROOTS*sizeof(lambda[0]));
 171   lambda[0] = 1;
 172
 173   if (no_eras > 0) {
 174     /* Init lambda to be the erasure locator polynomial */
 175     lambda[1] = alpha_to[modnn(PRIM*(NN-1-eras_pos[0]))];
 176     for (i = 1; i < no_eras; i++) {
 177       u = modnn(PRIM*(NN-1-eras_pos[i]));
 178       for (j = i+1; j > 0; j--) {
 179         tmp = index_of[lambda[j - 1]];
 180         if(tmp != A0)
 181           lambda[j] ^= alpha_to[modnn(u + tmp)];
 182       }
 183     }
 184   }
 185   for(i=0;i<NROOTS+1;i++)
 186     b[i] = index_of[lambda[i]];
 187
 188   /*
 189    * Begin Berlekamp-Massey algorithm to determine error+erasure
 190    * locator polynomial
 191    */
 192   r = no_eras;
 193   el = no_eras;
 194
 195   while (++r <= NROOTS) {       /* r is the step number */
 196     /* Compute discrepancy at the r-th step in poly-form */
 197     discr_r = 0;
 198     for (i = 0; i < r; i++){
 199       if ((lambda[i] != 0) && (s[r-i-1] != A0)) {
 200         discr_r ^= alpha_to[modnn(index_of[lambda[i]] + s[r-i-1])];
 201       }
 202     }
 203     discr_r = index_of[discr_r];        /* Index form */
 204     if (discr_r == A0) {
 205       /* 2 lines below: B(x) <-- x*B(x) */
 206       memmove(&b[1],b,NROOTS*sizeof(b[0]));
 207       b[0] = A0;
 208     } else {
 209       /* 7 lines below: T(x) <-- lambda(x) - discr_r*x*b(x) */
 210       t[0] = lambda[0];
 211       for (i = 0 ; i < NROOTS; i++) {
 212         if(b[i] != A0)
 213           t[i+1] = lambda[i+1] ^ alpha_to[modnn(discr_r + b[i])];
 214         else
 215           t[i+1] = lambda[i+1];
 216       }
 217       if (2 * el <= r + no_eras - 1) {
 218         el = r + no_eras - el;
 219         /*
 220          * 2 lines below: B(x) <-- inv(discr_r) *
 221          * lambda(x)
 222          */
 223         for (i = 0; i <= NROOTS; i++)
 224           b[i] = (lambda[i] == 0) ? A0 : modnn(index_of[lambda[i]] - discr_r + NN);
 225       } else {
 226         /* 2 lines below: B(x) <-- x*B(x) */
 227         memmove(&b[1],b,NROOTS*sizeof(b[0]));
 228         b[0] = A0;
 229       }
 230       memcpy(lambda,t,(NROOTS+1)*sizeof(t[0]));
 231     }
 232   }
 233
 234   /* Convert lambda to index form and compute deg(lambda(x)) */
 235   deg_lambda = 0;
 236   for(i=0;i<NROOTS+1;i++){
 237     lambda[i] = index_of[lambda[i]];
 238     if(lambda[i] != A0)
 239       deg_lambda = i;
 240   }
 241   /* Find roots of the error+erasure locator polynomial by Chien search */
 242   memcpy(&reg[1],&lambda[1],NROOTS*sizeof(reg[0]));
 243   count = 0;            /* Number of roots of lambda(x) */
 244   for (i = 1,k=IPRIM-1; i <= NN; i++,k = modnn(k+IPRIM)) {
 245     q = 1; /* lambda[0] is always 0 */
 246     for (j = deg_lambda; j > 0; j--){
 247       if (reg[j] != A0) {
 248         reg[j] = modnn(reg[j] + j);
 249         q ^= alpha_to[reg[j]];
 250       }
 251     }
 252     if (q != 0)
 253       continue; /* Not a root */
 254     /* store root (index-form) and error location number */
 255     root[count] = i;
 256     loc[count] = k;
 257     /* If we've already found max possible roots,
 258      * abort the search to save time
 259      */
 260     if(++count == deg_lambda)
 261       break;
 262   }
 263
 264   if (deg_lambda != count) {
 265     /*
 266      * deg(lambda) unequal to number of roots => uncorrectable
 267      * error detected
 268      */
 269     count = -1;
 270     goto finish;
 271   }
 272   /*
 273    * Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo
 274    * x**NROOTS). in index form. Also find deg(omega).
 275    */
 276   deg_omega = deg_lambda-1;
 277   for (i = 0; i <= deg_omega;i++){
 278     tmp = 0;
 279     for(j=i;j >= 0; j--){
 280       if ((s[i - j] != A0) && (lambda[j] != A0))
 281         tmp ^= alpha_to[modnn(s[i - j] + lambda[j])];
 282     }
 283     omega[i] = index_of[tmp];
 284   }
 285
 286   /*
 287    * Compute error values in poly-form. num1 = omega(inv(X(l))), num2 =
 288    * inv(X(l))**(FCR-1) and den = lambda_pr(inv(X(l))) all in poly-form
 289    */
 290   for (j = count-1; j >=0; j--) {
 291     num1 = 0;
 292     for (i = deg_omega; i >= 0; i--) {
 293       if (omega[i] != A0)
 294         num1  ^= alpha_to[modnn(omega[i] + i * root[j])];
 295     }
 296     num2 = alpha_to[modnn(root[j] * (FCR - 1) + NN)];
 297     den = 0;
 298
 299     /* lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i] */
 300     for (i = min(deg_lambda,NROOTS-1) & ~1; i >= 0; i -=2) {
 301       if(lambda[i+1] != A0)
 302         den ^= alpha_to[modnn(lambda[i+1] + i * root[j])];
 303     }
 304     /* Apply error to data */
 305     if (num1 != 0 && loc[j] >= PAD) {
 306       data[loc[j]-PAD] ^= alpha_to[modnn(index_of[num1] + index_of[num2] + NN - index_of[den])];
 307     }
 308   }
 309
 310  finish:
 311   if(eras_pos != NULL){
 312     for(i=0;i<count;i++)
 313       eras_pos[i] = loc[i];
 314   }
 315   return count;
 316 }