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
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
Generated on Tue Jul 12 2022 18:54:12 by 1.7.2