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