// bch.cpp
//
// C++ code by Kevin Harmon and Leonid Reyzin (reyzin@cs.bu.edu)
// for sublinear-time syndrome encoding and decoding of
// binary BCH codes. See pinsketch.txt for information
// on using this inside PinSketch, BCH-based secure sketches.
//
// Uses Victor Shoup's NTL (see http://www.shoup.net)
//
// Contains two public functions: BCHSyndromeCompute
// and BCHSyndromeDecode
//
// See Syndrome Encoding and Decoding of BCH Codes in Sublinear Time
// (Excerpted from "Fuzzy Extractors:
// How to Generate Strong Keys from Biometrics and Other Noisy Data,"
// SIAM Journal on Computing, 38(1):87-139, 2008,
// http://arxiv.org/abs/cs/0602007)
// by Yevgeniy Dodis, Rafail Ostrovsky, Leonid Reyzin and Adam Smith
// (file bch-excerpt.pdf) for the mathematics behind this.
//
// This code and explanatory notes are hosted at
// http://www.cs.bu.edu/~reyzin/code/fuzzy.html
//
//
#include "pinsketch.h"
///////////////////////////////////////////////////////////////////////////
// PURPOSE:
// Computes the syndrome of a sparse vector
// of the binary BCH code of design distance d.
// The vector is viewed as a vector of 0's and 1's
// being indexed by all nonzero elements of GF2E; because
// it is sparse, it is given as the set a of
// elements of GF2E where the coordinates of the vector are equal to 1.
// If used to compute the secure sketch, the sketch will
// tolerate symmetric difference of up to (d-1)/2
//
//
// ALGORITHM:
// The syndrome is computed as a vector of
// f(j) = (a_0)^j + (a_2)^j + ... + (a_s)^j
// for odd i from 1 do d-1, where a_i is the i-th component
// of the input vector A.
// (only the odd j are needed, because
// f(2j) is simply the square of f(j)).
// Because in C++ we number from 0, f(j) will reside
// in location (j-1)/2.
//
//
// ASSUMPTIONS:
// Let m=GF2E::degree() (i.e., the field is GF(2^m)).
// Assumes d is odd,
// greater than 1, and less than 2^m (else BCH codes don't make sense).
// Assumes the input set has no zeros (they will be ignored)
//
//
// RUNNING TIME:
// Takes time O(len*d) operations in GF(2^m),
// where len is the length of the input vector
//
void BCHSyndromeCompute(vec_GF2E & ss, const vec_GF2E & a, long d)
{
GF2E a_i_to_the_j, multiplier;
long i, j;
ss.SetLength((d-1)/2); // half the syndrome length,
// because even power not needed
// We will compute the fs in parallel: first add
// all the powers of a_1, then of a_2, ..., then of a_s
for (i = 0; i < a.length(); ++i)
{
a_i_to_the_j = a[i];
sqr(multiplier, a[i]); // multiplier = a[i]*a[i];
// special-case 0, because it doesn't need to be multiplied
// by the multiplier
ss[0] += a_i_to_the_j;
for (long j = 3; j < d; j+=2)
{
a_i_to_the_j *= multiplier;
ss[(j-1)/2] += a_i_to_the_j;
}
}
}
///////////////////////////////////////////////////////////////////////////
// Produces a vector res such that res[2i]=ss[i]
// and res[2i+1]=ss[i]*ss[i]
//
// Used to recover the redundant representation
// of the BCH syndrome (which includes even values of j)
// from the representation produced by BCHSyndromeCompute
// Because C++ indexes from 0, the j-th coordinate of the syndrome
// will end up in location j-1.
//
// Takes time O(d) operations in GF2E, where d is the output length
//
static
void InterpolateEvens(vec_GF2E & res, const vec_GF2E & ss)
{
// uses relation syn(j) = syn(j/2)^2 to recover syn from ss
long i;
res.SetLength(2*ss.length());
// odd coordinates (which, confusingly, means even i)
// are just copied from the input
for (i = 0; i < ss.length(); ++i)
res[2*i] = ss[i];
// even coordinates (odd i) are computed via squaring.
for (i = 1; i < res.length(); i+=2)
sqr(res[i], res[(i-1)/2]); // square
}
///////////////////////////////////////////////////////////////////////////
// PURPOSE:
// Returns true if f fully factors into distinct roots
// (i.e., if f is a product of distinct monic degree-1 polynomials
// times possibly a constant)
// and false otherwise.
// If f is zero, returns false.
//
// ALGORITHM:
// Let m=GF2E::degree() (i.e., the field is GF(2^m)).
// The check is accomplished by checking if f divides X^{2^m} - X,
// or equivalently if X^{2^m}-X is 0 mod f.
// X^{2^m} - X has 2^m distinct roots -- namely,
// every element of the field is a root. Hence, f divides it if and only
// if f has all its roots and they are all distinct.
//
// RUNNING TIME:
// Depends on NTL's implementation of FrobeniusMap, but for inputs of degree
// e that is relatively small compared m, should take e^{\log_2 3} m
// operations in GF(2^m). Note that \log_2 3 is about 1.585.
static
bool CheckIfDistinctRoots(const GF2EX & f)
{
if (IsZero(f))
return false;
// We hanlde degree 0 and degree 1 case separately, so that later
// we can assume X mod f is the same as X
if (deg(f) == 0 || deg(f) == 1)
return true;
GF2EXModulus F;
// F is the same as f, just more efficient modular operations
build(F, f);
GF2EX h;
FrobeniusMap (h, F); // h = X^{2^m} mod F
// If X^{2^m} - X = 0 mod F, then X^{2^m} mod F
// should be just X (because degree of F > 1)
return (IsX(h));
}
///////////////////////////////////////////////////////////////////////////
// PURPOSE:
// Given syndrome ssWithoutEvens of BCH code with design distance d,
// finds sparse vector (with no more than
// (d-1)/2 ones) with that syndrome
// (note that syndrome as well sparse vector
// representation are defined at BCHSyndromeCompute above).
// 'answer' returns positions of ones in the resulting vector.
// These positions are elements of GF2E (i.e., we view the vector
// as a vector whose positions are indexed by elements of GF2E).
// Returns false if no such vector exists, true otherwise
// The input syndrome is assumed to not have even powers, i.e.,
// has (d-1)/2 elements, such as the syndrome computed by BCHSyndromeCompute.
//
// ASSUMPTIONS:
// Let m=GF2E::degree() (i.e., the field is GF(2^m)).
// This algorithm assumes that d is odd, greater than 1, and less than 2^m.
// (else BCH codes don't make sense).
// Assumes input is of length (d-1)/2.
//
// ALGORITHM USED:
// Implements BCH decoding based on Euclidean algorithm;
// For the explanation of the algorithm, see
// Syndrome Encoding and Decoding of BCH Codes in Sublinear Time
// (Excerpted from Fuzzy Extractors:
// How to Generate Strong Keys from Biometrics and Other Noisy Data)
// by Yevgeniy Dodis, Rafail Ostrovsky, Leonid Reyzin and Adam Smith
// or Theorem 18.7 of Victor Shoup's "A Computational Introduction to
// Number Theory and Algebra" (first edition, 2005), or
// pp. 170-173 of "Introduction to Coding Theory" by Jurgen Bierbrauer.
//
//
// RUNNING TIME:
// If the output has e elements (i.e., the length of the output vector
// is e; note that e <= (d-1)/2), then
// the running time is O(d^2 + e^2 + e^{\log_2 3} m) operations in GF(2^m),
// each of which takes time O(m^{\log_2 3}) in NTL. Note that
// \log_2 3 is approximately 1.585.
//
bool BCHSyndromeDecode(vec_GF2E &answer, const vec_GF2E & ssWithoutEvens, long d)
{
long i;
vec_GF2E ss;
// This takes O(d) operation in GF(2^m)
InterpolateEvens(ss, ssWithoutEvens);
GF2EX r1, r2, r3, v1, v2, v3, q, temp;
GF2EX *Rold, *Rcur, *Rnew, *Vold, *Vcur, *Vnew, *tempPointer;
// Use pointers to avoid moving polynomials around
// An assignment of polynomials requires copying the coefficient vector;
// we will not assign polynomials, but will swap pointers instead
Rold = &r1;
Rcur = &r2;
Rnew = &r3;
Vold = &v1;
Vcur = &v2;
Vnew = &v3;
SetCoeff(*Rold, d-1, 1); // Rold holds z^{d-1}
// Rcur=S(z)/z where S is the syndrome poly, Rcur = \sum S_j z^{j-1}
// Note that because we index arrays from 0, S_j is stored in ss[j-1]
for (i=0; i= t) {
// Rold = Rcur*q + Rnew
DivRem(q, *Rnew, *Rold, *Rcur);
// Vnew = Vold - qVcur)
mul(temp, q, *Vcur);
sub (*Vnew, *Vold, temp);
// swap everything
tempPointer = Rold;
Rold = Rcur;
Rcur = Rnew;
Rnew = tempPointer;
tempPointer = Vold;
Vold = Vcur;
Vcur = Vnew;
Vnew = tempPointer;
}
// At the end of the loop, sigma(z) is Vcur
// (up to a constant factor, which doesn't matter,
// since we care about roots of sigma).
// The roots of sigma(z) are inverses of the points we
// are interested in.
// We will check that 0 is not
// a root of Vcur (else its inverse won't exist, and hence
// the right polynomial doesn't exist).
if (IsZero(ConstTerm(*Vcur)))
return false;
// Need sigma to be monic for FindRoots
MakeMonic(*Vcur);
// check if sigma(z) has distinct roots if not, return false
// this will take O(e^{\log_2 3} m) operations in GF(2^m),
// where e is the degree of sigma(z)
if (CheckIfDistinctRoots(*Vcur) == false)
return false;
// find roots of sigma(z)
// this will take O(e^2 + e^{\log_2 3} m) operations in GF(2^m),
// where e is the degree of sigma(z)
answer = FindRoots(*Vcur);
// take inverses of roots of sigma(z)
for (i = 0; i < answer.length(); ++i)
answer[i] = inv(answer[i]);
// It is now necessary to verify if the resulting vector
// has the correct syndrome: it is possible that it does
// not even though the polynomial sigma(z) factors
// completely
// This takes O(de) operations in GF(2^m)
vec_GF2E test;
BCHSyndromeCompute (test, answer, d);
return (test==ssWithoutEvens);
}