1 /*****************************
2 *
3 *
4 * Multiplication and Arithmetic on Galois Field GF(256)
5 *
6 * From Mee, Daniel, "Magnetic Recording, Volume III", Ch. 5 by Patel.
7 *
8 * (c) 1991 Henry Minsky
9 *
10 *
11 ******************************/
12
13
14 #include <stdio.h>
15 #include <ctype.h>
16 #include <gpac/tools.h>
17 #include <gpac/internal/reedsolomon.h>
18
19 #ifdef GPAC_ENABLE_MPE
20
21
22 /* This is one of 14 irreducible polynomials
23 * of degree 8 and cycle length 255. (Ch 5, pp. 275, Magnetic Recording)
24 * The high order 1 bit is implicit */
25 /* x^8 + x^4 + x^3 + x^2 + 1 */
26 #define PPOLY 0x1D
27
28
29 int gexp[512];
30 int glog[256];
31
32
33 static void init_exp_table (void);
34
35
36 void
init_galois_tables(void)37 init_galois_tables (void)
38 {
39 /* initialize the table of powers of alpha */
40 init_exp_table();
41 }
42
43
44 static void
init_exp_table(void)45 init_exp_table (void)
46 {
47 int i, z;
48 int p1,p2,p3,p4,p5,p6,p7,p8;
49
50 p2 = p3 = p4 = p5 = p6 = p7 = p8 = 0;
51 p1 = 1;
52
53 gexp[0] = 1;
54 gexp[255] = gexp[0];
55 glog[0] = 0; /* shouldn't log[0] be an error? */
56
57 for (i = 1; i < 256; i++) {
58 int pinit = p8;
59 p8 = p7;
60 p7 = p6;
61 p6 = p5;
62 p5 = p4 ^ pinit;
63 p4 = p3 ^ pinit;
64 p3 = p2 ^ pinit;
65 p2 = p1;
66 p1 = pinit;
67 gexp[i] = p1 + p2*2 + p3*4 + p4*8 + p5*16 + p6*32 + p7*64 + p8*128;
68 gexp[i+255] = gexp[i];
69 }
70
71 for (i = 1; i < 256; i++) {
72 for (z = 0; z < 256; z++) {
73 if (gexp[z] == i) {
74 glog[i] = z;
75 break;
76 }
77 }
78 }
79 }
80
81 /* multiplication using logarithms */
gmult(int a,int b)82 int gmult(int a, int b)
83 {
84 int i,j;
85 if (a==0 || b == 0) return (0);
86 i = glog[a];
87 j = glog[b];
88 return (gexp[i+j]);
89 }
90
91
ginv(int elt)92 int ginv (int elt)
93 {
94 return (gexp[255-glog[elt]]);
95 }
96
97
98
99 /***********************************************************************
100 * Berlekamp-Peterson and Berlekamp-Massey Algorithms for error-location
101 *
102 * From Cain, Clark, "Error-Correction Coding For Digital Communications", pp. 205.
103 *
104 * This finds the coefficients of the error locator polynomial.
105 *
106 * The roots are then found by looking for the values of a^n
107 * where evaluating the polynomial yields zero.
108 *
109 * Error correction is done using the error-evaluator equation on pp 207.
110 *
111 * hqm@ai.mit.edu Henry Minsky
112 */
113
114
115 /* The Error Locator Polynomial, also known as Lambda or Sigma. Lambda[0] == 1 */
116 static int Lambda[MAXDEG];
117
118 /* The Error Evaluator Polynomial */
119 static int Omega[MAXDEG];
120
121 /* local ANSI declarations */
122 static int compute_discrepancy(int lambda[], int S[], int L, int n);
123 static void init_gamma(int gamma[]);
124 static void compute_modified_omega (void);
125 static void mul_z_poly (int src[]);
126
127 /* error locations found using Chien's search*/
128 static int ErrorLocs[256];
129 static int NErrors;
130
131 /* erasure flags */
132 static int ErasureLocs[256];
133 static int NErasures;
134
135 /* From Cain, Clark, "Error-Correction Coding For Digital Communications", pp. 216. */
136 void
Modified_Berlekamp_Massey(void)137 Modified_Berlekamp_Massey (void)
138 {
139 int n, L, L2, k, i;
140 int psi[MAXDEG], psi2[MAXDEG], D[MAXDEG];
141 int gamma[MAXDEG];
142
143 /* initialize Gamma, the erasure locator polynomial */
144 init_gamma(gamma);
145
146 /* initialize to z */
147 copy_poly(D, gamma);
148 mul_z_poly(D);
149
150 copy_poly(psi, gamma);
151 k = -1;
152 L = NErasures;
153
154 for (n = NErasures; n < NPAR; n++) {
155
156 int d = compute_discrepancy(psi, synBytes, L, n);
157
158 if (d != 0) {
159
160 /* psi2 = psi - d*D */
161 for (i = 0; i < MAXDEG; i++) psi2[i] = psi[i] ^ gmult(d, D[i]);
162
163
164 if (L < (n-k)) {
165 L2 = n-k;
166 k = n-L;
167 /* D = scale_poly(ginv(d), psi); */
168 for (i = 0; i < MAXDEG; i++) D[i] = gmult(psi[i], ginv(d));
169 L = L2;
170 }
171
172 /* psi = psi2 */
173 for (i = 0; i < MAXDEG; i++) psi[i] = psi2[i];
174 }
175
176 mul_z_poly(D);
177 }
178
179 for(i = 0; i < MAXDEG; i++) Lambda[i] = psi[i];
180 compute_modified_omega();
181
182
183 }
184
185 /* given Psi (called Lambda in Modified_Berlekamp_Massey) and synBytes,
186 compute the combined erasure/error evaluator polynomial as
187 Psi*S mod z^4
188 */
189 void
compute_modified_omega()190 compute_modified_omega ()
191 {
192 int i;
193 int product[MAXDEG*2];
194
195 mult_polys(product, Lambda, synBytes);
196 zero_poly(Omega);
197 for(i = 0; i < NPAR; i++) Omega[i] = product[i];
198
199 }
200
201 /* polynomial multiplication */
202 void
mult_polys(int dst[],int p1[],int p2[])203 mult_polys (int dst[], int p1[], int p2[])
204 {
205 int i, j;
206 int tmp1[MAXDEG*2];
207
208 for (i=0; i < (MAXDEG*2); i++) dst[i] = 0;
209
210 for (i = 0; i < MAXDEG; i++) {
211 for(j=MAXDEG; j<(MAXDEG*2); j++) tmp1[j]=0;
212
213 /* scale tmp1 by p1[i] */
214 for(j=0; j<MAXDEG; j++) tmp1[j]=gmult(p2[j], p1[i]);
215 /* and mult (shift) tmp1 right by i */
216 for (j = (MAXDEG*2)-1; j >= i; j--) tmp1[j] = tmp1[j-i];
217 for (j = 0; j < i; j++) tmp1[j] = 0;
218
219 /* add into partial product */
220 for(j=0; j < (MAXDEG*2); j++) dst[j] ^= tmp1[j];
221 }
222 }
223
224
225
226 /* gamma = product (1-z*a^Ij) for erasure locs Ij */
227 void
init_gamma(int gamma[])228 init_gamma (int gamma[])
229 {
230 int e, tmp[MAXDEG];
231
232 zero_poly(gamma);
233 zero_poly(tmp);
234 gamma[0] = 1;
235
236 for (e = 0; e < NErasures; e++) {
237 copy_poly(tmp, gamma);
238 scale_poly(gexp[ErasureLocs[e]], tmp);
239 mul_z_poly(tmp);
240 add_polys(gamma, tmp);
241 }
242 }
243
244
245
246 void
compute_next_omega(int d,int A[],int dst[],int src[])247 compute_next_omega (int d, int A[], int dst[], int src[])
248 {
249 int i;
250 for ( i = 0; i < MAXDEG; i++) {
251 dst[i] = src[i] ^ gmult(d, A[i]);
252 }
253 }
254
255
256
257 int
compute_discrepancy(int lambda[],int S[],int L,int n)258 compute_discrepancy (int lambda[], int S[], int L, int n)
259 {
260 int i, sum=0;
261
262 for (i = 0; i <= L; i++)
263 sum ^= gmult(lambda[i], S[n-i]);
264 return (sum);
265 }
266
267 /********** polynomial arithmetic *******************/
268
add_polys(int dst[],int src[])269 void add_polys (int dst[], int src[])
270 {
271 int i;
272 for (i = 0; i < MAXDEG; i++) dst[i] ^= src[i];
273 }
274
copy_poly(int dst[],int src[])275 void copy_poly (int dst[], int src[])
276 {
277 int i;
278 for (i = 0; i < MAXDEG; i++) dst[i] = src[i];
279 }
280
scale_poly(int k,int poly[])281 void scale_poly (int k, int poly[])
282 {
283 int i;
284 for (i = 0; i < MAXDEG; i++) poly[i] = gmult(k, poly[i]);
285 }
286
287
zero_poly(int poly[])288 void zero_poly (int poly[])
289 {
290 int i;
291 for (i = 0; i < MAXDEG; i++) poly[i] = 0;
292 }
293
294
295 /* multiply by z, i.e., shift right by 1 */
mul_z_poly(int src[])296 static void mul_z_poly (int src[])
297 {
298 int i;
299 for (i = MAXDEG-1; i > 0; i--) src[i] = src[i-1];
300 src[0] = 0;
301 }
302
303
304 /* Finds all the roots of an error-locator polynomial with coefficients
305 * Lambda[j] by evaluating Lambda at successive values of alpha.
306 *
307 * This can be tested with the decoder's equations case.
308 */
309
310
311 void
Find_Roots(void)312 Find_Roots (void)
313 {
314 int r, k;
315 NErrors = 0;
316
317 for (r = 1; r < 256; r++) {
318 int sum = 0;
319 /* evaluate lambda at r */
320 for (k = 0; k < NPAR+1; k++) {
321 sum ^= gmult(gexp[(k*r)%255], Lambda[k]);
322 }
323 if (sum == 0)
324 {
325 ErrorLocs[NErrors] = (255-r);
326 NErrors++;
327 if (RS_DEBUG) fprintf(stderr, "Root found at r = %d, (255-r) = %d\n", r, (255-r));
328 }
329 }
330 }
331
332 /* Combined Erasure And Error Magnitude Computation
333 *
334 * Pass in the codeword, its size in bytes, as well as
335 * an array of any known erasure locations, along the number
336 * of these erasures.
337 *
338 * Evaluate Omega(actually Psi)/Lambda' at the roots
339 * alpha^(-i) for error locs i.
340 *
341 * Returns 1 if everything ok, or 0 if an out-of-bounds error is found
342 *
343 */
344
345 int
correct_errors_erasures(unsigned char codeword[],int csize,int nerasures,int erasures[])346 correct_errors_erasures (unsigned char codeword[],
347 int csize,
348 int nerasures,
349 int erasures[])
350 {
351 int i;
352
353 /* If you want to take advantage of erasure correction, be sure to
354 set NErasures and ErasureLocs[] with the locations of erasures.
355 */
356 NErasures = nerasures;
357 for (i = 0; i < NErasures; i++) ErasureLocs[i] = erasures[i];
358
359 Modified_Berlekamp_Massey();
360 Find_Roots();
361
362
363 if ((NErrors <= NPAR) && NErrors > 0) {
364 int r;
365
366 /* first check for illegal error locs */
367 for (r = 0; r < NErrors; r++) {
368 if (ErrorLocs[r] >= csize) {
369 if (RS_DEBUG) fprintf(stderr, "Error loc i=%d outside of codeword length %d\n", i, csize);
370 return(0);
371 }
372 }
373
374 for (r = 0; r < NErrors; r++) {
375 int num, denom, err, j;
376 i = ErrorLocs[r];
377 /* evaluate Omega at alpha^(-i) */
378
379 num = 0;
380 for (j = 0; j < MAXDEG; j++)
381 num ^= gmult(Omega[j], gexp[((255-i)*j)%255]);
382
383 /* evaluate Lambda' (derivative) at alpha^(-i) ; all odd powers disappear */
384 denom = 0;
385 for (j = 1; j < MAXDEG; j += 2) {
386 denom ^= gmult(Lambda[j], gexp[((255-i)*(j-1)) % 255]);
387 }
388
389 err = gmult(num, ginv(denom));
390 if (RS_DEBUG) fprintf(stderr, "Error magnitude %#x at loc %d\n", err, csize-i);
391
392 codeword[csize-i-1] ^= err;
393 }
394 return(1);
395 }
396 else {
397 if (RS_DEBUG && NErrors) fprintf(stderr, "Uncorrectable codeword\n");
398 return(0);
399 }
400 }
401
402
403
404 /*
405 * Reed Solomon Encoder/Decoder
406 *
407 * (c) Henry Minsky (hqm@ua.com), Universal Access 1991-1995
408 */
409
410 /* Encoder parity bytes */
411 int pBytes[MAXDEG];
412
413 /* Decoder syndrome bytes */
414 int synBytes[MAXDEG];
415
416 /* generator polynomial */
417 int genPoly[MAXDEG*2];
418
419 int RS_DEBUG = FALSE;
420
421 static void
422 compute_genpoly (int nbytes, int genpoly[]);
423
424 /* Initialize lookup tables, polynomials, etc. */
425 void
initialize_ecc()426 initialize_ecc ()
427 {
428 /* Initialize the galois field arithmetic tables */
429 init_galois_tables();
430
431 /* Compute the encoder generator polynomial */
432 compute_genpoly(NPAR, genPoly);
433 }
434
435 void
zero_fill_from(unsigned char buf[],int from,int to)436 zero_fill_from (unsigned char buf[], int from, int to)
437 {
438 int i;
439 for (i = from; i < to; i++) buf[i] = 0;
440 }
441
442 /* debugging routines */
443 void
print_parity(void)444 print_parity (void)
445 {
446 int i;
447 fprintf(stderr, "Parity Bytes: ");
448 for (i = 0; i < NPAR; i++)
449 fprintf(stderr, "[%d]:%x, ",i,pBytes[i]);
450 fprintf(stderr, "\n");
451 }
452
453
454 void
print_syndrome(void)455 print_syndrome (void)
456 {
457 int i;
458 fprintf(stderr, "Syndrome Bytes: ");
459 for (i = 0; i < NPAR; i++)
460 fprintf(stderr, "[%d]:%x, ",i,synBytes[i]);
461 fprintf(stderr, "\n");
462 }
463
464 /* Append the parity bytes onto the end of the message */
465 void
build_codeword(unsigned char msg[],int nbytes,unsigned char dst[])466 build_codeword (unsigned char msg[], int nbytes, unsigned char dst[])
467 {
468 int i;
469
470 for (i = 0; i < nbytes; i++) dst[i] = msg[i];
471
472 for (i = 0; i < NPAR; i++) {
473 dst[i+nbytes] = pBytes[NPAR-1-i];
474 }
475 }
476
477 /**********************************************************
478 * Reed Solomon Decoder
479 *
480 * Computes the syndrome of a codeword. Puts the results
481 * into the synBytes[] array.
482 */
483
484 void
decode_data(unsigned char data[],int nbytes)485 decode_data(unsigned char data[], int nbytes)
486 {
487 int i, j;
488 for (j = 0; j < NPAR; j++) {
489 int sum = 0;
490 for (i = 0; i < nbytes; i++) {
491 sum = data[i] ^ gmult(gexp[j+1], sum);
492 }
493 synBytes[j] = sum;
494 }
495 }
496
497
498 /* Check if the syndrome is zero */
499 int
check_syndrome(void)500 check_syndrome (void)
501 {
502 int i, nz = 0;
503 for (i =0 ; i < NPAR; i++) {
504 if (synBytes[i] != 0) nz = 1;
505 }
506 return nz;
507 }
508
509
510 void
debug_check_syndrome(void)511 debug_check_syndrome (void)
512 {
513 int i;
514
515 for (i = 0; i < 3; i++) {
516 fprintf(stderr, " inv log S[%d]/S[%d] = %d\n", i, i+1,
517 glog[gmult(synBytes[i], ginv(synBytes[i+1]))]);
518 }
519 }
520
521
522 /* Create a generator polynomial for an n byte RS code.
523 * The coefficients are returned in the genPoly arg.
524 * Make sure that the genPoly array which is passed in is
525 * at least n+1 bytes long.
526 */
527
528 static void
compute_genpoly(int nbytes,int genpoly[])529 compute_genpoly (int nbytes, int genpoly[])
530 {
531 int i, tp[256], tp1[256];
532
533 /* multiply (x + a^n) for n = 1 to nbytes */
534
535 zero_poly(tp1);
536 tp1[0] = 1;
537
538 for (i = 1; i <= nbytes; i++) {
539 zero_poly(tp);
540 tp[0] = gexp[i]; /* set up x+a^n */
541 tp[1] = 1;
542
543 mult_polys(genpoly, tp, tp1);
544 copy_poly(tp1, genpoly);
545 }
546 }
547
548 /* Simulate a LFSR with generator polynomial for n byte RS code.
549 * Pass in a pointer to the data array, and amount of data.
550 *
551 * The parity bytes are deposited into pBytes[], and the whole message
552 * and parity are copied to dest to make a codeword.
553 *
554 */
555
556 void
encode_data(unsigned char msg[],int nbytes,unsigned char dst[])557 encode_data (unsigned char msg[], int nbytes, unsigned char dst[])
558 {
559 int i, LFSR[NPAR+1],j;
560
561 for(i=0; i < NPAR+1; i++) LFSR[i]=0;
562
563 for (i = 0; i < nbytes; i++) {
564 int dbyte = msg[i] ^ LFSR[NPAR-1];
565 for (j = NPAR-1; j > 0; j--) {
566 LFSR[j] = LFSR[j-1] ^ gmult(genPoly[j], dbyte);
567 }
568 LFSR[0] = gmult(genPoly[0], dbyte);
569 }
570
571 for (i = 0; i < NPAR; i++)
572 pBytes[i] = LFSR[i];
573
574 build_codeword(msg, nbytes, dst);
575 }
576
577 #endif //GPAC_ENABLE_MPE
578