ZBar bar code reader . http://zbar.sourceforge.net/ ZBar is licensed under the GNU LGPL 2.1 to enable development of both open source and commercial projects.

Dependents:   GR-PEACH_Camera_in_barcode levkov_ov7670

Embed: (wiki syntax)

« Back to documentation index

Show/hide line numbers rs.c Source File

rs.c

00001 /*Copyright (C) 1991-1995  Henry Minsky (hqm@ua.com, hqm@ai.mit.edu)
00002   Copyright (C) 2008-2009  Timothy B. Terriberry (tterribe@xiph.org)
00003   You can redistribute this library and/or modify it under the terms of the
00004    GNU Lesser General Public License as published by the Free Software
00005    Foundation; either version 2.1 of the License, or (at your option) any later
00006    version.*/
00007 #include <stdlib.h>
00008 #include <string.h>
00009 #include "rs.h"
00010 
00011 /*Reed-Solomon encoder and decoder.
00012   Original implementation (C) Henry Minsky (hqm@ua.com, hqm@ai.mit.edu),
00013    Universal Access 1991-1995.
00014   Updates by Timothy B. Terriberry (C) 2008-2009:
00015    - Properly reject codes when error-locator polynomial has repeated roots or
00016       non-trivial irreducible factors.
00017    - Removed the hard-coded parity size and performed general API cleanup.
00018    - Allow multiple representations of GF(2**8), since different standards use
00019       different irreducible polynomials.
00020    - Allow different starting indices for the generator polynomial, since
00021       different standards use different values.
00022    - Greatly reduced the computation by eliminating unnecessary operations.
00023    - Explicitly solve for the roots of low-degree polynomials instead of using
00024       an exhaustive search.
00025      This is another major speed boost when there are few errors.*/
00026 
00027 
00028 /*Galois Field arithmetic in GF(2**8).*/
00029 
00030 void rs_gf256_init(rs_gf256 *_gf,unsigned _ppoly){
00031   unsigned p;
00032   int      i;
00033   /*Initialize the table of powers of a primtive root, alpha=0x02.*/
00034   p=1;
00035   for(i=0;i<256;i++){
00036     _gf->exp[i]=_gf->exp[i+255]=p;
00037     p=((p<<1)^(-(p>>7)&_ppoly))&0xFF;
00038   }
00039   /*Invert the table to recover the logs.*/
00040   for(i=0;i<255;i++)_gf->log[_gf->exp[i]]=i;
00041   /*Note that we rely on the fact that _gf->log[0]=0 below.*/
00042   _gf->log[0]=0;
00043 }
00044 
00045 /*Multiplication in GF(2**8) using logarithms.*/
00046 static unsigned rs_gmul(const rs_gf256 *_gf,unsigned _a,unsigned _b){
00047   return _a==0||_b==0?0:_gf->exp[_gf->log[_a]+_gf->log[_b]];
00048 }
00049 
00050 /*Division in GF(2**8) using logarithms.
00051   The result of division by zero is undefined.*/
00052 static unsigned rs_gdiv(const rs_gf256 *_gf,unsigned _a,unsigned _b){
00053   return _a==0?0:_gf->exp[_gf->log[_a]+255-_gf->log[_b]];
00054 }
00055 
00056 /*Multiplication in GF(2**8) when one of the numbers is known to be non-zero
00057    (proven by representing it by its logarithm).*/
00058 static unsigned rs_hgmul(const rs_gf256 *_gf,unsigned _a,unsigned _logb){
00059   return _a==0?0:_gf->exp[_gf->log[_a]+_logb];
00060 }
00061 
00062 /*Square root in GF(2**8) using logarithms.*/
00063 static unsigned rs_gsqrt(const rs_gf256 *_gf,unsigned _a){
00064   unsigned loga;
00065   if(!_a)return 0;
00066   loga=_gf->log[_a];
00067   return _gf->exp[loga+(255&-(loga&1))>>1];
00068 }
00069 
00070 /*Polynomial root finding in GF(2**8).
00071   Each routine returns a list of the distinct roots (i.e., with duplicates
00072    removed).*/
00073 
00074 /*Solve a quadratic equation x**2 + _b*x + _c in GF(2**8) using the method
00075    of~\cite{Wal99}.
00076   Returns the number of distinct roots.
00077   ARTICLE{Wal99,
00078     author="C. Wayne Walker",
00079     title="New Formulas for Solving Quadratic Equations over Certain Finite
00080      Fields",
00081     journal="{IEEE} Transactions on Information Theory",
00082     volume=45,
00083     number=1,
00084     pages="283--284",
00085     month=Jan,
00086     year=1999
00087   }*/
00088 static int rs_quadratic_solve(const rs_gf256 *_gf,unsigned _b,unsigned _c,
00089  unsigned char _x[2]){
00090   unsigned b;
00091   unsigned logb;
00092   unsigned logb2;
00093   unsigned logb4;
00094   unsigned logb8;
00095   unsigned logb12;
00096   unsigned logb14;
00097   unsigned logc;
00098   unsigned logc2;
00099   unsigned logc4;
00100   unsigned c8;
00101   unsigned g3;
00102   unsigned z3;
00103   unsigned l3;
00104   unsigned c0;
00105   unsigned g2;
00106   unsigned l2;
00107   unsigned z2;
00108   int      inc;
00109   /*If _b is zero, all we need is a square root.*/
00110   if(!_b){
00111     _x[0]=rs_gsqrt(_gf,_c);
00112     return 1;
00113   }
00114   /*If _c is zero, 0 and _b are the roots.*/
00115   if(!_c){
00116     _x[0]=0;
00117     _x[1]=_b;
00118     return 2;
00119   }
00120   logb=_gf->log[_b];
00121   logc=_gf->log[_c];
00122   /*If _b lies in GF(2**4), scale x to move it out.*/
00123   inc=logb%(255/15)==0;
00124   if(inc){
00125     b=_gf->exp[logb+254];
00126     logb=_gf->log[b];
00127     _c=_gf->exp[logc+253];
00128     logc=_gf->log[_c];
00129   }
00130   else b=_b;
00131   logb2=_gf->log[_gf->exp[logb<<1]];
00132   logb4=_gf->log[_gf->exp[logb2<<1]];
00133   logb8=_gf->log[_gf->exp[logb4<<1]];
00134   logb12=_gf->log[_gf->exp[logb4+logb8]];
00135   logb14=_gf->log[_gf->exp[logb2+logb12]];
00136   logc2=_gf->log[_gf->exp[logc<<1]];
00137   logc4=_gf->log[_gf->exp[logc2<<1]];
00138   c8=_gf->exp[logc4<<1];
00139   g3=rs_hgmul(_gf,
00140    _gf->exp[logb14+logc]^_gf->exp[logb12+logc2]^_gf->exp[logb8+logc4]^c8,logb);
00141   /*If g3 doesn't lie in GF(2**4), then our roots lie in an extension field.
00142     Note that we rely on the fact that _gf->log[0]==0 here.*/
00143   if(_gf->log[g3]%(255/15)!=0)return 0;
00144   /*Construct the corresponding quadratic in GF(2**4):
00145     x**2 + x/alpha**(255/15) + l3/alpha**(2*(255/15))*/
00146   z3=rs_gdiv(_gf,g3,_gf->exp[logb8<<1]^b);
00147   l3=rs_hgmul(_gf,rs_gmul(_gf,z3,z3)^rs_hgmul(_gf,z3,logb)^_c,255-logb2);
00148   c0=rs_hgmul(_gf,l3,255-2*(255/15));
00149   /*Construct the corresponding quadratic in GF(2**2):
00150     x**2 + x/alpha**(255/3) + l2/alpha**(2*(255/3))*/
00151   g2=rs_hgmul(_gf,
00152    rs_hgmul(_gf,c0,255-2*(255/15))^rs_gmul(_gf,c0,c0),255-255/15);
00153   z2=rs_gdiv(_gf,g2,_gf->exp[255-(255/15)*4]^_gf->exp[255-(255/15)]);
00154   l2=rs_hgmul(_gf,
00155    rs_gmul(_gf,z2,z2)^rs_hgmul(_gf,z2,255-(255/15))^c0,2*(255/15));
00156   /*Back substitute to the solution in the original field.*/
00157   _x[0]=_gf->exp[_gf->log[z3^rs_hgmul(_gf,
00158    rs_hgmul(_gf,l2,255/3)^rs_hgmul(_gf,z2,255/15),logb)]+inc];
00159   _x[1]=_x[0]^_b;
00160   return 2;
00161 }
00162 
00163 /*Solve a cubic equation x**3 + _a*x**2 + _b*x + _c in GF(2**8).
00164   Returns the number of distinct roots.*/
00165 static int rs_cubic_solve(const rs_gf256 *_gf,
00166  unsigned _a,unsigned _b,unsigned _c,unsigned char _x[3]){
00167   unsigned k;
00168   unsigned logd;
00169   unsigned d2;
00170   unsigned logd2;
00171   unsigned logw;
00172   int      nroots;
00173   /*If _c is zero, factor out the 0 root.*/
00174   if(!_c){
00175     nroots=rs_quadratic_solve(_gf,_a,_b,_x);
00176     if(_b)_x[nroots++]=0;
00177     return nroots;
00178   }
00179   /*Substitute x=_a+y*sqrt(_a**2+_b) to get y**3 + y + k == 0,
00180      k = (_a*_b+c)/(_a**2+b)**(3/2).*/
00181   k=rs_gmul(_gf,_a,_b)^_c;
00182   d2=rs_gmul(_gf,_a,_a)^_b;
00183   if(!d2){
00184     int logx;
00185     if(!k){
00186       /*We have a triple root.*/
00187       _x[0]=_a;
00188       return 1;
00189     }
00190     logx=_gf->log[k];
00191     if(logx%3!=0)return 0;
00192     logx/=3;
00193     _x[0]=_a^_gf->exp[logx];
00194     _x[1]=_a^_gf->exp[logx+255/3];
00195     _x[2]=_a^_x[0]^_x[1];
00196     return 3;
00197   }
00198   logd2=_gf->log[d2];
00199   logd=logd2+(255&-(logd2&1))>>1;
00200   k=rs_gdiv(_gf,k,_gf->exp[logd+logd2]);
00201   /*Substitute y=w+1/w and z=w**3 to get z**2 + k*z + 1 == 0.*/
00202   nroots=rs_quadratic_solve(_gf,k,1,_x);
00203   if(nroots<1){
00204     /*The Reed-Solomon code is only valid if we can find 3 distinct roots in
00205        GF(2**8), so if we know there's only one, we don't actually need to find
00206        it.
00207       Note that we're also called by the quartic solver, but if we contain a
00208        non-trivial irreducible factor, than so does the original
00209        quartic~\cite{LW72}, and failing to return a root here actually saves us
00210        some work there, also.*/
00211     return 0;
00212   }
00213   /*Recover w from z.*/
00214   logw=_gf->log[_x[0]];
00215   if(logw){
00216     if(logw%3!=0)return 0;
00217     logw/=3;
00218     /*Recover x from w.*/
00219     _x[0]=_gf->exp[_gf->log[_gf->exp[logw]^_gf->exp[255-logw]]+logd]^_a;
00220     logw+=255/3;
00221     _x[1]=_gf->exp[_gf->log[_gf->exp[logw]^_gf->exp[255-logw]]+logd]^_a;
00222     _x[2]=_x[0]^_x[1]^_a;
00223     return 3;
00224   }
00225   else{
00226     _x[0]=_a;
00227     /*In this case _x[1] is a double root, so we know the Reed-Solomon code is
00228        invalid.
00229       Note that we still have to return at least one root, because if we're
00230        being called by the quartic solver, the quartic might still have 4
00231        distinct roots.
00232       But we don't need more than one root, so we can avoid computing the
00233        expensive one.*/
00234     /*_x[1]=_gf->exp[_gf->log[_gf->exp[255/3]^_gf->exp[2*(255/3)]]+logd]^_a;*/
00235     return 1;
00236   }
00237 }
00238 
00239 /*Solve a quartic equation x**4 + _a*x**3 + _b*x**2 + _c*x + _d in GF(2**8) by
00240    decomposing it into the cases given by~\cite{LW72}.
00241   Returns the number of distinct roots.
00242   @ARTICLE{LW72,
00243     author="Philip A. Leonard and Kenneth S. Williams",
00244     title="Quartics over $GF(2^n)$",
00245     journal="Proceedings of the American Mathematical Society",
00246     volume=36,
00247     number=2,
00248     pages="347--450",
00249     month=Dec,
00250     year=1972
00251   }*/
00252 static int rs_quartic_solve(const rs_gf256 *_gf,
00253  unsigned _a,unsigned _b,unsigned _c,unsigned _d,unsigned char _x[3]){
00254   unsigned r;
00255   unsigned s;
00256   unsigned t;
00257   unsigned b;
00258   int      nroots;
00259   int      i;
00260   /*If _d is zero, factor out the 0 root.*/
00261   if(!_d){
00262     nroots=rs_cubic_solve(_gf,_a,_b,_c,_x);
00263     if(_c)_x[nroots++]=0;
00264     return nroots;
00265   }
00266   if(_a){
00267     unsigned loga;
00268     /*Substitute x=(1/y) + sqrt(_c/_a) to eliminate the cubic term.*/
00269     loga=_gf->log[_a];
00270     r=rs_hgmul(_gf,_c,255-loga);
00271     s=rs_gsqrt(_gf,r);
00272     t=_d^rs_gmul(_gf,_b,r)^rs_gmul(_gf,r,r);
00273     if(t){
00274       unsigned logti;
00275       logti=255-_gf->log[t];
00276       /*The result is still quartic, but with no cubic term.*/
00277       nroots=rs_quartic_solve(_gf,0,rs_hgmul(_gf,_b^rs_hgmul(_gf,s,loga),logti),
00278        _gf->exp[loga+logti],_gf->exp[logti],_x);
00279       for(i=0;i<nroots;i++)_x[i]=_gf->exp[255-_gf->log[_x[i]]]^s;
00280     }
00281     else{
00282       /*s must be a root~\cite{LW72}, and is in fact a double-root~\cite{CCO69}.
00283         Thus we're left with only a quadratic to solve.
00284         @ARTICLE{CCO69,
00285           author="Robert T. Chien and B. D. Cunningham and I. B. Oldham",
00286           title="Hybrid Methods for Finding Roots of a Polynomial---With
00287            Applications to {BCH} Decoding",
00288           journal="{IEEE} Transactions on Information Theory",
00289           volume=15,
00290           number=2,
00291           pages="329--335",
00292           month=Mar,
00293           year=1969
00294         }*/
00295       nroots=rs_quadratic_solve(_gf,_a,_b^r,_x);
00296       /*s may be a triple root if s=_b/_a, but not quadruple, since _a!=0.*/
00297       if(nroots!=2||_x[0]!=s&&_x[1]!=s)_x[nroots++]=s;
00298     }
00299     return nroots;
00300   }
00301   /*If there are no odd powers, it's really just a quadratic in disguise.*/
00302   if(!_c)return rs_quadratic_solve(_gf,rs_gsqrt(_gf,_b),rs_gsqrt(_gf,_d),_x);
00303   /*Factor into (x**2 + r*x + s)*(x**2 + r*x + t) by solving for r, which can
00304      be shown to satisfy r**3 + _b*r + _c == 0.*/
00305   nroots=rs_cubic_solve(_gf,0,_b,_c,_x);
00306   if(nroots<1){
00307     /*The Reed-Solomon code is only valid if we can find 4 distinct roots in
00308        GF(2**8).
00309       If the cubic does not factor into 3 (possibly duplicate) roots, then we
00310        know that the quartic must have a non-trivial irreducible factor.*/
00311     return 0;
00312   }
00313   r=_x[0];
00314   /*Now solve for s and t.*/
00315   b=rs_gdiv(_gf,_c,r);
00316   nroots=rs_quadratic_solve(_gf,b,_d,_x);
00317   if(nroots<2)return 0;
00318   s=_x[0];
00319   t=_x[1];
00320   /*_c=r*(s^t) was non-zero, so s and t must be distinct.
00321     But if z is a root of z**2 ^ r*z ^ s, then so is (z^r), and s = z*(z^r).
00322     Hence if z is also a root of z**2 ^ r*z ^ t, then t = s, a contradiction.
00323     Thus all four roots are distinct, if they exist.*/
00324   nroots=rs_quadratic_solve(_gf,r,s,_x);
00325   return nroots+rs_quadratic_solve(_gf,r,t,_x+nroots);
00326 }
00327 
00328 /*Polynomial arithmetic with coefficients in GF(2**8).*/
00329 
00330 static void rs_poly_zero(unsigned char *_p,int _dp1){
00331   memset(_p,0,_dp1*sizeof(*_p));
00332 }
00333 
00334 static void rs_poly_copy(unsigned char *_p,const unsigned char *_q,int _dp1){
00335   memcpy(_p,_q,_dp1*sizeof(*_p));
00336 }
00337 
00338 /*Multiply the polynomial by the free variable, x (shift the coefficients).
00339   The number of coefficients, _dp1, must be non-zero.*/
00340 static void rs_poly_mul_x(unsigned char *_p,const unsigned char *_q,int _dp1){
00341   memmove(_p+1,_q,(_dp1-1)*sizeof(*_p));
00342   _p[0]=0;
00343 }
00344 
00345 /*Divide the polynomial by the free variable, x (shift the coefficients).
00346   The number of coefficients, _dp1, must be non-zero.*/
00347 static void rs_poly_div_x(unsigned char *_p,const unsigned char *_q,int _dp1){
00348   memmove(_p,_q+1,(_dp1-1)*sizeof(*_p));
00349   _p[_dp1-1]=0;
00350 }
00351 
00352 /*Compute the first (d+1) coefficients of the product of a degree e and a
00353    degree f polynomial.*/
00354 static void rs_poly_mult(const rs_gf256 *_gf,unsigned char *_p,int _dp1,
00355  const unsigned char *_q,int _ep1,const unsigned char *_r,int _fp1){
00356   int m;
00357   int i;
00358   rs_poly_zero(_p,_dp1);
00359   m=_ep1<_dp1?_ep1:_dp1;
00360   for(i=0;i<m;i++)if(_q[i]!=0){
00361     unsigned logqi;
00362     int      n;
00363     int      j;
00364     n=_dp1-i<_fp1?_dp1-i:_fp1;
00365     logqi=_gf->log[_q[i]];
00366     for(j=0;j<n;j++)_p[i+j]^=rs_hgmul(_gf,_r[j],logqi);
00367   }
00368 }
00369 
00370 /*Decoding.*/
00371 
00372 /*Computes the syndrome of a codeword.*/
00373 static void rs_calc_syndrome(const rs_gf256 *_gf,int _m0,
00374  unsigned char *_s,int _npar,const unsigned char *_data,int _ndata){
00375   int i;
00376   int j;
00377   for(j=0;j<_npar;j++){
00378     unsigned alphaj;
00379     unsigned sj;
00380     sj=0;
00381     alphaj=_gf->log[_gf->exp[j+_m0]];
00382     for(i=0;i<_ndata;i++)sj=_data[i]^rs_hgmul(_gf,sj,alphaj);
00383     _s[j]=sj;
00384   }
00385 }
00386 
00387 /*Berlekamp-Peterson and Berlekamp-Massey Algorithms for error-location,
00388    modified to handle known erasures, from \cite{CC81}, p. 205.
00389   This finds the coefficients of the error locator polynomial.
00390   The roots are then found by looking for the values of alpha**n where
00391    evaluating the polynomial yields zero.
00392   Error correction is done using the error-evaluator equation on p. 207.
00393   @BOOK{CC81,
00394     author="George C. Clark, Jr and J. Bibb Cain",
00395     title="Error-Correction Coding for Digitial Communications",
00396     series="Applications of Communications Theory",
00397     publisher="Springer",
00398     address="New York, NY",
00399     month=Jun,
00400     year=1981
00401   }*/
00402 
00403 /*Initialize lambda to the product of (1-x*alpha**e[i]) for erasure locations
00404    e[i].
00405   Note that the user passes in array indices counting from the beginning of the
00406    data, while our polynomial indexes starting from the end, so
00407    e[i]=(_ndata-1)-_erasures[i].*/
00408 static void rs_init_lambda(const rs_gf256 *_gf,unsigned char *_lambda,int _npar,
00409  const unsigned char *_erasures,int _nerasures,int _ndata){
00410   int i;
00411   int j;
00412   rs_poly_zero(_lambda,(_npar<4?4:_npar)+1);
00413   _lambda[0]=1;
00414   for(i=0;i<_nerasures;i++)for(j=i+1;j>0;j--){
00415     _lambda[j]^=rs_hgmul(_gf,_lambda[j-1],_ndata-1-_erasures[i]);
00416   }
00417 }
00418 
00419 /*From \cite{CC81}, p. 216.
00420   Returns the number of errors detected (degree of _lambda).*/
00421 static int rs_modified_berlekamp_massey(const rs_gf256 *_gf,
00422  unsigned char *_lambda,const unsigned char *_s,unsigned char *_omega,int _npar,
00423  const unsigned char *_erasures,int _nerasures,int _ndata){
00424   unsigned char tt[256];
00425   int           n;
00426   int           l;
00427   int           k;
00428   int           i;
00429   /*Initialize _lambda, the error locator-polynomial, with the location of
00430      known erasures.*/
00431   rs_init_lambda(_gf,_lambda,_npar,_erasures,_nerasures,_ndata);
00432   rs_poly_copy(tt,_lambda,_npar+1);
00433   l=_nerasures;
00434   k=0;
00435   for(n=_nerasures+1;n<=_npar;n++){
00436     unsigned d;
00437     rs_poly_mul_x(tt,tt,n-k+1);
00438     d=0;
00439     for(i=0;i<=l;i++)d^=rs_gmul(_gf,_lambda[i],_s[n-1-i]);
00440     if(d!=0){
00441       unsigned logd;
00442       logd=_gf->log[d];
00443       if(l<n-k){
00444         int t;
00445         for(i=0;i<=n-k;i++){
00446           unsigned tti;
00447           tti=tt[i];
00448           tt[i]=rs_hgmul(_gf,_lambda[i],255-logd);
00449           _lambda[i]=_lambda[i]^rs_hgmul(_gf,tti,logd);
00450         }
00451         t=n-k;
00452         k=n-l;
00453         l=t;
00454       }
00455       else for(i=0;i<=l;i++)_lambda[i]=_lambda[i]^rs_hgmul(_gf,tt[i],logd);
00456     }
00457   }
00458   rs_poly_mult(_gf,_omega,_npar,_lambda,l+1,_s,_npar);
00459   return l;
00460 }
00461 
00462 /*Finds all the roots of an error-locator polynomial _lambda by evaluating it
00463    at successive values of alpha, and returns the positions of the associated
00464    errors in _epos.
00465   Returns the number of valid roots identified.*/
00466 static int rs_find_roots(const rs_gf256 *_gf,unsigned char *_epos,
00467  const unsigned char *_lambda,int _nerrors,int _ndata){
00468   unsigned alpha;
00469   int      nroots;
00470   int      i;
00471   nroots=0;
00472   if(_nerrors<=4){
00473     /*Explicit solutions for higher degrees are possible.
00474       Chien uses large lookup tables to solve quintics, and Truong et al. give
00475        special algorithms for degree up through 11, though they use exhaustive
00476        search (with reduced complexity) for some portions.
00477       Quartics are good enough for reading CDs, and represent a reasonable code
00478        complexity trade-off without requiring any extra tables.
00479       Note that _lambda[0] is always 1.*/
00480     _nerrors=rs_quartic_solve(_gf,_lambda[1],_lambda[2],_lambda[3],_lambda[4],
00481      _epos);
00482     for(i=0;i<_nerrors;i++)if(_epos[i]){
00483       alpha=_gf->log[_epos[i]];
00484       if((int)alpha<_ndata)_epos[nroots++]=alpha;
00485     }
00486     return nroots;
00487   }
00488   else for(alpha=0;(int)alpha<_ndata;alpha++){
00489     unsigned alphai;
00490     unsigned sum;
00491     sum=0;
00492     alphai=0;
00493     for(i=0;i<=_nerrors;i++){
00494       sum^=rs_hgmul(_gf,_lambda[_nerrors-i],alphai);
00495       alphai=_gf->log[_gf->exp[alphai+alpha]];
00496     }
00497     if(!sum)_epos[nroots++]=alpha;
00498   }
00499   return nroots;
00500 }
00501 
00502 /*Corrects a codeword with _ndata<256 bytes, of which the last _npar are parity
00503    bytes.
00504   Known locations of errors can be passed in the _erasures array.
00505   Twice as many (up to _npar) errors with a known location can be corrected
00506    compared to errors with an unknown location.
00507   Returns the number of errors corrected if successful, or a negative number if
00508    the message could not be corrected because too many errors were detected.*/
00509 int rs_correct(const rs_gf256 *_gf,int _m0,unsigned char *_data,int _ndata,
00510  int _npar,const unsigned char *_erasures,int _nerasures){
00511   /*lambda must have storage for at least five entries to avoid special cases
00512      in the low-degree polynomial solver.*/
00513   unsigned char lambda[256];
00514   unsigned char omega[256];
00515   unsigned char epos[256];
00516   unsigned char s[256];
00517   int           i;
00518   /*If we already have too many erasures, we can't possibly succeed.*/
00519   if(_nerasures>_npar)return -1;
00520   /*Compute the syndrome values.*/
00521   rs_calc_syndrome(_gf,_m0,s,_npar,_data,_ndata);
00522   /*Check for a non-zero value.*/
00523   for(i=0;i<_npar;i++)if(s[i]){
00524     int nerrors;
00525     int j;
00526     /*Construct the error locator polynomial.*/
00527     nerrors=rs_modified_berlekamp_massey(_gf,lambda,s,omega,_npar,
00528      _erasures,_nerasures,_ndata);
00529     /*If we can't locate any errors, we can't force the syndrome values to
00530        zero, and must have a decoding error.
00531       Conversely, if we have too many errors, there's no reason to even attempt
00532        the root search.*/
00533     if(nerrors<=0||nerrors-_nerasures>_npar-_nerasures>>1)return -1;
00534     /*Compute the locations of the errors.
00535       If they are not all distinct, or some of them were outside the valid
00536        range for our block size, we have a decoding error.*/
00537     if(rs_find_roots(_gf,epos,lambda,nerrors,_ndata)<nerrors)return -1;
00538     /*Now compute the error magnitudes.*/
00539     for(i=0;i<nerrors;i++){
00540       unsigned a;
00541       unsigned b;
00542       unsigned alpha;
00543       unsigned alphan1;
00544       unsigned alphan2;
00545       unsigned alphanj;
00546       alpha=epos[i];
00547       /*Evaluate omega at alpha**-1.*/
00548       a=0;
00549       alphan1=255-alpha;
00550       alphanj=0;
00551       for(j=0;j<_npar;j++){
00552         a^=rs_hgmul(_gf,omega[j],alphanj);
00553         alphanj=_gf->log[_gf->exp[alphanj+alphan1]];
00554       }
00555       /*Evaluate the derivative of lambda at alpha**-1
00556         All the odd powers vanish.*/
00557       b=0;
00558       alphan2=_gf->log[_gf->exp[alphan1<<1]];
00559       alphanj=alphan1+_m0*alpha%255;
00560       for(j=1;j<=_npar;j+=2){
00561         b^=rs_hgmul(_gf,lambda[j],alphanj);
00562         alphanj=_gf->log[_gf->exp[alphanj+alphan2]];
00563       }
00564       /*Apply the correction.*/
00565       _data[_ndata-1-alpha]^=rs_gdiv(_gf,a,b);
00566     }
00567     return nerrors;
00568   }
00569   return 0;
00570 }
00571 
00572 /*Encoding.*/
00573 
00574 /*Create an _npar-coefficient generator polynomial for a Reed-Solomon code
00575    with _npar<256 parity bytes.*/
00576 void rs_compute_genpoly(const rs_gf256 *_gf,int _m0,
00577  unsigned char *_genpoly,int _npar){
00578   int i;
00579   if(_npar<=0)return;
00580   rs_poly_zero(_genpoly,_npar);
00581   _genpoly[0]=1;
00582   /*Multiply by (x+alpha^i) for i = 1 ... _ndata.*/
00583   for(i=0;i<_npar;i++){
00584     unsigned alphai;
00585     int      n;
00586     int      j;
00587     n=i+1<_npar-1?i+1:_npar-1;
00588     alphai=_gf->log[_gf->exp[_m0+i]];
00589     for(j=n;j>0;j--)_genpoly[j]=_genpoly[j-1]^rs_hgmul(_gf,_genpoly[j],alphai);
00590     _genpoly[0]=rs_hgmul(_gf,_genpoly[0],alphai);
00591   }
00592 }
00593 
00594 /*Adds _npar<=_ndata parity bytes to an _ndata-_npar byte message.
00595   _data must contain room for _ndata<256 bytes.*/
00596 void rs_encode(const rs_gf256 *_gf,unsigned char *_data,int _ndata,
00597  const unsigned char *_genpoly,int _npar){
00598   unsigned char *lfsr;
00599   unsigned       d;
00600   int            i;
00601   int            j;
00602   if(_npar<=0)return;
00603   lfsr=_data+_ndata-_npar;
00604   rs_poly_zero(lfsr,_npar);
00605   for(i=0;i<_ndata-_npar;i++){
00606     d=_data[i]^lfsr[0];
00607     if(d){
00608       unsigned logd;
00609       logd=_gf->log[d];
00610       for(j=0;j<_npar-1;j++){
00611         lfsr[j]=lfsr[j+1]^rs_hgmul(_gf,_genpoly[_npar-1-j],logd);
00612       }
00613       lfsr[_npar-1]=rs_hgmul(_gf,_genpoly[0],logd);
00614     }
00615     else rs_poly_div_x(lfsr,lfsr,_npar);
00616   }
00617 }
00618 
00619 #if defined(RS_TEST_ENC)
00620 #include <stdio.h>
00621 #include <stdlib.h>
00622 
00623 int main(void){
00624   rs_gf256 gf;
00625   int      k;
00626   rs_gf256_init(&gf,QR_PPOLY);
00627   srand(0);
00628   for(k=0;k<64*1024;k++){
00629     unsigned char genpoly[256];
00630     unsigned char data[256];
00631     unsigned char epos[256];
00632     int           ndata;
00633     int           npar;
00634     int           nerrors;
00635     int           i;
00636     ndata=rand()&0xFF;
00637     npar=ndata>0?rand()%ndata:0;
00638     for(i=0;i<ndata-npar;i++)data[i]=rand()&0xFF;
00639     rs_compute_genpoly(&gf,QR_M0,genpoly,npar);
00640     rs_encode(&gf,QR_M0,data,ndata,genpoly,npar);
00641     /*Write a clean version of the codeword.*/
00642     printf("%i %i",ndata,npar);
00643     for(i=0;i<ndata;i++)printf(" %i",data[i]);
00644     printf(" 0\n");
00645     /*Write the correct output to compare the decoder against.*/
00646     fprintf(stderr,"Success!\n",nerrors);
00647     for(i=0;i<ndata;i++)fprintf(stderr,"%i%s",data[i],i+1<ndata?" ":"\n");
00648     if(npar>0){
00649       /*Corrupt it.*/
00650       nerrors=rand()%(npar+1);
00651       if(nerrors>0){
00652         /*This test is not quite correct: there could be so many errors it
00653            comes within (npar>>1) errors of another valid codeword.
00654           I don't know a simple way to test for that without trying to decode
00655            the corrupt codeword, though, which is the very code we're trying to
00656            test.*/
00657         if(nerrors<=npar>>1){
00658           fprintf(stderr,"Success!\n",nerrors);
00659           for(i=0;i<ndata;i++){
00660             fprintf(stderr,"%i%s",data[i],i+1<ndata?" ":"\n");
00661           }
00662         }
00663         else fprintf(stderr,"Failure.\n");
00664         fprintf(stderr,"Success!\n",nerrors);
00665         for(i=0;i<ndata;i++)fprintf(stderr,"%i%s",data[i],i+1<ndata?" ":"\n");
00666         for(i=0;i<ndata;i++)epos[i]=i;
00667         for(i=0;i<nerrors;i++){
00668           unsigned char e;
00669           int           ei;
00670           ei=rand()%(ndata-i)+i;
00671           e=epos[ei];
00672           epos[ei]=epos[i];
00673           epos[i]=e;
00674           data[e]^=rand()%255+1;
00675         }
00676         /*First with no erasure locations.*/
00677         printf("%i %i",ndata,npar);
00678         for(i=0;i<ndata;i++)printf(" %i",data[i]);
00679         printf(" 0\n");
00680         /*Now with erasure locations.*/
00681         printf("%i %i",ndata,npar);
00682         for(i=0;i<ndata;i++)printf(" %i",data[i]);
00683         printf(" %i",nerrors);
00684         for(i=0;i<nerrors;i++)printf(" %i",epos[i]);
00685         printf("\n");
00686       }
00687     }
00688   }
00689   return 0;
00690 }
00691 #endif
00692 
00693 #if defined(RS_TEST_DEC)
00694 #include <stdio.h>
00695 
00696 int main(void){
00697   rs_gf256 gf;
00698   rs_gf256_init(&gf,QR_PPOLY);
00699   for(;;){
00700     unsigned char data[255];
00701     unsigned char erasures[255];
00702     int idata[255];
00703     int ierasures[255];
00704     int ndata;
00705     int npar;
00706     int nerasures;
00707     int nerrors;
00708     int i;
00709     if(scanf("%i",&ndata)<1||ndata<0||ndata>255||
00710      scanf("%i",&npar)<1||npar<0||npar>ndata)break;
00711     for(i=0;i<ndata;i++){
00712       if(scanf("%i",idata+i)<1||idata[i]<0||idata[i]>255)break;
00713       data[i]=idata[i];
00714     }
00715     if(i<ndata)break;
00716     if(scanf("%i",&nerasures)<1||nerasures<0||nerasures>ndata)break;
00717     for(i=0;i<nerasures;i++){
00718       if(scanf("%i",ierasures+i)<1||ierasures[i]<0||ierasures[i]>=ndata)break;
00719       erasures[i]=ierasures[i];
00720     }
00721     nerrors=rs_correct(&gf,QR_M0,data,ndata,npar,erasures,nerasures);
00722     if(nerrors>=0){
00723       unsigned char data2[255];
00724       unsigned char genpoly[255];
00725       for(i=0;i<ndata-npar;i++)data2[i]=data[i];
00726       rs_compute_genpoly(&gf,QR_M0,genpoly,npar);
00727       rs_encode(&gf,QR_M0,data2,ndata,genpoly,npar);
00728       for(i=ndata-npar;i<ndata;i++)if(data[i]!=data2[i]){
00729         printf("Abject failure! %i!=%i\n",data[i],data2[i]);
00730       }
00731       printf("Success!\n",nerrors);
00732       for(i=0;i<ndata;i++)printf("%i%s",data[i],i+1<ndata?" ":"\n");
00733     }
00734     else printf("Failure.\n");
00735   }
00736   return 0;
00737 }
00738 #endif
00739 
00740 #if defined(RS_TEST_ROOTS)
00741 #include <stdio.h>
00742 
00743 /*Exhaustively test the root finder.*/
00744 int main(void){
00745   rs_gf256 gf;
00746   int      a;
00747   int      b;
00748   int      c;
00749   int      d;
00750   rs_gf256_init(&gf,QR_PPOLY);
00751   for(a=0;a<256;a++)for(b=0;b<256;b++)for(c=0;c<256;c++)for(d=0;d<256;d++){
00752     unsigned char x[4];
00753     unsigned char r[4];
00754     unsigned      x2;
00755     unsigned      e[5];
00756     int           nroots;
00757     int           mroots;
00758     int           i;
00759     int           j;
00760     nroots=rs_quartic_solve(&gf,a,b,c,d,x);
00761     for(i=0;i<nroots;i++){
00762       x2=rs_gmul(&gf,x[i],x[i]);
00763       e[0]=rs_gmul(&gf,x2,x2)^rs_gmul(&gf,a,rs_gmul(&gf,x[i],x2))^
00764        rs_gmul(&gf,b,x2)^rs_gmul(&gf,c,x[i])^d;
00765       if(e[0]){
00766         printf("Invalid root: (0x%02X)**4 ^ 0x%02X*(0x%02X)**3 ^ "
00767          "0x%02X*(0x%02X)**2 ^ 0x%02X(0x%02X) ^ 0x%02X = 0x%02X\n",
00768          x[i],a,x[i],b,x[i],c,x[i],d,e[0]);
00769       }
00770       for(j=0;j<i;j++)if(x[i]==x[j]){
00771         printf("Repeated root %i=%i: (0x%02X)**4 ^ 0x%02X*(0x%02X)**3 ^ "
00772          "0x%02X*(0x%02X)**2 ^ 0x%02X(0x%02X) ^ 0x%02X = 0x%02X\n",
00773          i,j,x[i],a,x[i],b,x[i],c,x[i],d,e[0]);
00774       }
00775     }
00776     mroots=0;
00777     for(j=1;j<256;j++){
00778       int logx;
00779       int logx2;
00780       logx=gf.log[j];
00781       logx2=gf.log[gf.exp[logx<<1]];
00782       e[mroots]=d^rs_hgmul(&gf,c,logx)^rs_hgmul(&gf,b,logx2)^
00783        rs_hgmul(&gf,a,gf.log[gf.exp[logx+logx2]])^gf.exp[logx2<<1];
00784       if(!e[mroots])r[mroots++]=j;
00785     }
00786     /*We only care about missing roots if the quartic had 4 distinct, non-zero
00787        roots.*/
00788     if(mroots==4)for(j=0;j<mroots;j++){
00789       for(i=0;i<nroots;i++)if(x[i]==r[j])break;
00790       if(i>=nroots){
00791         printf("Missing root: (0x%02X)**4 ^ 0x%02X*(0x%02X)**3 ^ "
00792          "0x%02X*(0x%02X)**2 ^ 0x%02X(0x%02X) ^ 0x%02X = 0x%02X\n",
00793          r[j],a,r[j],b,r[j],c,r[j],d,e[j]);
00794       }
00795     }
00796   }
00797   return 0;
00798 }
00799 #endif
00800