1 /*Copyright (C) 1991-1995 Henry Minsky (hqm@ua.com, hqm@ai.mit.edu)
2 Copyright (C) 2008-2009 Timothy B. Terriberry (tterribe@xiph.org)
3 You can redistribute this library and/or modify it under the terms of the
4 GNU Lesser General Public License as published by the Free Software
5 Foundation; either version 2.1 of the License, or (at your option) any later
6 version.*/
7 #include <stdlib.h>
8 #include <string.h>
9 #include "rs.h"
10
11 /*Reed-Solomon encoder and decoder.
12 Original implementation (C) Henry Minsky (hqm@ua.com, hqm@ai.mit.edu),
13 Universal Access 1991-1995.
14 Updates by Timothy B. Terriberry (C) 2008-2009:
15 - Properly reject codes when error-locator polynomial has repeated roots or
16 non-trivial irreducible factors.
17 - Removed the hard-coded parity size and performed general API cleanup.
18 - Allow multiple representations of GF(2**8), since different standards use
19 different irreducible polynomials.
20 - Allow different starting indices for the generator polynomial, since
21 different standards use different values.
22 - Greatly reduced the computation by eliminating unnecessary operations.
23 - Explicitly solve for the roots of low-degree polynomials instead of using
24 an exhaustive search.
25 This is another major speed boost when there are few errors.*/
26
27
28 /*Galois Field arithmetic in GF(2**8).*/
29
rs_gf256_init(rs_gf256 * _gf,unsigned _ppoly)30 void rs_gf256_init(rs_gf256 *_gf,unsigned _ppoly){
31 unsigned p;
32 int i;
33 /*Initialize the table of powers of a primtive root, alpha=0x02.*/
34 p=1;
35 for(i=0;i<256;i++){
36 _gf->exp[i]=_gf->exp[i+255]=p;
37 p=((p<<1)^(-(p>>7)&_ppoly))&0xFF;
38 }
39 /*Invert the table to recover the logs.*/
40 for(i=0;i<255;i++)_gf->log[_gf->exp[i]]=i;
41 /*Note that we rely on the fact that _gf->log[0]=0 below.*/
42 _gf->log[0]=0;
43 }
44
45 /*Multiplication in GF(2**8) using logarithms.*/
rs_gmul(const rs_gf256 * _gf,unsigned _a,unsigned _b)46 static unsigned rs_gmul(const rs_gf256 *_gf,unsigned _a,unsigned _b){
47 return _a==0||_b==0?0:_gf->exp[_gf->log[_a]+_gf->log[_b]];
48 }
49
50 /*Division in GF(2**8) using logarithms.
51 The result of division by zero is undefined.*/
rs_gdiv(const rs_gf256 * _gf,unsigned _a,unsigned _b)52 static unsigned rs_gdiv(const rs_gf256 *_gf,unsigned _a,unsigned _b){
53 return _a==0?0:_gf->exp[_gf->log[_a]+255-_gf->log[_b]];
54 }
55
56 /*Multiplication in GF(2**8) when one of the numbers is known to be non-zero
57 (proven by representing it by its logarithm).*/
rs_hgmul(const rs_gf256 * _gf,unsigned _a,unsigned _logb)58 static unsigned rs_hgmul(const rs_gf256 *_gf,unsigned _a,unsigned _logb){
59 return _a==0?0:_gf->exp[_gf->log[_a]+_logb];
60 }
61
62 /*Square root in GF(2**8) using logarithms.*/
rs_gsqrt(const rs_gf256 * _gf,unsigned _a)63 static unsigned rs_gsqrt(const rs_gf256 *_gf,unsigned _a){
64 unsigned loga;
65 if(!_a)return 0;
66 loga=_gf->log[_a];
67 return _gf->exp[loga+(255&-(loga&1))>>1];
68 }
69
70 /*Polynomial root finding in GF(2**8).
71 Each routine returns a list of the distinct roots (i.e., with duplicates
72 removed).*/
73
74 /*Solve a quadratic equation x**2 + _b*x + _c in GF(2**8) using the method
75 of~\cite{Wal99}.
76 Returns the number of distinct roots.
77 ARTICLE{Wal99,
78 author="C. Wayne Walker",
79 title="New Formulas for Solving Quadratic Equations over Certain Finite
80 Fields",
81 journal="{IEEE} Transactions on Information Theory",
82 volume=45,
83 number=1,
84 pages="283--284",
85 month=Jan,
86 year=1999
87 }*/
rs_quadratic_solve(const rs_gf256 * _gf,unsigned _b,unsigned _c,unsigned char _x[2])88 static int rs_quadratic_solve(const rs_gf256 *_gf,unsigned _b,unsigned _c,
89 unsigned char _x[2]){
90 unsigned b;
91 unsigned logb;
92 unsigned logb2;
93 unsigned logb4;
94 unsigned logb8;
95 unsigned logb12;
96 unsigned logb14;
97 unsigned logc;
98 unsigned logc2;
99 unsigned logc4;
100 unsigned c8;
101 unsigned g3;
102 unsigned z3;
103 unsigned l3;
104 unsigned c0;
105 unsigned g2;
106 unsigned l2;
107 unsigned z2;
108 int inc;
109 /*If _b is zero, all we need is a square root.*/
110 if(!_b){
111 _x[0]=rs_gsqrt(_gf,_c);
112 return 1;
113 }
114 /*If _c is zero, 0 and _b are the roots.*/
115 if(!_c){
116 _x[0]=0;
117 _x[1]=_b;
118 return 2;
119 }
120 logb=_gf->log[_b];
121 logc=_gf->log[_c];
122 /*If _b lies in GF(2**4), scale x to move it out.*/
123 inc=logb%(255/15)==0;
124 if(inc){
125 b=_gf->exp[logb+254];
126 logb=_gf->log[b];
127 _c=_gf->exp[logc+253];
128 logc=_gf->log[_c];
129 }
130 else b=_b;
131 logb2=_gf->log[_gf->exp[logb<<1]];
132 logb4=_gf->log[_gf->exp[logb2<<1]];
133 logb8=_gf->log[_gf->exp[logb4<<1]];
134 logb12=_gf->log[_gf->exp[logb4+logb8]];
135 logb14=_gf->log[_gf->exp[logb2+logb12]];
136 logc2=_gf->log[_gf->exp[logc<<1]];
137 logc4=_gf->log[_gf->exp[logc2<<1]];
138 c8=_gf->exp[logc4<<1];
139 g3=rs_hgmul(_gf,
140 _gf->exp[logb14+logc]^_gf->exp[logb12+logc2]^_gf->exp[logb8+logc4]^c8,logb);
141 /*If g3 doesn't lie in GF(2**4), then our roots lie in an extension field.
142 Note that we rely on the fact that _gf->log[0]==0 here.*/
143 if(_gf->log[g3]%(255/15)!=0)return 0;
144 /*Construct the corresponding quadratic in GF(2**4):
145 x**2 + x/alpha**(255/15) + l3/alpha**(2*(255/15))*/
146 z3=rs_gdiv(_gf,g3,_gf->exp[logb8<<1]^b);
147 l3=rs_hgmul(_gf,rs_gmul(_gf,z3,z3)^rs_hgmul(_gf,z3,logb)^_c,255-logb2);
148 c0=rs_hgmul(_gf,l3,255-2*(255/15));
149 /*Construct the corresponding quadratic in GF(2**2):
150 x**2 + x/alpha**(255/3) + l2/alpha**(2*(255/3))*/
151 g2=rs_hgmul(_gf,
152 rs_hgmul(_gf,c0,255-2*(255/15))^rs_gmul(_gf,c0,c0),255-255/15);
153 z2=rs_gdiv(_gf,g2,_gf->exp[255-(255/15)*4]^_gf->exp[255-(255/15)]);
154 l2=rs_hgmul(_gf,
155 rs_gmul(_gf,z2,z2)^rs_hgmul(_gf,z2,255-(255/15))^c0,2*(255/15));
156 /*Back substitute to the solution in the original field.*/
157 _x[0]=_gf->exp[_gf->log[z3^rs_hgmul(_gf,
158 rs_hgmul(_gf,l2,255/3)^rs_hgmul(_gf,z2,255/15),logb)]+inc];
159 _x[1]=_x[0]^_b;
160 return 2;
161 }
162
163 /*Solve a cubic equation x**3 + _a*x**2 + _b*x + _c in GF(2**8).
164 Returns the number of distinct roots.*/
rs_cubic_solve(const rs_gf256 * _gf,unsigned _a,unsigned _b,unsigned _c,unsigned char _x[3])165 static int rs_cubic_solve(const rs_gf256 *_gf,
166 unsigned _a,unsigned _b,unsigned _c,unsigned char _x[3]){
167 unsigned k;
168 unsigned logd;
169 unsigned d2;
170 unsigned logd2;
171 unsigned logw;
172 int nroots;
173 /*If _c is zero, factor out the 0 root.*/
174 if(!_c){
175 nroots=rs_quadratic_solve(_gf,_a,_b,_x);
176 if(_b)_x[nroots++]=0;
177 return nroots;
178 }
179 /*Substitute x=_a+y*sqrt(_a**2+_b) to get y**3 + y + k == 0,
180 k = (_a*_b+c)/(_a**2+b)**(3/2).*/
181 k=rs_gmul(_gf,_a,_b)^_c;
182 d2=rs_gmul(_gf,_a,_a)^_b;
183 if(!d2){
184 int logx;
185 if(!k){
186 /*We have a triple root.*/
187 _x[0]=_a;
188 return 1;
189 }
190 logx=_gf->log[k];
191 if(logx%3!=0)return 0;
192 logx/=3;
193 _x[0]=_a^_gf->exp[logx];
194 _x[1]=_a^_gf->exp[logx+255/3];
195 _x[2]=_a^_x[0]^_x[1];
196 return 3;
197 }
198 logd2=_gf->log[d2];
199 logd=logd2+(255&-(logd2&1))>>1;
200 k=rs_gdiv(_gf,k,_gf->exp[logd+logd2]);
201 /*Substitute y=w+1/w and z=w**3 to get z**2 + k*z + 1 == 0.*/
202 nroots=rs_quadratic_solve(_gf,k,1,_x);
203 if(nroots<1){
204 /*The Reed-Solomon code is only valid if we can find 3 distinct roots in
205 GF(2**8), so if we know there's only one, we don't actually need to find
206 it.
207 Note that we're also called by the quartic solver, but if we contain a
208 non-trivial irreducible factor, than so does the original
209 quartic~\cite{LW72}, and failing to return a root here actually saves us
210 some work there, also.*/
211 return 0;
212 }
213 /*Recover w from z.*/
214 logw=_gf->log[_x[0]];
215 if(logw){
216 if(logw%3!=0)return 0;
217 logw/=3;
218 /*Recover x from w.*/
219 _x[0]=_gf->exp[_gf->log[_gf->exp[logw]^_gf->exp[255-logw]]+logd]^_a;
220 logw+=255/3;
221 _x[1]=_gf->exp[_gf->log[_gf->exp[logw]^_gf->exp[255-logw]]+logd]^_a;
222 _x[2]=_x[0]^_x[1]^_a;
223 return 3;
224 }
225 else{
226 _x[0]=_a;
227 /*In this case _x[1] is a double root, so we know the Reed-Solomon code is
228 invalid.
229 Note that we still have to return at least one root, because if we're
230 being called by the quartic solver, the quartic might still have 4
231 distinct roots.
232 But we don't need more than one root, so we can avoid computing the
233 expensive one.*/
234 /*_x[1]=_gf->exp[_gf->log[_gf->exp[255/3]^_gf->exp[2*(255/3)]]+logd]^_a;*/
235 return 1;
236 }
237 }
238
239 /*Solve a quartic equation x**4 + _a*x**3 + _b*x**2 + _c*x + _d in GF(2**8) by
240 decomposing it into the cases given by~\cite{LW72}.
241 Returns the number of distinct roots.
242 @ARTICLE{LW72,
243 author="Philip A. Leonard and Kenneth S. Williams",
244 title="Quartics over $GF(2^n)$",
245 journal="Proceedings of the American Mathematical Society",
246 volume=36,
247 number=2,
248 pages="347--450",
249 month=Dec,
250 year=1972
251 }*/
rs_quartic_solve(const rs_gf256 * _gf,unsigned _a,unsigned _b,unsigned _c,unsigned _d,unsigned char _x[3])252 static int rs_quartic_solve(const rs_gf256 *_gf,
253 unsigned _a,unsigned _b,unsigned _c,unsigned _d,unsigned char _x[3]){
254 unsigned r;
255 unsigned s;
256 unsigned t;
257 unsigned b;
258 int nroots;
259 int i;
260 /*If _d is zero, factor out the 0 root.*/
261 if(!_d){
262 nroots=rs_cubic_solve(_gf,_a,_b,_c,_x);
263 if(_c)_x[nroots++]=0;
264 return nroots;
265 }
266 if(_a){
267 unsigned loga;
268 /*Substitute x=(1/y) + sqrt(_c/_a) to eliminate the cubic term.*/
269 loga=_gf->log[_a];
270 r=rs_hgmul(_gf,_c,255-loga);
271 s=rs_gsqrt(_gf,r);
272 t=_d^rs_gmul(_gf,_b,r)^rs_gmul(_gf,r,r);
273 if(t){
274 unsigned logti;
275 logti=255-_gf->log[t];
276 /*The result is still quartic, but with no cubic term.*/
277 nroots=rs_quartic_solve(_gf,0,rs_hgmul(_gf,_b^rs_hgmul(_gf,s,loga),logti),
278 _gf->exp[loga+logti],_gf->exp[logti],_x);
279 for(i=0;i<nroots;i++)_x[i]=_gf->exp[255-_gf->log[_x[i]]]^s;
280 }
281 else{
282 /*s must be a root~\cite{LW72}, and is in fact a double-root~\cite{CCO69}.
283 Thus we're left with only a quadratic to solve.
284 @ARTICLE{CCO69,
285 author="Robert T. Chien and B. D. Cunningham and I. B. Oldham",
286 title="Hybrid Methods for Finding Roots of a Polynomial---With
287 Applications to {BCH} Decoding",
288 journal="{IEEE} Transactions on Information Theory",
289 volume=15,
290 number=2,
291 pages="329--335",
292 month=Mar,
293 year=1969
294 }*/
295 nroots=rs_quadratic_solve(_gf,_a,_b^r,_x);
296 /*s may be a triple root if s=_b/_a, but not quadruple, since _a!=0.*/
297 if(nroots!=2||_x[0]!=s&&_x[1]!=s)_x[nroots++]=s;
298 }
299 return nroots;
300 }
301 /*If there are no odd powers, it's really just a quadratic in disguise.*/
302 if(!_c)return rs_quadratic_solve(_gf,rs_gsqrt(_gf,_b),rs_gsqrt(_gf,_d),_x);
303 /*Factor into (x**2 + r*x + s)*(x**2 + r*x + t) by solving for r, which can
304 be shown to satisfy r**3 + _b*r + _c == 0.*/
305 nroots=rs_cubic_solve(_gf,0,_b,_c,_x);
306 if(nroots<1){
307 /*The Reed-Solomon code is only valid if we can find 4 distinct roots in
308 GF(2**8).
309 If the cubic does not factor into 3 (possibly duplicate) roots, then we
310 know that the quartic must have a non-trivial irreducible factor.*/
311 return 0;
312 }
313 r=_x[0];
314 /*Now solve for s and t.*/
315 b=rs_gdiv(_gf,_c,r);
316 nroots=rs_quadratic_solve(_gf,b,_d,_x);
317 if(nroots<2)return 0;
318 s=_x[0];
319 t=_x[1];
320 /*_c=r*(s^t) was non-zero, so s and t must be distinct.
321 But if z is a root of z**2 ^ r*z ^ s, then so is (z^r), and s = z*(z^r).
322 Hence if z is also a root of z**2 ^ r*z ^ t, then t = s, a contradiction.
323 Thus all four roots are distinct, if they exist.*/
324 nroots=rs_quadratic_solve(_gf,r,s,_x);
325 return nroots+rs_quadratic_solve(_gf,r,t,_x+nroots);
326 }
327
328 /*Polynomial arithmetic with coefficients in GF(2**8).*/
329
rs_poly_zero(unsigned char * _p,int _dp1)330 static void rs_poly_zero(unsigned char *_p,int _dp1){
331 memset(_p,0,_dp1*sizeof(*_p));
332 }
333
rs_poly_copy(unsigned char * _p,const unsigned char * _q,int _dp1)334 static void rs_poly_copy(unsigned char *_p,const unsigned char *_q,int _dp1){
335 memcpy(_p,_q,_dp1*sizeof(*_p));
336 }
337
338 /*Multiply the polynomial by the free variable, x (shift the coefficients).
339 The number of coefficients, _dp1, must be non-zero.*/
rs_poly_mul_x(unsigned char * _p,const unsigned char * _q,int _dp1)340 static void rs_poly_mul_x(unsigned char *_p,const unsigned char *_q,int _dp1){
341 memmove(_p+1,_q,(_dp1-1)*sizeof(*_p));
342 _p[0]=0;
343 }
344
345 /*Divide the polynomial by the free variable, x (shift the coefficients).
346 The number of coefficients, _dp1, must be non-zero.*/
rs_poly_div_x(unsigned char * _p,const unsigned char * _q,int _dp1)347 static void rs_poly_div_x(unsigned char *_p,const unsigned char *_q,int _dp1){
348 memmove(_p,_q+1,(_dp1-1)*sizeof(*_p));
349 _p[_dp1-1]=0;
350 }
351
352 /*Compute the first (d+1) coefficients of the product of a degree e and a
353 degree f polynomial.*/
rs_poly_mult(const rs_gf256 * _gf,unsigned char * _p,int _dp1,const unsigned char * _q,int _ep1,const unsigned char * _r,int _fp1)354 static void rs_poly_mult(const rs_gf256 *_gf,unsigned char *_p,int _dp1,
355 const unsigned char *_q,int _ep1,const unsigned char *_r,int _fp1){
356 int m;
357 int i;
358 rs_poly_zero(_p,_dp1);
359 m=_ep1<_dp1?_ep1:_dp1;
360 for(i=0;i<m;i++)if(_q[i]!=0){
361 unsigned logqi;
362 int n;
363 int j;
364 n=_dp1-i<_fp1?_dp1-i:_fp1;
365 logqi=_gf->log[_q[i]];
366 for(j=0;j<n;j++)_p[i+j]^=rs_hgmul(_gf,_r[j],logqi);
367 }
368 }
369
370 /*Decoding.*/
371
372 /*Computes the syndrome of a codeword.*/
rs_calc_syndrome(const rs_gf256 * _gf,int _m0,unsigned char * _s,int _npar,const unsigned char * _data,int _ndata)373 static void rs_calc_syndrome(const rs_gf256 *_gf,int _m0,
374 unsigned char *_s,int _npar,const unsigned char *_data,int _ndata){
375 int i;
376 int j;
377 for(j=0;j<_npar;j++){
378 unsigned alphaj;
379 unsigned sj;
380 sj=0;
381 alphaj=_gf->log[_gf->exp[j+_m0]];
382 for(i=0;i<_ndata;i++)sj=_data[i]^rs_hgmul(_gf,sj,alphaj);
383 _s[j]=sj;
384 }
385 }
386
387 /*Berlekamp-Peterson and Berlekamp-Massey Algorithms for error-location,
388 modified to handle known erasures, from \cite{CC81}, p. 205.
389 This finds the coefficients of the error locator polynomial.
390 The roots are then found by looking for the values of alpha**n where
391 evaluating the polynomial yields zero.
392 Error correction is done using the error-evaluator equation on p. 207.
393 @BOOK{CC81,
394 author="George C. Clark, Jr and J. Bibb Cain",
395 title="Error-Correction Coding for Digital Communications",
396 series="Applications of Communications Theory",
397 publisher="Springer",
398 address="New York, NY",
399 month=Jun,
400 year=1981
401 }*/
402
403 /*Initialize lambda to the product of (1-x*alpha**e[i]) for erasure locations
404 e[i].
405 Note that the user passes in array indices counting from the beginning of the
406 data, while our polynomial indexes starting from the end, so
407 e[i]=(_ndata-1)-_erasures[i].*/
rs_init_lambda(const rs_gf256 * _gf,unsigned char * _lambda,int _npar,const unsigned char * _erasures,int _nerasures,int _ndata)408 static void rs_init_lambda(const rs_gf256 *_gf,unsigned char *_lambda,int _npar,
409 const unsigned char *_erasures,int _nerasures,int _ndata){
410 int i;
411 int j;
412 rs_poly_zero(_lambda,(_npar<4?4:_npar)+1);
413 _lambda[0]=1;
414 for(i=0;i<_nerasures;i++)for(j=i+1;j>0;j--){
415 _lambda[j]^=rs_hgmul(_gf,_lambda[j-1],_ndata-1-_erasures[i]);
416 }
417 }
418
419 /*From \cite{CC81}, p. 216.
420 Returns the number of errors detected (degree of _lambda).*/
rs_modified_berlekamp_massey(const rs_gf256 * _gf,unsigned char * _lambda,const unsigned char * _s,unsigned char * _omega,int _npar,const unsigned char * _erasures,int _nerasures,int _ndata)421 static int rs_modified_berlekamp_massey(const rs_gf256 *_gf,
422 unsigned char *_lambda,const unsigned char *_s,unsigned char *_omega,int _npar,
423 const unsigned char *_erasures,int _nerasures,int _ndata){
424 unsigned char tt[256];
425 int n;
426 int l;
427 int k;
428 int i;
429 /*Initialize _lambda, the error locator-polynomial, with the location of
430 known erasures.*/
431 rs_init_lambda(_gf,_lambda,_npar,_erasures,_nerasures,_ndata);
432 rs_poly_copy(tt,_lambda,_npar+1);
433 l=_nerasures;
434 k=0;
435 for(n=_nerasures+1;n<=_npar;n++){
436 unsigned d;
437 rs_poly_mul_x(tt,tt,n-k+1);
438 d=0;
439 for(i=0;i<=l;i++)d^=rs_gmul(_gf,_lambda[i],_s[n-1-i]);
440 if(d!=0){
441 unsigned logd;
442 logd=_gf->log[d];
443 if(l<n-k){
444 int t;
445 for(i=0;i<=n-k;i++){
446 unsigned tti;
447 tti=tt[i];
448 tt[i]=rs_hgmul(_gf,_lambda[i],255-logd);
449 _lambda[i]=_lambda[i]^rs_hgmul(_gf,tti,logd);
450 }
451 t=n-k;
452 k=n-l;
453 l=t;
454 }
455 else for(i=0;i<=l;i++)_lambda[i]=_lambda[i]^rs_hgmul(_gf,tt[i],logd);
456 }
457 }
458 rs_poly_mult(_gf,_omega,_npar,_lambda,l+1,_s,_npar);
459 return l;
460 }
461
462 /*Finds all the roots of an error-locator polynomial _lambda by evaluating it
463 at successive values of alpha, and returns the positions of the associated
464 errors in _epos.
465 Returns the number of valid roots identified.*/
rs_find_roots(const rs_gf256 * _gf,unsigned char * _epos,const unsigned char * _lambda,int _nerrors,int _ndata)466 static int rs_find_roots(const rs_gf256 *_gf,unsigned char *_epos,
467 const unsigned char *_lambda,int _nerrors,int _ndata){
468 unsigned alpha;
469 int nroots;
470 int i;
471 nroots=0;
472 if(_nerrors<=4){
473 /*Explicit solutions for higher degrees are possible.
474 Chien uses large lookup tables to solve quintics, and Truong et al. give
475 special algorithms for degree up through 11, though they use exhaustive
476 search (with reduced complexity) for some portions.
477 Quartics are good enough for reading CDs, and represent a reasonable code
478 complexity trade-off without requiring any extra tables.
479 Note that _lambda[0] is always 1.*/
480 _nerrors=rs_quartic_solve(_gf,_lambda[1],_lambda[2],_lambda[3],_lambda[4],
481 _epos);
482 for(i=0;i<_nerrors;i++)if(_epos[i]){
483 alpha=_gf->log[_epos[i]];
484 if((int)alpha<_ndata)_epos[nroots++]=alpha;
485 }
486 return nroots;
487 }
488 else for(alpha=0;(int)alpha<_ndata;alpha++){
489 unsigned alphai;
490 unsigned sum;
491 sum=0;
492 alphai=0;
493 for(i=0;i<=_nerrors;i++){
494 sum^=rs_hgmul(_gf,_lambda[_nerrors-i],alphai);
495 alphai=_gf->log[_gf->exp[alphai+alpha]];
496 }
497 if(!sum)_epos[nroots++]=alpha;
498 }
499 return nroots;
500 }
501
502 /*Corrects a codeword with _ndata<256 bytes, of which the last _npar are parity
503 bytes.
504 Known locations of errors can be passed in the _erasures array.
505 Twice as many (up to _npar) errors with a known location can be corrected
506 compared to errors with an unknown location.
507 Returns the number of errors corrected if successful, or a negative number if
508 the message could not be corrected because too many errors were detected.*/
rs_correct(const rs_gf256 * _gf,int _m0,unsigned char * _data,int _ndata,int _npar,const unsigned char * _erasures,int _nerasures)509 int rs_correct(const rs_gf256 *_gf,int _m0,unsigned char *_data,int _ndata,
510 int _npar,const unsigned char *_erasures,int _nerasures){
511 /*lambda must have storage for at least five entries to avoid special cases
512 in the low-degree polynomial solver.*/
513 unsigned char lambda[256];
514 unsigned char omega[256];
515 unsigned char epos[256];
516 unsigned char s[256];
517 int i;
518 /*If we already have too many erasures, we can't possibly succeed.*/
519 if(_nerasures>_npar)return -1;
520 /*Compute the syndrome values.*/
521 rs_calc_syndrome(_gf,_m0,s,_npar,_data,_ndata);
522 /*Check for a non-zero value.*/
523 for(i=0;i<_npar;i++)if(s[i]){
524 int nerrors;
525 int j;
526 /*Construct the error locator polynomial.*/
527 nerrors=rs_modified_berlekamp_massey(_gf,lambda,s,omega,_npar,
528 _erasures,_nerasures,_ndata);
529 /*If we can't locate any errors, we can't force the syndrome values to
530 zero, and must have a decoding error.
531 Conversely, if we have too many errors, there's no reason to even attempt
532 the root search.*/
533 if(nerrors<=0||nerrors-_nerasures>_npar-_nerasures>>1)return -1;
534 /*Compute the locations of the errors.
535 If they are not all distinct, or some of them were outside the valid
536 range for our block size, we have a decoding error.*/
537 if(rs_find_roots(_gf,epos,lambda,nerrors,_ndata)<nerrors)return -1;
538 /*Now compute the error magnitudes.*/
539 for(i=0;i<nerrors;i++){
540 unsigned a;
541 unsigned b;
542 unsigned alpha;
543 unsigned alphan1;
544 unsigned alphan2;
545 unsigned alphanj;
546 alpha=epos[i];
547 /*Evaluate omega at alpha**-1.*/
548 a=0;
549 alphan1=255-alpha;
550 alphanj=0;
551 for(j=0;j<_npar;j++){
552 a^=rs_hgmul(_gf,omega[j],alphanj);
553 alphanj=_gf->log[_gf->exp[alphanj+alphan1]];
554 }
555 /*Evaluate the derivative of lambda at alpha**-1
556 All the odd powers vanish.*/
557 b=0;
558 alphan2=_gf->log[_gf->exp[alphan1<<1]];
559 alphanj=alphan1+_m0*alpha%255;
560 for(j=1;j<=_npar;j+=2){
561 b^=rs_hgmul(_gf,lambda[j],alphanj);
562 alphanj=_gf->log[_gf->exp[alphanj+alphan2]];
563 }
564 /*Apply the correction.*/
565 _data[_ndata-1-alpha]^=rs_gdiv(_gf,a,b);
566 }
567 return nerrors;
568 }
569 return 0;
570 }
571
572 /*Encoding.*/
573
574 /*Create an _npar-coefficient generator polynomial for a Reed-Solomon code
575 with _npar<256 parity bytes.*/
rs_compute_genpoly(const rs_gf256 * _gf,int _m0,unsigned char * _genpoly,int _npar)576 void rs_compute_genpoly(const rs_gf256 *_gf,int _m0,
577 unsigned char *_genpoly,int _npar){
578 int i;
579 if(_npar<=0)return;
580 rs_poly_zero(_genpoly,_npar);
581 _genpoly[0]=1;
582 /*Multiply by (x+alpha^i) for i = 1 ... _ndata.*/
583 for(i=0;i<_npar;i++){
584 unsigned alphai;
585 int n;
586 int j;
587 n=i+1<_npar-1?i+1:_npar-1;
588 alphai=_gf->log[_gf->exp[_m0+i]];
589 for(j=n;j>0;j--)_genpoly[j]=_genpoly[j-1]^rs_hgmul(_gf,_genpoly[j],alphai);
590 _genpoly[0]=rs_hgmul(_gf,_genpoly[0],alphai);
591 }
592 }
593
594 /*Adds _npar<=_ndata parity bytes to an _ndata-_npar byte message.
595 _data must contain room for _ndata<256 bytes.*/
rs_encode(const rs_gf256 * _gf,unsigned char * _data,int _ndata,const unsigned char * _genpoly,int _npar)596 void rs_encode(const rs_gf256 *_gf,unsigned char *_data,int _ndata,
597 const unsigned char *_genpoly,int _npar){
598 unsigned char *lfsr;
599 unsigned d;
600 int i;
601 int j;
602 if(_npar<=0)return;
603 lfsr=_data+_ndata-_npar;
604 rs_poly_zero(lfsr,_npar);
605 for(i=0;i<_ndata-_npar;i++){
606 d=_data[i]^lfsr[0];
607 if(d){
608 unsigned logd;
609 logd=_gf->log[d];
610 for(j=0;j<_npar-1;j++){
611 lfsr[j]=lfsr[j+1]^rs_hgmul(_gf,_genpoly[_npar-1-j],logd);
612 }
613 lfsr[_npar-1]=rs_hgmul(_gf,_genpoly[0],logd);
614 }
615 else rs_poly_div_x(lfsr,lfsr,_npar);
616 }
617 }
618
619 #if defined(RS_TEST_ENC)
620 #include <stdio.h>
621 #include <stdlib.h>
622
main(void)623 int main(void){
624 rs_gf256 gf;
625 int k;
626 rs_gf256_init(&gf,QR_PPOLY);
627 srand(0);
628 for(k=0;k<64*1024;k++){
629 unsigned char genpoly[256];
630 unsigned char data[256];
631 unsigned char epos[256];
632 int ndata;
633 int npar;
634 int nerrors;
635 int i;
636 ndata=rand()&0xFF;
637 npar=ndata>0?rand()%ndata:0;
638 for(i=0;i<ndata-npar;i++)data[i]=rand()&0xFF;
639 rs_compute_genpoly(&gf,QR_M0,genpoly,npar);
640 rs_encode(&gf,QR_M0,data,ndata,genpoly,npar);
641 /*Write a clean version of the codeword.*/
642 printf("%i %i",ndata,npar);
643 for(i=0;i<ndata;i++)printf(" %i",data[i]);
644 printf(" 0\n");
645 /*Write the correct output to compare the decoder against.*/
646 fprintf(stderr,"Success!\n",nerrors);
647 for(i=0;i<ndata;i++)fprintf(stderr,"%i%s",data[i],i+1<ndata?" ":"\n");
648 if(npar>0){
649 /*Corrupt it.*/
650 nerrors=rand()%(npar+1);
651 if(nerrors>0){
652 /*This test is not quite correct: there could be so many errors it
653 comes within (npar>>1) errors of another valid codeword.
654 I don't know a simple way to test for that without trying to decode
655 the corrupt codeword, though, which is the very code we're trying to
656 test.*/
657 if(nerrors<=npar>>1){
658 fprintf(stderr,"Success!\n",nerrors);
659 for(i=0;i<ndata;i++){
660 fprintf(stderr,"%i%s",data[i],i+1<ndata?" ":"\n");
661 }
662 }
663 else fprintf(stderr,"Failure.\n");
664 fprintf(stderr,"Success!\n",nerrors);
665 for(i=0;i<ndata;i++)fprintf(stderr,"%i%s",data[i],i+1<ndata?" ":"\n");
666 for(i=0;i<ndata;i++)epos[i]=i;
667 for(i=0;i<nerrors;i++){
668 unsigned char e;
669 int ei;
670 ei=rand()%(ndata-i)+i;
671 e=epos[ei];
672 epos[ei]=epos[i];
673 epos[i]=e;
674 data[e]^=rand()%255+1;
675 }
676 /*First with no erasure locations.*/
677 printf("%i %i",ndata,npar);
678 for(i=0;i<ndata;i++)printf(" %i",data[i]);
679 printf(" 0\n");
680 /*Now with erasure locations.*/
681 printf("%i %i",ndata,npar);
682 for(i=0;i<ndata;i++)printf(" %i",data[i]);
683 printf(" %i",nerrors);
684 for(i=0;i<nerrors;i++)printf(" %i",epos[i]);
685 printf("\n");
686 }
687 }
688 }
689 return 0;
690 }
691 #endif
692
693 #if defined(RS_TEST_DEC)
694 #include <stdio.h>
695
main(void)696 int main(void){
697 rs_gf256 gf;
698 rs_gf256_init(&gf,QR_PPOLY);
699 for(;;){
700 unsigned char data[255];
701 unsigned char erasures[255];
702 int idata[255];
703 int ierasures[255];
704 int ndata;
705 int npar;
706 int nerasures;
707 int nerrors;
708 int i;
709 if(scanf("%i",&ndata)<1||ndata<0||ndata>255||
710 scanf("%i",&npar)<1||npar<0||npar>ndata)break;
711 for(i=0;i<ndata;i++){
712 if(scanf("%i",idata+i)<1||idata[i]<0||idata[i]>255)break;
713 data[i]=idata[i];
714 }
715 if(i<ndata)break;
716 if(scanf("%i",&nerasures)<1||nerasures<0||nerasures>ndata)break;
717 for(i=0;i<nerasures;i++){
718 if(scanf("%i",ierasures+i)<1||ierasures[i]<0||ierasures[i]>=ndata)break;
719 erasures[i]=ierasures[i];
720 }
721 nerrors=rs_correct(&gf,QR_M0,data,ndata,npar,erasures,nerasures);
722 if(nerrors>=0){
723 unsigned char data2[255];
724 unsigned char genpoly[255];
725 for(i=0;i<ndata-npar;i++)data2[i]=data[i];
726 rs_compute_genpoly(&gf,QR_M0,genpoly,npar);
727 rs_encode(&gf,QR_M0,data2,ndata,genpoly,npar);
728 for(i=ndata-npar;i<ndata;i++)if(data[i]!=data2[i]){
729 printf("Abject failure! %i!=%i\n",data[i],data2[i]);
730 }
731 printf("Success!\n",nerrors);
732 for(i=0;i<ndata;i++)printf("%i%s",data[i],i+1<ndata?" ":"\n");
733 }
734 else printf("Failure.\n");
735 }
736 return 0;
737 }
738 #endif
739
740 #if defined(RS_TEST_ROOTS)
741 #include <stdio.h>
742
743 /*Exhaustively test the root finder.*/
main(void)744 int main(void){
745 rs_gf256 gf;
746 int a;
747 int b;
748 int c;
749 int d;
750 rs_gf256_init(&gf,QR_PPOLY);
751 for(a=0;a<256;a++)for(b=0;b<256;b++)for(c=0;c<256;c++)for(d=0;d<256;d++){
752 unsigned char x[4];
753 unsigned char r[4];
754 unsigned x2;
755 unsigned e[5];
756 int nroots;
757 int mroots;
758 int i;
759 int j;
760 nroots=rs_quartic_solve(&gf,a,b,c,d,x);
761 for(i=0;i<nroots;i++){
762 x2=rs_gmul(&gf,x[i],x[i]);
763 e[0]=rs_gmul(&gf,x2,x2)^rs_gmul(&gf,a,rs_gmul(&gf,x[i],x2))^
764 rs_gmul(&gf,b,x2)^rs_gmul(&gf,c,x[i])^d;
765 if(e[0]){
766 printf("Invalid root: (0x%02X)**4 ^ 0x%02X*(0x%02X)**3 ^ "
767 "0x%02X*(0x%02X)**2 ^ 0x%02X(0x%02X) ^ 0x%02X = 0x%02X\n",
768 x[i],a,x[i],b,x[i],c,x[i],d,e[0]);
769 }
770 for(j=0;j<i;j++)if(x[i]==x[j]){
771 printf("Repeated root %i=%i: (0x%02X)**4 ^ 0x%02X*(0x%02X)**3 ^ "
772 "0x%02X*(0x%02X)**2 ^ 0x%02X(0x%02X) ^ 0x%02X = 0x%02X\n",
773 i,j,x[i],a,x[i],b,x[i],c,x[i],d,e[0]);
774 }
775 }
776 mroots=0;
777 for(j=1;j<256;j++){
778 int logx;
779 int logx2;
780 logx=gf.log[j];
781 logx2=gf.log[gf.exp[logx<<1]];
782 e[mroots]=d^rs_hgmul(&gf,c,logx)^rs_hgmul(&gf,b,logx2)^
783 rs_hgmul(&gf,a,gf.log[gf.exp[logx+logx2]])^gf.exp[logx2<<1];
784 if(!e[mroots])r[mroots++]=j;
785 }
786 /*We only care about missing roots if the quartic had 4 distinct, non-zero
787 roots.*/
788 if(mroots==4)for(j=0;j<mroots;j++){
789 for(i=0;i<nroots;i++)if(x[i]==r[j])break;
790 if(i>=nroots){
791 printf("Missing root: (0x%02X)**4 ^ 0x%02X*(0x%02X)**3 ^ "
792 "0x%02X*(0x%02X)**2 ^ 0x%02X(0x%02X) ^ 0x%02X = 0x%02X\n",
793 r[j],a,r[j],b,r[j],c,r[j],d,e[j]);
794 }
795 }
796 }
797 return 0;
798 }
799 #endif
800