flight/libraries/rscode/berlekamp.c

   1 /***********************************************************************
   2  * Copyright Henry Minsky (hqm@alum.mit.edu) 1991-2009
   3  *
   4  * This software library is licensed under terms of the GNU GENERAL
   5  * PUBLIC LICENSE
   6  *
   7  *
   8  * RSCODE is free software: you can redistribute it and/or modify
   9  * it under the terms of the GNU General Public License as published by
  10  * the Free Software Foundation, either version 3 of the License, or
  11  * (at your option) any later version.
  12  *
  13  * RSCODE is distributed in the hope that it will be useful,
  14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  * GNU General Public License for more details.
  17  *
  18  * You should have received a copy of the GNU General Public License
  19  * along with Rscode.  If not, see <http://www.gnu.org/licenses/>.
  20  *
  21  * Commercial licensing is available under a separate license, please
  22  * contact author for details.
  23  *
  24  * Source code is available at http://rscode.sourceforge.net
  25  * Berlekamp-Peterson and Berlekamp-Massey Algorithms for error-location
  26  *
  27  * From Cain, Clark, "Error-Correction Coding For Digital Communications", pp. 205.
  28  *
  29  * This finds the coefficients of the error locator polynomial.
  30  *
  31  * The roots are then found by looking for the values of a^n
  32  * where evaluating the polynomial yields zero.
  33  *
  34  * Error correction is done using the error-evaluator equation  on pp 207.
  35  *
  36  */
  37
  38 #include <stdio.h>
  39 #include "ecc.h"
  40
  41 /* The Error Locator Polynomial, also known as Lambda or Sigma. Lambda[0] == 1 */
  42 static int Lambda[MAXDEG];
  43
  44 /* The Error Evaluator Polynomial */
  45 static int Omega[MAXDEG];
  46
  47 /* local ANSI declarations */
  48 static int compute_discrepancy(int lambda[], int S[], int L, int n);
  49 static void init_gamma(int gamma[]);
  50 static void compute_modified_omega (void);
  51 static void mul_z_poly (int src[]);
  52
  53 /* error locations found using Chien's search*/
  54 static int ErrorLocs[256];
  55 static int NErrors;
  56
  57 /* erasure flags */
  58 static int ErasureLocs[256];
  59 static int NErasures;
  60
  61 /* From  Cain, Clark, "Error-Correction Coding For Digital Communications", pp. 216. */
  62 void
  63 Modified_Berlekamp_Massey (void)
  64 {
  65   int n, L, L2, k, d, i;
  66   int psi[MAXDEG], psi2[MAXDEG], D[MAXDEG];
  67   int gamma[MAXDEG];
  68
  69   /* initialize Gamma, the erasure locator polynomial */
  70   init_gamma(gamma);
  71
  72   /* initialize to z */
  73   copy_poly(D, gamma);
  74   mul_z_poly(D);
  75
  76   copy_poly(psi, gamma);
  77   k = -1; L = NErasures;
  78
  79   for (n = NErasures; n < RS_ECC_NPARITY; n++) {
  80
  81     d = compute_discrepancy(psi, synBytes, L, n);
  82
  83     if (d != 0) {
  84
  85       /* psi2 = psi - d*D */
  86       for (i = 0; i < MAXDEG; i++) psi2[i] = psi[i] ^ gmult(d, D[i]);
  87
  88
  89       if (L < (n-k)) {
  90         L2 = n-k;
  91         k = n-L;
  92         /* D = scale_poly(ginv(d), psi); */
  93         for (i = 0; i < MAXDEG; i++) D[i] = gmult(psi[i], ginv(d));
  94         L = L2;
  95       }
  96
  97       /* psi = psi2 */
  98       for (i = 0; i < MAXDEG; i++) psi[i] = psi2[i];
  99     }
 100
 101     mul_z_poly(D);
 102   }
 103
 104   for(i = 0; i < MAXDEG; i++) Lambda[i] = psi[i];
 105   compute_modified_omega();
 106
 107
 108 }
 109
 110 /* given Psi (called Lambda in Modified_Berlekamp_Massey) and synBytes,
 111    compute the combined erasure/error evaluator polynomial as
 112    Psi*S mod z^4
 113   */
 114 void
 115 compute_modified_omega ()
 116 {
 117   int i;
 118   int product[MAXDEG*2];
 119
 120   mult_polys(product, Lambda, synBytes);
 121   zero_poly(Omega);
 122   for(i = 0; i < RS_ECC_NPARITY; i++) Omega[i] = product[i];
 123
 124 }
 125
 126 /* polynomial multiplication */
 127 void
 128 mult_polys (int dst[], int p1[], int p2[])
 129 {
 130   int i, j;
 131   int tmp1[MAXDEG*2];
 132
 133   for (i=0; i < (MAXDEG*2); i++) dst[i] = 0;
 134
 135   for (i = 0; i < MAXDEG; i++) {
 136     for(j=MAXDEG; j<(MAXDEG*2); j++) tmp1[j]=0;
 137
 138     /* scale tmp1 by p1[i] */
 139     for(j=0; j<MAXDEG; j++) tmp1[j]=gmult(p2[j], p1[i]);
 140     /* and mult (shift) tmp1 right by i */
 141     for (j = (MAXDEG*2)-1; j >= i; j--) tmp1[j] = tmp1[j-i];
 142     for (j = 0; j < i; j++) tmp1[j] = 0;
 143
 144     /* add into partial product */
 145     for(j=0; j < (MAXDEG*2); j++) dst[j] ^= tmp1[j];
 146   }
 147 }
 148
 149
 150
 151 /* gamma = product (1-z*a^Ij) for erasure locs Ij */
 152 void
 153 init_gamma (int gamma[])
 154 {
 155   int e, tmp[MAXDEG];
 156
 157   zero_poly(gamma);
 158   zero_poly(tmp);
 159   gamma[0] = 1;
 160
 161   for (e = 0; e < NErasures; e++) {
 162     copy_poly(tmp, gamma);
 163     scale_poly(gexp[ErasureLocs[e]], tmp);
 164     mul_z_poly(tmp);
 165     add_polys(gamma, tmp);
 166   }
 167 }
 168
 169
 170
 171 void
 172 compute_next_omega (int d, int A[], int dst[], int src[])
 173 {
 174   int i;
 175   for ( i = 0; i < MAXDEG;  i++) {
 176     dst[i] = src[i] ^ gmult(d, A[i]);
 177   }
 178 }
 179
 180
 181
 182 int
 183 compute_discrepancy (int lambda[], int S[], int L, int n)
 184 {
 185   int i, sum=0;
 186
 187   for (i = 0; i <= L; i++)
 188     sum ^= gmult(lambda[i], S[n-i]);
 189   return (sum);
 190 }
 191
 192 /********** polynomial arithmetic *******************/
 193
 194 void add_polys (int dst[], int src[])
 195 {
 196   int i;
 197   for (i = 0; i < MAXDEG; i++) dst[i] ^= src[i];
 198 }
 199
 200 void copy_poly (int dst[], int src[])
 201 {
 202   int i;
 203   for (i = 0; i < MAXDEG; i++) dst[i] = src[i];
 204 }
 205
 206 void scale_poly (int k, int poly[])
 207 {
 208   int i;
 209   for (i = 0; i < MAXDEG; i++) poly[i] = gmult(k, poly[i]);
 210 }
 211
 212
 213 void zero_poly (int poly[])
 214 {
 215   int i;
 216   for (i = 0; i < MAXDEG; i++) poly[i] = 0;
 217 }
 218
 219
 220 /* multiply by z, i.e., shift right by 1 */
 221 static void mul_z_poly (int src[])
 222 {
 223   int i;
 224   for (i = MAXDEG-1; i > 0; i--) src[i] = src[i-1];
 225   src[0] = 0;
 226 }
 227
 228
 229 /* Finds all the roots of an error-locator polynomial with coefficients
 230  * Lambda[j] by evaluating Lambda at successive values of alpha.
 231  *
 232  * This can be tested with the decoder's equations case.
 233  */
 234
 235
 236 void
 237 Find_Roots (void)
 238 {
 239   int sum, r, k;
 240   NErrors = 0;
 241
 242   for (r = 1; r < 256; r++) {
 243     sum = 0;
 244     /* evaluate lambda at r */
 245     for (k = 0; k < RS_ECC_NPARITY+1; k++) {
 246       sum ^= gmult(gexp[(k*r)%255], Lambda[k]);
 247     }
 248     if (sum == 0)
 249       {
 250         ErrorLocs[NErrors] = (255-r); NErrors++;
 251         //if (DEBUG) fprintf(stderr, "Root found at r = %d, (255-r) = %d\n", r, (255-r));
 252       }
 253   }
 254 }
 255
 256 /* Combined Erasure And Error Magnitude Computation
 257  *
 258  * Pass in the codeword, its size in bytes, as well as
 259  * an array of any known erasure locations, along the number
 260  * of these erasures.
 261  *
 262  * Evaluate Omega(actually Psi)/Lambda' at the roots
 263  * alpha^(-i) for error locs i.
 264  *
 265  * Returns 1 if everything ok, or 0 if an out-of-bounds error is found
 266  *
 267  */
 268
 269 int
 270 correct_errors_erasures (unsigned char codeword[],
 271                          int csize,
 272                          int nerasures,
 273                          int erasures[])
 274 {
 275   int r, i, j, err;
 276
 277   /* If you want to take advantage of erasure correction, be sure to
 278      set NErasures and ErasureLocs[] with the locations of erasures.
 279      */
 280   NErasures = nerasures;
 281   for (i = 0; i < NErasures; i++) ErasureLocs[i] = erasures[i];
 282
 283   Modified_Berlekamp_Massey();
 284   Find_Roots();
 285
 286
 287   if ((NErrors <= RS_ECC_NPARITY) && NErrors > 0) {
 288
 289     /* first check for illegal error locs */
 290     for (r = 0; r < NErrors; r++) {
 291       if (ErrorLocs[r] >= csize) {
 292                                 //if (DEBUG) fprintf(stderr, "Error loc i=%d outside of codeword length %d\n", i, csize);
 293         return(0);
 294       }
 295     }
 296
 297     for (r = 0; r < NErrors; r++) {
 298       int num, denom;
 299       i = ErrorLocs[r];
 300       /* evaluate Omega at alpha^(-i) */
 301
 302       num = 0;
 303       for (j = 0; j < MAXDEG; j++)
 304         num ^= gmult(Omega[j], gexp[((255-i)*j)%255]);
 305
 306       /* evaluate Lambda' (derivative) at alpha^(-i) ; all odd powers disappear */
 307       denom = 0;
 308       for (j = 1; j < MAXDEG; j += 2) {
 309         denom ^= gmult(Lambda[j], gexp[((255-i)*(j-1)) % 255]);
 310       }
 311
 312       err = gmult(num, ginv(denom));
 313       //if (DEBUG) fprintf(stderr, "Error magnitude %#x at loc %d\n", err, csize-i);
 314
 315       codeword[csize-i-1] ^= err;
 316     }
 317     return(1);
 318   }
 319   else {
 320     //if (DEBUG && NErrors) fprintf(stderr, "Uncorrectable codeword\n");
 321     return(0);
 322   }
 323 }
 324