#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