summaryrefslogtreecommitdiffstats
path: root/zbar/qrcode/rs.c
diff options
context:
space:
mode:
Diffstat (limited to 'zbar/qrcode/rs.c')
-rw-r--r--zbar/qrcode/rs.c907
1 files changed, 907 insertions, 0 deletions
diff --git a/zbar/qrcode/rs.c b/zbar/qrcode/rs.c
new file mode 100644
index 0000000..bb195e9
--- /dev/null
+++ b/zbar/qrcode/rs.c
@@ -0,0 +1,907 @@
+/*Copyright (C) 1991-1995 Henry Minsky (hqm@ua.com, hqm@ai.mit.edu)
+ Copyright (C) 2008-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.*/
+#include "rs.h"
+#include <stdlib.h>
+#include <string.h>
+
+/*Reed-Solomon encoder and decoder.
+ Original implementation (C) Henry Minsky (hqm@ua.com, hqm@ai.mit.edu),
+ Universal Access 1991-1995.
+ Updates by Timothy B. Terriberry (C) 2008-2009:
+ - Properly reject codes when error-locator polynomial has repeated roots or
+ non-trivial irreducible factors.
+ - Removed the hard-coded parity size and performed general API cleanup.
+ - Allow multiple representations of GF(2**8), since different standards use
+ different irreducible polynomials.
+ - Allow different starting indices for the generator polynomial, since
+ different standards use different values.
+ - Greatly reduced the computation by eliminating unnecessary operations.
+ - Explicitly solve for the roots of low-degree polynomials instead of using
+ an exhaustive search.
+ This is another major speed boost when there are few errors.*/
+
+/*Galois Field arithmetic in GF(2**8).*/
+
+void rs_gf256_init(rs_gf256 *_gf, unsigned _ppoly)
+{
+ unsigned p;
+ int i;
+ /*Initialize the table of powers of a primtive root, alpha=0x02.*/
+ p = 1;
+ for (i = 0; i < 256; i++) {
+ _gf->exp[i] = _gf->exp[i + 255] = p;
+ p = ((p << 1) ^ (-(p >> 7) & _ppoly)) & 0xFF;
+ }
+ /*Invert the table to recover the logs.*/
+ for (i = 0; i < 255; i++)
+ _gf->log[_gf->exp[i]] = i;
+ /*Note that we rely on the fact that _gf->log[0]=0 below.*/
+ _gf->log[0] = 0;
+}
+
+/*Multiplication in GF(2**8) using logarithms.*/
+static unsigned rs_gmul(const rs_gf256 *_gf, unsigned _a, unsigned _b)
+{
+ return _a == 0 || _b == 0 ? 0 : _gf->exp[_gf->log[_a] + _gf->log[_b]];
+}
+
+/*Division in GF(2**8) using logarithms.
+ The result of division by zero is undefined.*/
+static unsigned rs_gdiv(const rs_gf256 *_gf, unsigned _a, unsigned _b)
+{
+ return _a == 0 ? 0 : _gf->exp[_gf->log[_a] + 255 - _gf->log[_b]];
+}
+
+/*Multiplication in GF(2**8) when one of the numbers is known to be non-zero
+ (proven by representing it by its logarithm).*/
+static unsigned rs_hgmul(const rs_gf256 *_gf, unsigned _a, unsigned _logb)
+{
+ return _a == 0 ? 0 : _gf->exp[_gf->log[_a] + _logb];
+}
+
+/*Square root in GF(2**8) using logarithms.*/
+static unsigned rs_gsqrt(const rs_gf256 *_gf, unsigned _a)
+{
+ unsigned loga;
+ if (!_a)
+ return 0;
+ loga = _gf->log[_a];
+ return _gf->exp[loga + (255 & -(loga & 1)) >> 1];
+}
+
+/*Polynomial root finding in GF(2**8).
+ Each routine returns a list of the distinct roots (i.e., with duplicates
+ removed).*/
+
+/*Solve a quadratic equation x**2 + _b*x + _c in GF(2**8) using the method
+ of~\cite{Wal99}.
+ Returns the number of distinct roots.
+ ARTICLE{Wal99,
+ author="C. Wayne Walker",
+ title="New Formulas for Solving Quadratic Equations over Certain Finite
+ Fields",
+ journal="{IEEE} Transactions on Information Theory",
+ volume=45,
+ number=1,
+ pages="283--284",
+ month=Jan,
+ year=1999
+ }*/
+static int rs_quadratic_solve(const rs_gf256 *_gf, unsigned _b, unsigned _c,
+ unsigned char _x[2])
+{
+ unsigned b;
+ unsigned logb;
+ unsigned logb2;
+ unsigned logb4;
+ unsigned logb8;
+ unsigned logb12;
+ unsigned logb14;
+ unsigned logc;
+ unsigned logc2;
+ unsigned logc4;
+ unsigned c8;
+ unsigned g3;
+ unsigned z3;
+ unsigned l3;
+ unsigned c0;
+ unsigned g2;
+ unsigned l2;
+ unsigned z2;
+ int inc;
+ /*If _b is zero, all we need is a square root.*/
+ if (!_b) {
+ _x[0] = rs_gsqrt(_gf, _c);
+ return 1;
+ }
+ /*If _c is zero, 0 and _b are the roots.*/
+ if (!_c) {
+ _x[0] = 0;
+ _x[1] = _b;
+ return 2;
+ }
+ logb = _gf->log[_b];
+ logc = _gf->log[_c];
+ /*If _b lies in GF(2**4), scale x to move it out.*/
+ inc = logb % (255 / 15) == 0;
+ if (inc) {
+ b = _gf->exp[logb + 254];
+ logb = _gf->log[b];
+ _c = _gf->exp[logc + 253];
+ logc = _gf->log[_c];
+ } else
+ b = _b;
+ logb2 = _gf->log[_gf->exp[logb << 1]];
+ logb4 = _gf->log[_gf->exp[logb2 << 1]];
+ logb8 = _gf->log[_gf->exp[logb4 << 1]];
+ logb12 = _gf->log[_gf->exp[logb4 + logb8]];
+ logb14 = _gf->log[_gf->exp[logb2 + logb12]];
+ logc2 = _gf->log[_gf->exp[logc << 1]];
+ logc4 = _gf->log[_gf->exp[logc2 << 1]];
+ c8 = _gf->exp[logc4 << 1];
+ g3 = rs_hgmul(_gf,
+ _gf->exp[logb14 + logc] ^ _gf->exp[logb12 + logc2] ^
+ _gf->exp[logb8 + logc4] ^ c8,
+ logb);
+ /*If g3 doesn't lie in GF(2**4), then our roots lie in an extension field.
+ Note that we rely on the fact that _gf->log[0]==0 here.*/
+ if (_gf->log[g3] % (255 / 15) != 0)
+ return 0;
+ /*Construct the corresponding quadratic in GF(2**4):
+ x**2 + x/alpha**(255/15) + l3/alpha**(2*(255/15))*/
+ z3 = rs_gdiv(_gf, g3, _gf->exp[logb8 << 1] ^ b);
+ l3 = rs_hgmul(_gf, rs_gmul(_gf, z3, z3) ^ rs_hgmul(_gf, z3, logb) ^ _c,
+ 255 - logb2);
+ c0 = rs_hgmul(_gf, l3, 255 - 2 * (255 / 15));
+ /*Construct the corresponding quadratic in GF(2**2):
+ x**2 + x/alpha**(255/3) + l2/alpha**(2*(255/3))*/
+ g2 =
+ rs_hgmul(_gf,
+ rs_hgmul(_gf, c0, 255 - 2 * (255 / 15)) ^ rs_gmul(_gf, c0, c0),
+ 255 - 255 / 15);
+ z2 = rs_gdiv(_gf, g2,
+ _gf->exp[255 - (255 / 15) * 4] ^ _gf->exp[255 - (255 / 15)]);
+ l2 = rs_hgmul(_gf,
+ rs_gmul(_gf, z2, z2) ^ rs_hgmul(_gf, z2, 255 - (255 / 15)) ^
+ c0,
+ 2 * (255 / 15));
+ /*Back substitute to the solution in the original field.*/
+ _x[0] = _gf->exp[_gf->log[z3 ^ rs_hgmul(_gf,
+ rs_hgmul(_gf, l2, 255 / 3) ^
+ rs_hgmul(_gf, z2, 255 / 15),
+ logb)] +
+ inc];
+ _x[1] = _x[0] ^ _b;
+ return 2;
+}
+
+/*Solve a cubic equation x**3 + _a*x**2 + _b*x + _c in GF(2**8).
+ Returns the number of distinct roots.*/
+static int rs_cubic_solve(const rs_gf256 *_gf, unsigned _a, unsigned _b,
+ unsigned _c, unsigned char _x[3])
+{
+ unsigned k;
+ unsigned logd;
+ unsigned d2;
+ unsigned logd2;
+ unsigned logw;
+ int nroots;
+ /*If _c is zero, factor out the 0 root.*/
+ if (!_c) {
+ nroots = rs_quadratic_solve(_gf, _a, _b, _x);
+ if (_b)
+ _x[nroots++] = 0;
+ return nroots;
+ }
+ /*Substitute x=_a+y*sqrt(_a**2+_b) to get y**3 + y + k == 0,
+ k = (_a*_b+c)/(_a**2+b)**(3/2).*/
+ k = rs_gmul(_gf, _a, _b) ^ _c;
+ d2 = rs_gmul(_gf, _a, _a) ^ _b;
+ if (!d2) {
+ int logx;
+ if (!k) {
+ /*We have a triple root.*/
+ _x[0] = _a;
+ return 1;
+ }
+ logx = _gf->log[k];
+ if (logx % 3 != 0)
+ return 0;
+ logx /= 3;
+ _x[0] = _a ^ _gf->exp[logx];
+ _x[1] = _a ^ _gf->exp[logx + 255 / 3];
+ _x[2] = _a ^ _x[0] ^ _x[1];
+ return 3;
+ }
+ logd2 = _gf->log[d2];
+ logd = logd2 + (255 & -(logd2 & 1)) >> 1;
+ k = rs_gdiv(_gf, k, _gf->exp[logd + logd2]);
+ /*Substitute y=w+1/w and z=w**3 to get z**2 + k*z + 1 == 0.*/
+ nroots = rs_quadratic_solve(_gf, k, 1, _x);
+ if (nroots < 1) {
+ /*The Reed-Solomon code is only valid if we can find 3 distinct roots in
+ GF(2**8), so if we know there's only one, we don't actually need to find
+ it.
+ Note that we're also called by the quartic solver, but if we contain a
+ non-trivial irreducible factor, than so does the original
+ quartic~\cite{LW72}, and failing to return a root here actually saves us
+ some work there, also.*/
+ return 0;
+ }
+ /*Recover w from z.*/
+ logw = _gf->log[_x[0]];
+ if (logw) {
+ if (logw % 3 != 0)
+ return 0;
+ logw /= 3;
+ /*Recover x from w.*/
+ _x[0] =
+ _gf->exp[_gf->log[_gf->exp[logw] ^ _gf->exp[255 - logw]] + logd] ^
+ _a;
+ logw += 255 / 3;
+ _x[1] =
+ _gf->exp[_gf->log[_gf->exp[logw] ^ _gf->exp[255 - logw]] + logd] ^
+ _a;
+ _x[2] = _x[0] ^ _x[1] ^ _a;
+ return 3;
+ } else {
+ _x[0] = _a;
+ /*In this case _x[1] is a double root, so we know the Reed-Solomon code is
+ invalid.
+ Note that we still have to return at least one root, because if we're
+ being called by the quartic solver, the quartic might still have 4
+ distinct roots.
+ But we don't need more than one root, so we can avoid computing the
+ expensive one.*/
+ /*_x[1]=_gf->exp[_gf->log[_gf->exp[255/3]^_gf->exp[2*(255/3)]]+logd]^_a;*/
+ return 1;
+ }
+}
+
+/*Solve a quartic equation x**4 + _a*x**3 + _b*x**2 + _c*x + _d in GF(2**8) by
+ decomposing it into the cases given by~\cite{LW72}.
+ Returns the number of distinct roots.
+ @ARTICLE{LW72,
+ author="Philip A. Leonard and Kenneth S. Williams",
+ title="Quartics over $GF(2^n)$",
+ journal="Proceedings of the American Mathematical Society",
+ volume=36,
+ number=2,
+ pages="347--450",
+ month=Dec,
+ year=1972
+ }*/
+static int rs_quartic_solve(const rs_gf256 *_gf, unsigned _a, unsigned _b,
+ unsigned _c, unsigned _d, unsigned char _x[3])
+{
+ unsigned r;
+ unsigned s;
+ unsigned t;
+ unsigned b;
+ int nroots;
+ int i;
+ /*If _d is zero, factor out the 0 root.*/
+ if (!_d) {
+ nroots = rs_cubic_solve(_gf, _a, _b, _c, _x);
+ if (_c)
+ _x[nroots++] = 0;
+ return nroots;
+ }
+ if (_a) {
+ unsigned loga;
+ /*Substitute x=(1/y) + sqrt(_c/_a) to eliminate the cubic term.*/
+ loga = _gf->log[_a];
+ r = rs_hgmul(_gf, _c, 255 - loga);
+ s = rs_gsqrt(_gf, r);
+ t = _d ^ rs_gmul(_gf, _b, r) ^ rs_gmul(_gf, r, r);
+ if (t) {
+ unsigned logti;
+ logti = 255 - _gf->log[t];
+ /*The result is still quartic, but with no cubic term.*/
+ nroots = rs_quartic_solve(
+ _gf, 0, rs_hgmul(_gf, _b ^ rs_hgmul(_gf, s, loga), logti),
+ _gf->exp[loga + logti], _gf->exp[logti], _x);
+ for (i = 0; i < nroots; i++)
+ _x[i] = _gf->exp[255 - _gf->log[_x[i]]] ^ s;
+ } else {
+ /*s must be a root~\cite{LW72}, and is in fact a double-root~\cite{CCO69}.
+ Thus we're left with only a quadratic to solve.
+ @ARTICLE{CCO69,
+ author="Robert T. Chien and B. D. Cunningham and I. B. Oldham",
+ title="Hybrid Methods for Finding Roots of a Polynomial---With
+ Applications to {BCH} Decoding",
+ journal="{IEEE} Transactions on Information Theory",
+ volume=15,
+ number=2,
+ pages="329--335",
+ month=Mar,
+ year=1969
+ }*/
+ nroots = rs_quadratic_solve(_gf, _a, _b ^ r, _x);
+ /*s may be a triple root if s=_b/_a, but not quadruple, since _a!=0.*/
+ if (nroots != 2 || _x[0] != s && _x[1] != s)
+ _x[nroots++] = s;
+ }
+ return nroots;
+ }
+ /*If there are no odd powers, it's really just a quadratic in disguise.*/
+ if (!_c)
+ return rs_quadratic_solve(_gf, rs_gsqrt(_gf, _b), rs_gsqrt(_gf, _d),
+ _x);
+ /*Factor into (x**2 + r*x + s)*(x**2 + r*x + t) by solving for r, which can
+ be shown to satisfy r**3 + _b*r + _c == 0.*/
+ nroots = rs_cubic_solve(_gf, 0, _b, _c, _x);
+ if (nroots < 1) {
+ /*The Reed-Solomon code is only valid if we can find 4 distinct roots in
+ GF(2**8).
+ If the cubic does not factor into 3 (possibly duplicate) roots, then we
+ know that the quartic must have a non-trivial irreducible factor.*/
+ return 0;
+ }
+ r = _x[0];
+ /*Now solve for s and t.*/
+ b = rs_gdiv(_gf, _c, r);
+ nroots = rs_quadratic_solve(_gf, b, _d, _x);
+ if (nroots < 2)
+ return 0;
+ s = _x[0];
+ t = _x[1];
+ /*_c=r*(s^t) was non-zero, so s and t must be distinct.
+ But if z is a root of z**2 ^ r*z ^ s, then so is (z^r), and s = z*(z^r).
+ Hence if z is also a root of z**2 ^ r*z ^ t, then t = s, a contradiction.
+ Thus all four roots are distinct, if they exist.*/
+ nroots = rs_quadratic_solve(_gf, r, s, _x);
+ return nroots + rs_quadratic_solve(_gf, r, t, _x + nroots);
+}
+
+/*Polynomial arithmetic with coefficients in GF(2**8).*/
+
+static void rs_poly_zero(unsigned char *_p, int _dp1)
+{
+ memset(_p, 0, _dp1 * sizeof(*_p));
+}
+
+static void rs_poly_copy(unsigned char *_p, const unsigned char *_q, int _dp1)
+{
+ memcpy(_p, _q, _dp1 * sizeof(*_p));
+}
+
+/*Multiply the polynomial by the free variable, x (shift the coefficients).
+ The number of coefficients, _dp1, must be non-zero.*/
+static void rs_poly_mul_x(unsigned char *_p, const unsigned char *_q, int _dp1)
+{
+ memmove(_p + 1, _q, (_dp1 - 1) * sizeof(*_p));
+ _p[0] = 0;
+}
+
+/*Divide the polynomial by the free variable, x (shift the coefficients).
+ The number of coefficients, _dp1, must be non-zero.*/
+static void rs_poly_div_x(unsigned char *_p, const unsigned char *_q, int _dp1)
+{
+ memmove(_p, _q + 1, (_dp1 - 1) * sizeof(*_p));
+ _p[_dp1 - 1] = 0;
+}
+
+/*Compute the first (d+1) coefficients of the product of a degree e and a
+ degree f polynomial.*/
+static void rs_poly_mult(const rs_gf256 *_gf, unsigned char *_p, int _dp1,
+ const unsigned char *_q, int _ep1,
+ const unsigned char *_r, int _fp1)
+{
+ int m;
+ int i;
+ rs_poly_zero(_p, _dp1);
+ m = _ep1 < _dp1 ? _ep1 : _dp1;
+ for (i = 0; i < m; i++)
+ if (_q[i] != 0) {
+ unsigned logqi;
+ int n;
+ int j;
+ n = _dp1 - i < _fp1 ? _dp1 - i : _fp1;
+ logqi = _gf->log[_q[i]];
+ for (j = 0; j < n; j++)
+ _p[i + j] ^= rs_hgmul(_gf, _r[j], logqi);
+ }
+}
+
+/*Decoding.*/
+
+/*Computes the syndrome of a codeword.*/
+static void rs_calc_syndrome(const rs_gf256 *_gf, int _m0, unsigned char *_s,
+ int _npar, const unsigned char *_data, int _ndata)
+{
+ int i;
+ int j;
+ for (j = 0; j < _npar; j++) {
+ unsigned alphaj;
+ unsigned sj;
+ sj = 0;
+ alphaj = _gf->log[_gf->exp[j + _m0]];
+ for (i = 0; i < _ndata; i++)
+ sj = _data[i] ^ rs_hgmul(_gf, sj, alphaj);
+ _s[j] = sj;
+ }
+}
+
+/*Berlekamp-Peterson and Berlekamp-Massey Algorithms for error-location,
+ modified to handle known erasures, from \cite{CC81}, p. 205.
+ This finds the coefficients of the error locator polynomial.
+ The roots are then found by looking for the values of alpha**n where
+ evaluating the polynomial yields zero.
+ Error correction is done using the error-evaluator equation on p. 207.
+ @BOOK{CC81,
+ author="George C. Clark, Jr and J. Bibb Cain",
+ title="Error-Correction Coding for Digital Communications",
+ series="Applications of Communications Theory",
+ publisher="Springer",
+ address="New York, NY",
+ month=Jun,
+ year=1981
+ }*/
+
+/*Initialize lambda to the product of (1-x*alpha**e[i]) for erasure locations
+ e[i].
+ Note that the user passes in array indices counting from the beginning of the
+ data, while our polynomial indexes starting from the end, so
+ e[i]=(_ndata-1)-_erasures[i].*/
+static void rs_init_lambda(const rs_gf256 *_gf, unsigned char *_lambda,
+ int _npar, const unsigned char *_erasures,
+ int _nerasures, int _ndata)
+{
+ int i;
+ int j;
+ rs_poly_zero(_lambda, (_npar < 4 ? 4 : _npar) + 1);
+ _lambda[0] = 1;
+ for (i = 0; i < _nerasures; i++)
+ for (j = i + 1; j > 0; j--) {
+ _lambda[j] ^=
+ rs_hgmul(_gf, _lambda[j - 1], _ndata - 1 - _erasures[i]);
+ }
+}
+
+/*From \cite{CC81}, p. 216.
+ Returns the number of errors detected (degree of _lambda).*/
+static int rs_modified_berlekamp_massey(const rs_gf256 *_gf,
+ unsigned char *_lambda,
+ const unsigned char *_s,
+ unsigned char *_omega, int _npar,
+ const unsigned char *_erasures,
+ int _nerasures, int _ndata)
+{
+ unsigned char tt[256];
+ int n;
+ int l;
+ int k;
+ int i;
+ /*Initialize _lambda, the error locator-polynomial, with the location of
+ known erasures.*/
+ rs_init_lambda(_gf, _lambda, _npar, _erasures, _nerasures, _ndata);
+ rs_poly_copy(tt, _lambda, _npar + 1);
+ l = _nerasures;
+ k = 0;
+ for (n = _nerasures + 1; n <= _npar; n++) {
+ unsigned d;
+ rs_poly_mul_x(tt, tt, n - k + 1);
+ d = 0;
+ for (i = 0; i <= l; i++)
+ d ^= rs_gmul(_gf, _lambda[i], _s[n - 1 - i]);
+ if (d != 0) {
+ unsigned logd;
+ logd = _gf->log[d];
+ if (l < n - k) {
+ int t;
+ for (i = 0; i <= n - k; i++) {
+ unsigned tti;
+ tti = tt[i];
+ tt[i] = rs_hgmul(_gf, _lambda[i], 255 - logd);
+ _lambda[i] = _lambda[i] ^ rs_hgmul(_gf, tti, logd);
+ }
+ t = n - k;
+ k = n - l;
+ l = t;
+ } else
+ for (i = 0; i <= l; i++)
+ _lambda[i] = _lambda[i] ^ rs_hgmul(_gf, tt[i], logd);
+ }
+ }
+ rs_poly_mult(_gf, _omega, _npar, _lambda, l + 1, _s, _npar);
+ return l;
+}
+
+/*Finds all the roots of an error-locator polynomial _lambda by evaluating it
+ at successive values of alpha, and returns the positions of the associated
+ errors in _epos.
+ Returns the number of valid roots identified.*/
+static int rs_find_roots(const rs_gf256 *_gf, unsigned char *_epos,
+ const unsigned char *_lambda, int _nerrors, int _ndata)
+{
+ unsigned alpha;
+ int nroots;
+ int i;
+ nroots = 0;
+ if (_nerrors <= 4) {
+ /*Explicit solutions for higher degrees are possible.
+ Chien uses large lookup tables to solve quintics, and Truong et al. give
+ special algorithms for degree up through 11, though they use exhaustive
+ search (with reduced complexity) for some portions.
+ Quartics are good enough for reading CDs, and represent a reasonable code
+ complexity trade-off without requiring any extra tables.
+ Note that _lambda[0] is always 1.*/
+ _nerrors = rs_quartic_solve(_gf, _lambda[1], _lambda[2], _lambda[3],
+ _lambda[4], _epos);
+ for (i = 0; i < _nerrors; i++)
+ if (_epos[i]) {
+ alpha = _gf->log[_epos[i]];
+ if ((int)alpha < _ndata)
+ _epos[nroots++] = alpha;
+ }
+ return nroots;
+ } else
+ for (alpha = 0; (int)alpha < _ndata; alpha++) {
+ unsigned alphai;
+ unsigned sum;
+ sum = 0;
+ alphai = 0;
+ for (i = 0; i <= _nerrors; i++) {
+ sum ^= rs_hgmul(_gf, _lambda[_nerrors - i], alphai);
+ alphai = _gf->log[_gf->exp[alphai + alpha]];
+ }
+ if (!sum)
+ _epos[nroots++] = alpha;
+ }
+ return nroots;
+}
+
+/*Corrects a codeword with _ndata<256 bytes, of which the last _npar are parity
+ bytes.
+ Known locations of errors can be passed in the _erasures array.
+ Twice as many (up to _npar) errors with a known location can be corrected
+ compared to errors with an unknown location.
+ Returns the number of errors corrected if successful, or a negative number if
+ the message could not be corrected because too many errors were detected.*/
+int rs_correct(const rs_gf256 *_gf, int _m0, unsigned char *_data, int _ndata,
+ int _npar, const unsigned char *_erasures, int _nerasures)
+{
+ /*lambda must have storage for at least five entries to avoid special cases
+ in the low-degree polynomial solver.*/
+ unsigned char lambda[256];
+ unsigned char omega[256];
+ unsigned char epos[256];
+ unsigned char s[256];
+ int i;
+ /*If we already have too many erasures, we can't possibly succeed.*/
+ if (_nerasures > _npar)
+ return -1;
+ /*Compute the syndrome values.*/
+ rs_calc_syndrome(_gf, _m0, s, _npar, _data, _ndata);
+ /*Check for a non-zero value.*/
+ for (i = 0; i < _npar; i++)
+ if (s[i]) {
+ int nerrors;
+ int j;
+ /*Construct the error locator polynomial.*/
+ nerrors = rs_modified_berlekamp_massey(_gf, lambda, s, omega, _npar,
+ _erasures, _nerasures,
+ _ndata);
+ /*If we can't locate any errors, we can't force the syndrome values to
+ zero, and must have a decoding error.
+ Conversely, if we have too many errors, there's no reason to even attempt
+ the root search.*/
+ if (nerrors <= 0 || nerrors - _nerasures > _npar - _nerasures >> 1)
+ return -1;
+ /*Compute the locations of the errors.
+ If they are not all distinct, or some of them were outside the valid
+ range for our block size, we have a decoding error.*/
+ if (rs_find_roots(_gf, epos, lambda, nerrors, _ndata) < nerrors)
+ return -1;
+ /*Now compute the error magnitudes.*/
+ for (i = 0; i < nerrors; i++) {
+ unsigned a;
+ unsigned b;
+ unsigned alpha;
+ unsigned alphan1;
+ unsigned alphan2;
+ unsigned alphanj;
+ alpha = epos[i];
+ /*Evaluate omega at alpha**-1.*/
+ a = 0;
+ alphan1 = 255 - alpha;
+ alphanj = 0;
+ for (j = 0; j < _npar; j++) {
+ a ^= rs_hgmul(_gf, omega[j], alphanj);
+ alphanj = _gf->log[_gf->exp[alphanj + alphan1]];
+ }
+ /*Evaluate the derivative of lambda at alpha**-1
+ All the odd powers vanish.*/
+ b = 0;
+ alphan2 = _gf->log[_gf->exp[alphan1 << 1]];
+ alphanj = alphan1 + _m0 * alpha % 255;
+ for (j = 1; j <= _npar; j += 2) {
+ b ^= rs_hgmul(_gf, lambda[j], alphanj);
+ alphanj = _gf->log[_gf->exp[alphanj + alphan2]];
+ }
+ /*Apply the correction.*/
+ _data[_ndata - 1 - alpha] ^= rs_gdiv(_gf, a, b);
+ }
+ return nerrors;
+ }
+ return 0;
+}
+
+/*Encoding.*/
+
+/*Create an _npar-coefficient generator polynomial for a Reed-Solomon code
+ with _npar<256 parity bytes.*/
+void rs_compute_genpoly(const rs_gf256 *_gf, int _m0, unsigned char *_genpoly,
+ int _npar)
+{
+ int i;
+ if (_npar <= 0)
+ return;
+ rs_poly_zero(_genpoly, _npar);
+ _genpoly[0] = 1;
+ /*Multiply by (x+alpha^i) for i = 1 ... _ndata.*/
+ for (i = 0; i < _npar; i++) {
+ unsigned alphai;
+ int n;
+ int j;
+ n = i + 1 < _npar - 1 ? i + 1 : _npar - 1;
+ alphai = _gf->log[_gf->exp[_m0 + i]];
+ for (j = n; j > 0; j--)
+ _genpoly[j] = _genpoly[j - 1] ^ rs_hgmul(_gf, _genpoly[j], alphai);
+ _genpoly[0] = rs_hgmul(_gf, _genpoly[0], alphai);
+ }
+}
+
+/*Adds _npar<=_ndata parity bytes to an _ndata-_npar byte message.
+ _data must contain room for _ndata<256 bytes.*/
+void rs_encode(const rs_gf256 *_gf, unsigned char *_data, int _ndata,
+ const unsigned char *_genpoly, int _npar)
+{
+ unsigned char *lfsr;
+ unsigned d;
+ int i;
+ int j;
+ if (_npar <= 0)
+ return;
+ lfsr = _data + _ndata - _npar;
+ rs_poly_zero(lfsr, _npar);
+ for (i = 0; i < _ndata - _npar; i++) {
+ d = _data[i] ^ lfsr[0];
+ if (d) {
+ unsigned logd;
+ logd = _gf->log[d];
+ for (j = 0; j < _npar - 1; j++) {
+ lfsr[j] = lfsr[j + 1] ^
+ rs_hgmul(_gf, _genpoly[_npar - 1 - j], logd);
+ }
+ lfsr[_npar - 1] = rs_hgmul(_gf, _genpoly[0], logd);
+ } else
+ rs_poly_div_x(lfsr, lfsr, _npar);
+ }
+}
+
+#if defined(RS_TEST_ENC)
+#include <stdio.h>
+#include <stdlib.h>
+
+int main(void)
+{
+ rs_gf256 gf;
+ int k;
+ rs_gf256_init(&gf, QR_PPOLY);
+ srand(0);
+ for (k = 0; k < 64 * 1024; k++) {
+ unsigned char genpoly[256];
+ unsigned char data[256];
+ unsigned char epos[256];
+ int ndata;
+ int npar;
+ int nerrors;
+ int i;
+ ndata = rand() & 0xFF;
+ npar = ndata > 0 ? rand() % ndata : 0;
+ for (i = 0; i < ndata - npar; i++)
+ data[i] = rand() & 0xFF;
+ rs_compute_genpoly(&gf, QR_M0, genpoly, npar);
+ rs_encode(&gf, QR_M0, data, ndata, genpoly, npar);
+ /*Write a clean version of the codeword.*/
+ printf("%i %i", ndata, npar);
+ for (i = 0; i < ndata; i++)
+ printf(" %i", data[i]);
+ printf(" 0\n");
+ /*Write the correct output to compare the decoder against.*/
+ fprintf(stderr, "Success!\n", nerrors);
+ for (i = 0; i < ndata; i++)
+ fprintf(stderr, "%i%s", data[i], i + 1 < ndata ? " " : "\n");
+ if (npar > 0) {
+ /*Corrupt it.*/
+ nerrors = rand() % (npar + 1);
+ if (nerrors > 0) {
+ /*This test is not quite correct: there could be so many errors it
+ comes within (npar>>1) errors of another valid codeword.
+ I don't know a simple way to test for that without trying to decode
+ the corrupt codeword, though, which is the very code we're trying to
+ test.*/
+ if (nerrors <= npar >> 1) {
+ fprintf(stderr, "Success!\n", nerrors);
+ for (i = 0; i < ndata; i++) {
+ fprintf(stderr, "%i%s", data[i],
+ i + 1 < ndata ? " " : "\n");
+ }
+ } else
+ fprintf(stderr, "Failure.\n");
+ fprintf(stderr, "Success!\n", nerrors);
+ for (i = 0; i < ndata; i++)
+ fprintf(stderr, "%i%s", data[i],
+ i + 1 < ndata ? " " : "\n");
+ for (i = 0; i < ndata; i++)
+ epos[i] = i;
+ for (i = 0; i < nerrors; i++) {
+ unsigned char e;
+ int ei;
+ ei = rand() % (ndata - i) + i;
+ e = epos[ei];
+ epos[ei] = epos[i];
+ epos[i] = e;
+ data[e] ^= rand() % 255 + 1;
+ }
+ /*First with no erasure locations.*/
+ printf("%i %i", ndata, npar);
+ for (i = 0; i < ndata; i++)
+ printf(" %i", data[i]);
+ printf(" 0\n");
+ /*Now with erasure locations.*/
+ printf("%i %i", ndata, npar);
+ for (i = 0; i < ndata; i++)
+ printf(" %i", data[i]);
+ printf(" %i", nerrors);
+ for (i = 0; i < nerrors; i++)
+ printf(" %i", epos[i]);
+ printf("\n");
+ }
+ }
+ }
+ return 0;
+}
+#endif
+
+#if defined(RS_TEST_DEC)
+#include <stdio.h>
+
+int main(void)
+{
+ rs_gf256 gf;
+ rs_gf256_init(&gf, QR_PPOLY);
+ for (;;) {
+ unsigned char data[255];
+ unsigned char erasures[255];
+ int idata[255];
+ int ierasures[255];
+ int ndata;
+ int npar;
+ int nerasures;
+ int nerrors;
+ int i;
+ if (scanf("%i", &ndata) < 1 || ndata < 0 || ndata > 255 ||
+ scanf("%i", &npar) < 1 || npar < 0 || npar > ndata)
+ break;
+ for (i = 0; i < ndata; i++) {
+ if (scanf("%i", idata + i) < 1 || idata[i] < 0 || idata[i] > 255)
+ break;
+ data[i] = idata[i];
+ }
+ if (i < ndata)
+ break;
+ if (scanf("%i", &nerasures) < 1 || nerasures < 0 || nerasures > ndata)
+ break;
+ for (i = 0; i < nerasures; i++) {
+ if (scanf("%i", ierasures + i) < 1 || ierasures[i] < 0 ||
+ ierasures[i] >= ndata)
+ break;
+ erasures[i] = ierasures[i];
+ }
+ nerrors =
+ rs_correct(&gf, QR_M0, data, ndata, npar, erasures, nerasures);
+ if (nerrors >= 0) {
+ unsigned char data2[255];
+ unsigned char genpoly[255];
+ for (i = 0; i < ndata - npar; i++)
+ data2[i] = data[i];
+ rs_compute_genpoly(&gf, QR_M0, genpoly, npar);
+ rs_encode(&gf, QR_M0, data2, ndata, genpoly, npar);
+ for (i = ndata - npar; i < ndata; i++)
+ if (data[i] != data2[i]) {
+ printf("Abject failure! %i!=%i\n", data[i], data2[i]);
+ }
+ printf("Success!\n", nerrors);
+ for (i = 0; i < ndata; i++)
+ printf("%i%s", data[i], i + 1 < ndata ? " " : "\n");
+ } else
+ printf("Failure.\n");
+ }
+ return 0;
+}
+#endif
+
+#if defined(RS_TEST_ROOTS)
+#include <stdio.h>
+
+/*Exhaustively test the root finder.*/
+int main(void)
+{
+ rs_gf256 gf;
+ int a;
+ int b;
+ int c;
+ int d;
+ rs_gf256_init(&gf, QR_PPOLY);
+ for (a = 0; a < 256; a++)
+ for (b = 0; b < 256; b++)
+ for (c = 0; c < 256; c++)
+ for (d = 0; d < 256; d++) {
+ unsigned char x[4];
+ unsigned char r[4];
+ unsigned x2;
+ unsigned e[5];
+ int nroots;
+ int mroots;
+ int i;
+ int j;
+ nroots = rs_quartic_solve(&gf, a, b, c, d, x);
+ for (i = 0; i < nroots; i++) {
+ x2 = rs_gmul(&gf, x[i], x[i]);
+ e[0] = rs_gmul(&gf, x2, x2) ^
+ rs_gmul(&gf, a, rs_gmul(&gf, x[i], x2)) ^
+ rs_gmul(&gf, b, x2) ^ rs_gmul(&gf, c, x[i]) ^ d;
+ if (e[0]) {
+ printf(
+ "Invalid root: (0x%02X)**4 ^ 0x%02X*(0x%02X)**3 ^ "
+ "0x%02X*(0x%02X)**2 ^ 0x%02X(0x%02X) ^ 0x%02X = 0x%02X\n",
+ x[i], a, x[i], b, x[i], c, x[i], d, e[0]);
+ }
+ for (j = 0; j < i; j++)
+ if (x[i] == x[j]) {
+ printf(
+ "Repeated root %i=%i: (0x%02X)**4 ^ 0x%02X*(0x%02X)**3 ^ "
+ "0x%02X*(0x%02X)**2 ^ 0x%02X(0x%02X) ^ 0x%02X = 0x%02X\n",
+ i, j, x[i], a, x[i], b, x[i], c, x[i], d,
+ e[0]);
+ }
+ }
+ mroots = 0;
+ for (j = 1; j < 256; j++) {
+ int logx;
+ int logx2;
+ logx = gf.log[j];
+ logx2 = gf.log[gf.exp[logx << 1]];
+ e[mroots] =
+ d ^ rs_hgmul(&gf, c, logx) ^
+ rs_hgmul(&gf, b, logx2) ^
+ rs_hgmul(&gf, a, gf.log[gf.exp[logx + logx2]]) ^
+ gf.exp[logx2 << 1];
+ if (!e[mroots])
+ r[mroots++] = j;
+ }
+ /*We only care about missing roots if the quartic had 4 distinct, non-zero
+ roots.*/
+ if (mroots == 4)
+ for (j = 0; j < mroots; j++) {
+ for (i = 0; i < nroots; i++)
+ if (x[i] == r[j])
+ break;
+ if (i >= nroots) {
+ printf(
+ "Missing root: (0x%02X)**4 ^ 0x%02X*(0x%02X)**3 ^ "
+ "0x%02X*(0x%02X)**2 ^ 0x%02X(0x%02X) ^ 0x%02X = 0x%02X\n",
+ r[j], a, r[j], b, r[j], c, r[j], d, e[j]);
+ }
+ }
+ }
+ return 0;
+}
+#endif
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-07 06:43:54 +0000