#include "NTL/mat_ZZ.h"
#include "NTL/ZZ_pX.h"
NTL_CLIENT
static bool LinearSolve (vec_ZZ_p & solution, mat_ZZ_p & A);
/*
* WelchBerlekamp.cpp
* Authors: Soren Johnson and Leonid Reyzin (reyzin@cs.bu.edu)
*
* This code and explanatory notes are hosted at
* http://www.cs.bu.edu/~reyzin/code/fuzzy.html
*
* Uses Victor Shoup's NTL, see http://www.shoup.net
*
* This function decodes [n,n-2t,2t+1] Reed-Solomon codes over GF(p) using
* the Welch-Berlekamp algorithm, as described in Section 3 of Notes for
* Lecture 10 of Madhu Sudan's 2001 course "Algorithmic Introduction to
* Coding Theory" http://people.csail.mit.edu/madhu/FT01/ (that description
* is based, in turn, on a description in "Highly Resilient Correctors for
* Multivariate Polynomials," by Gemmel and Sudan, Information Processing
* Letters, 43(4):169-174, 1992).
*
* Namely, given vectors X and Y of length n each, it finds a polynomial f
* over GF(p) of degree at most n-2t-1 such that f(X[i])=Y[i] for at least
* n-t values of i. Such a polynomial f, if it exists, is unique by the
* funamental theorem of algebra (if there were two such polynomials f_1,
* f_2, then f_1-f_2 would have at least n-2t roots, but degree at most
* n-2t-1).
*
* Assumes (does not check!) that no value appears in X more than once,
* that X and Y are of the same length, and that ZZ_p::init has been run
* and that ZZ_p::modulus() is a prime.
*
* The polynomial f is placed into the "result" argument. Returns true if
* successful, false if such a polynomial does not exist.
*
* This algorithm takes Theta(n^3) operations over GF(p). It is not the
* best known algorithm for this problem; faster algorithms are referenced
* in the lecture notes mentioned above.
*
* NOTE for those using this code as part of Improved Juels-Sudan Secure
* Sketch: the meaning of n and t in this function is different from the
* meaning of n and t in ijs.cpp and ijsio.cpp, in order to be consistent
* with standard coding notation.
*/
bool WelchBerlekamp (ZZ_pX & result, const vec_ZZ_p & X, const vec_ZZ_p & Y, int t)
{
int n = X.length();
if (t<0)
return false;
if (t==0) { // special case t=0: then no errors are allowed,
interpolate(result, X, Y); // and result is just an interpolation of X, Y.
return true; // We special-case it here to avoid
} // creating a special case later in the code
if (n<2*t) // degree n-2t-1<-1 doesn't make sense if n<2t.
return false;
// Note that the case of n=2t is possible: the result should be the 0
// polynomial if at least half of the Y values are 0, and false otherwise
// The decoding algorithm handles that correctly below
int i, j;
mat_ZZ_p A; // the matrix to hold the linear system we'll be solving
vec_ZZ_p NE; // coefficients of N and E as one vector
A.SetDims(n, n + 1);
/* We know that N has degree at most n-t-1,
E has degree exactly t and is monic,
and N(X[i])=Y[i]E(X[i]) for all i. Thus, we need to solve
the following matrix equation to find the coefficients of N and E:
[1 x0 x0^2 ... x0^(n-t-1) -y0 -y0x0 -y0x0^2 ... -y0x0^(t-1)] [N0 ] [y0x0^t]
[1 x1 x1^2 ... x1^(n-t-1) -y1 -y1x1 -y1x1^2 ... -y1x1^(t-1)] [N1 ] [y1x1^t]
... ... ...
... [N(n-t-1)] ...
... [E0 ] = ...
... [E1 ] ...
... ... ...
[1 xm xm^2 ... xm^(n-t-1) -ym -ymxm -ymxm^2 ... -ymxm^(t-1)] [E(t-1) ] [ymxm^t]
where m = n - 1 (it just got too messy to write (n-1) at the bottom).
*/
int numNCoeffs=n-t;
// It is important that numNCoeffs=0; row--) {
// remove all the known variables from the equation
t = A[row][n];
int col;
for (int col = n-1; col>=firstDeterminedValue; col--)
t-=(A[row][col]*solution[col]);
// now we need to find the first nonzero coefficient in this row
// if it exists before firstDeterminedValue
// because the matrix is in row echelon form, the first nonzero
// coefficient cannot be before the diagonal
for (col=row; col