rsdecoder.h

// Percy++ Copyright 2007,2012,2013,2014
// Ian Goldberg <iang@cs.uwaterloo.ca>,
// Casey Devet <cjdevet@uwaterloo.ca>,
// Ann Yang <y242yang@uwaterloo.ca>
// 
// This program is free software; you can redistribute it and/or modify
// it under the terms of version 2 of the GNU General Public License as
// published by the Free Software Foundation.
// 
// This program 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 General Public License for more details.
// 
// There is a copy of the GNU General Public License in the COPYING file
// packaged with this plugin; if you cannot find it, write to the Free
// Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// 02110-1301 USA

#ifndef __GSDECODER_H__
#define __GSDECODER_H__

#include <vector>
#include <map>
#include <set>
#include <string>
#include <NTL/mat_ZZ_p.h>
#include <NTL/mat_GF2E.h>
#include "percyresult.h"
#include "FXY.h"
#include "portfolio.h"

// An implementation of the Guruswami-Sudan decoder.  This decoder will
// produce all possible codewords of degree at most t that match the
// given codeword with k non-erasures in more than sqrt(k*t) places.
//
// It runs in polynomial time, but in the worst case, it's O(k^12).
// Luckily, the worst case is when t=k-2, which is easily solved by
// brute force.  (You will need to do this brute force outside this
// class.)
//
// This implementation is based on the K\"otter and Roth-Ruckenstein
// algorithms, as described in:  R. J. McEliece. The Guruswami-Sudan
// Decoding Algorithm for Reed-Solomon Codes. IPN Progress Report
// 42-153, May 15, 2003.
// http://citeseer.ist.psu.edu/mceliece03guruswamisudan.html

NTL_CLIENT

// Note: In most places in this file, F (which is usually a field ZZ_p
// or GF2E) may also be a ring ZZ_p modulo a product of two large
// primes.

template <class FX>
struct RecoveryPoly {
    RecoveryPoly(vector<nservers_t> G, FX phi) : G(G), phi(phi) {}

    vector<nservers_t> G;
    FX phi;
};

template <class FX>
struct RecoveryPolyMulti {
    RecoveryPolyMulti(vector<nservers_t> G, vector<FX> phis) : G(G), phis(phis) {}

    vector<nservers_t> G;
    vector<FX> phis;
};

template <class F>
struct DecoderResult {
    DecoderResult (vector<nservers_t> G, vector<map<dbsize_t, F> > recovered) : 
        G(G), recovered(recovered) {}
    
    vector<nservers_t> G;
    vector<map<dbsize_t, F> > recovered;
};


// A class for caching the results of C() computations.  Note that you
// need to keep a separate Ccache around for each different field you
// work in.

typedef pair<unsigned int, unsigned int> pairint;
#define Ccache_t map<pairint, F>

template <class F>
static F C(Ccache_t &Ccache, unsigned int n, unsigned int k)
{
    // Do we already know the answer to this one?
    pairint nk;
    nk.first = n;
    nk.second = k;
    typename Ccache_t::const_iterator Ci = Ccache.find(nk);
    if (Ci != Ccache.end()) {
    // Yes; return it
    return Ci->second;
    }
    
    // Calculate in ZZ in case we're working in a low-characteristic
    // field.
    ZZ num, dem;
    num = 1;
    dem = 1;
    unsigned int i;

    for (i=0;i<k;++i) {
    num *= (n-i);
    dem *= (k-i);
    }

    num /= dem;

    // Now convert it into the right field, cache the answer, and return
    // it
    F res;
    conv(res, num);
    Ccache[nk] = res;
    return res;
}


// Comparison functions used to let polynomials be keys for maps
bool cmp (const GF2EX& a, const GF2EX& b);
bool cmp (const ZZ_pX& a, const ZZ_pX& b);


// It's been tested on ZZ_p and GF2E.  I haven't tried it over
// ZZ_pE, but I don't know why it wouldn't work there too; it should
// only be necessary to specialize DT (see below) appropriately.
template <class F, class vec_F, class FX, class FXY, class mat_F>
class RSDecoder {
    public:
    RSDecoder() {}

    RSDecoder(const ZZ &p1, const ZZ &p2): p1(p1), p2(p2) {}

    bool Recover (dbsize_t bytes_per_word, nservers_t t, nservers_t h,
        const vector<nservers_t> &goodservers, const vector<vector<vec_F> > &values,
        const vec_F &server_indices, vector<vector<DecoderResult<F> > > &results,
        vector<std::set<dbsize_t> > &decoded,
        const vector<vec_F> &vec_interp_indices,
        const nqueries_t multi_only = 0);

    static string append(const string &s, const F &wz,
        unsigned int bytes_per_word);

    static void test_interpolate(unsigned short t,
        const vec_F &values, const vec_F &indices,
        const vector<unsigned short> &I,
        const vector<unsigned short> &G,
        unsigned short &numagree, unsigned short &numdisagree,
        vector<unsigned short> &vecagree, FX &phi);

#ifdef TEST_RR
    public:
#else
    private:
#endif
    // Return a list of roots for y of the bivariate polynomial
    // P(x,y).  If degreebound >= 0, only return those roots with
    // degree <= degreebound.  This routine handles the case where F
    // is the integers mod p1*p2 as well as fields.  This routine
    // may also return some spurious values, so you need to validate
    // the answers it gives.
    //
    // This is an implementation of the Roth-Ruckenstein algorithm
    // for finding roots of bivariate polynomials, as explained in
    // section IX of:  R. J. McEliece. The Guruswami-Sudan Decoding
    // Algorithm for Reed-Solomon Codes. IPN Progress Report 42-153,
    // May 15, 2003.
    // http://citeseer.ist.psu.edu/mceliece03guruswamisudan.html
    vector<FX> findroots(const FXY &P, int degreebound);

#if defined(TEST_FINDPOLYS) || defined(TIME_FINDPOLYS)
    public:
#else
    private:
#endif

    // findpolys wrapper function
    vector< RecoveryPoly<FX> > findpolys(int k,
        unsigned int t, const vector<unsigned short> &goodservers,
        const vec_F &indices, const vec_F &shares, TestType &testType, DPType
        &dpType, int &gord);
    vector< RecoveryPoly<FX> > findpolys(int k,
        unsigned int t, const vector<unsigned short> &goodservers,
        const vec_F &indices, const vec_F &shares, TestType &testType)
    {
        DPType dptype; int gord;
        return findpolys(k, t, goodservers, indices, shares, testType, dptype,
                gord);
    }
    vector< RecoveryPoly<FX> > findpolys(int k,
        unsigned int t, const vector<unsigned short> &goodservers,
        const vec_F &indices, const vec_F &shares)
    {
        TestType testtype; DPType dptype; int gord;
        return findpolys(k, t, goodservers, indices, shares, testtype, dptype,
                gord);
    }

    // Portfolio algorithm for single-polynomial decoding
    vector< RecoveryPoly<FX> > findpolys_best(unsigned int n, int ell,
        unsigned int h, const vector<unsigned short> &goodservers,
        const vec_F &indices, const vec_F &shares, TestType &testType,
        DPType &dpType, int &gord);

    // Wrapper for multi-polynomial decoding
    vector< RecoveryPolyMulti<FX> > findpolys_multi(unsigned int k,
        unsigned int t, const vector<unsigned short> &goodservers,
        const vec_F &indices, const vector<vec_F> &shares, TestType testType = CH_TK1);

    // Use the K\"otter algorithm (below) to construct an
    // appropriate  bivariate polynomial and use the
    // Roth-Ruckenstein algorithm (above) to extract its roots.
    // Filter the results of the above to end up with just the
    // answers we're looking for.
    vector< RecoveryPoly<FX> > findpolys_gs(unsigned int k, unsigned int t,
        const vector<unsigned short> &goodservers,
        const vec_F& indices, const vec_F& shares, TestType testType = KOTTER);

    // Find all polynomials of degree at most k that agree with the
    // given (index,share) pairs at at least t of the n points.
    // This version simply uses brute force and interpolates all
    // subsets of size k+1 of the n points.  Note that in general,
    // this version does *not* run in polynomial time, but for some
    // cases (with n-k small, for example) it is faster than
    // Guruswami-Sudan.
    //
    vector<RecoveryPoly<FX> > findpolys_brute(
        int k_signed, unsigned int t,
        const vector<unsigned short> &goodservers,
        const vec_F &indices, const vec_F &shares);

    // Find a polynomial of degree ell that agrees with at least h shares.  Uses
    // the Berlekamp-Welch algorithm.  This algorithm is only guaranteed to work
    // when 2h > n + ell.
    vector<RecoveryPoly<FX> > findpolys_bw(
            unsigned int ell, unsigned int h,
            const vector<unsigned short> &goodservers,
            const vec_F &indices, const vec_F &shares);

    // A dynamic programming approach to decoding.  We can assume that some
    // servers are correct, some are wrong or a combination of the two
    vector<RecoveryPoly<FX> > findpolys_dp (unsigned int ell, unsigned int h,
            const vector<unsigned short> &goodservers,
            const vec_F &indices, const vec_F &shares,
            const DPType, unsigned int gord);

#ifdef INCLUDE_GENERAL_MULTI
    // (Needed for findpolys_ch_multi)
    vector<map<unsigned int, FX> > solve_partial_solution (
            unsigned int ell, unsigned int h, map<unsigned int, FX> partial,
            const vector<vector<FX> >& rows, vector<unsigned int> row_order,
            const vector<vector<unsigned int> >& order_of_expts);

    // The general Cohn-Heninger multi-polynomial decoding algorithm as
    // described in "Approximate common divisors via lattices." (2011)
    vector<RecoveryPolyMulti<FX> > findpolys_ch_multi(
            unsigned int ell, unsigned int h, unsigned int m, unsigned int& t,
            unsigned int& k,
            const vector<unsigned short> &goodservers,
            const vec_F &indices, const vector<vec_F> &shares);
#endif

    // The linear version of the Cohn-Heninger multi-polynomial decoding
    // algorithm as described in "Approximate common divisors via lattices."
    // (2011).  This is the general algorithm with t=k=1.
    vector<RecoveryPolyMulti<FX> > findpolys_ch_tk1(
            unsigned int ell, unsigned int h, unsigned int m,
            const vector<unsigned short> &goodservers,
            const vec_F &indices, const vector<vec_F> &shares);

    // Solve the linear system of equations Ax=b
    vector<FX> solve_linsys_FX (vector<vector<FX> >& A, vector<FX>& b,
            unsigned int modulus = 0);

    private:
    // A structure for holding derived types.  It needs to be
    // manually specialized for every kind of base field (ZZ_p,
    // GF2E, etc.) we may want to work in.
    struct DT {};

    // Construct a bivariate polynomial P(x,y) such that, for any
    // polynomial f(x) of degree at most v that agrees with at least
    // t of the given points, (y-f(x)) is a factor of P(x,y).  This
    // version is K"otter's Interpolation Algorithm, as described in
    // Section VII of R. J.  McEliece. The Guruswami-Sudan Decoding
    // Algorithm for Reed-Solomon Codes. IPN Progress Report 42-153,
    // May 15, 2003.
    // http://citeseer.ist.psu.edu/mceliece03guruswamisudan.html
    FXY interpolate_kotter(unsigned int v, unsigned int t,
        const vector<unsigned short> &goodservers,
        const vec_F &indices, const vec_F &shares);

private:
    // Construct a bivariate polynomial Q(x,y) such that, for any
    // polynomial g(x) of degree at most ell that agrees with at least
    // h of the given points, (y-g(x)) is a factor of Q(x,y).  This
    // version is Cohn and Heninger's improvement on the Guruswami-Sudan
    // Decoding Algorithm as described in sections 1.2, 2.2 and 4 of
    // Cohn, H. and Heninger, N.  Ideal forms of Coppersmith's theorem and
    // Guruswami-Sudan list decoding.  August 6, 2010.
    FXY interpolate_cohn_heninger(unsigned int v, unsigned int h,
        const vector<unsigned short> &goodservers,
        const vec_F &indices, const vec_F &shares);

    // Perform lattice basis reduction on a lattice with m basis vectors in F[z].
    // The lattice is represented by a vector 'lattice' in (F[z][x])^m such that
    // the coefficient vectors of each element in 'lattice' is a basis for the
    // lattice.  The result is a polynomial in F[z][x] whose coefficient vector
    // is the smallest vector in the lattice.  This implementation uses the
    // algorithm found in Mulders, T. and Storjohann, A.  On Lattice Reduction for
    // Polynomial Matrices.  December 2000.
    vector<vector<FX> > reduce_lattice_MS(vector<vector<FX> > lattice,
            // The auxiliary matrix must have the same number of rows as lattice
            // If not, no operations will be performed on aux.
            vector<vector<FX> >& aux,
            const int reduceTo = -1, const unsigned int reduceNumber = 1);

    FXY reduce_lattice_opt(vector<FXY> lattice, const unsigned int n,
            const unsigned int ell, const unsigned int k);

    // Evaluate the (r,s)th Hasse mixed partial derivative of g at
    // the point (alpha, beta), which is:
    //   \sum_{i,j} C(i,r) C(j,s) a_{i,j} alpha^{i-r} beta^{j-s}
    // = \sum_j C(j,s) beta^{j-s} \sum_i C(i,r) a_{i,j} alpha^{i-r}
    // where g(x,y) = \sum_{i,j} a_{i,j} x^i y^j
    // See page 14 of R. J. McEliece. The Guruswami-Sudan Decoding
    // Algorithm for Reed-Solomon Codes. IPN Progress Report 42-153,
    // May 15, 2003.
    // http://citeseer.ist.psu.edu/mceliece03guruswamisudan.html
    F evalhasse(const FXY &g, unsigned int r, unsigned int s,
        F alpha, F beta, Ccache_t &Ccache);

    // Depth-first search on the tree of coefficients; used in
    // Roth-Ruckenstein.
    void dfs(vector<FX> &res,
        int u, 
        vector<int> &pi, 
        vector<int> &Deg,
        vector<F> &Coeff, 
        vector<FXY> &Q, 
        int &t, 
        int degreebound);

    // Return a list of roots for y of the bivariate polynomial
    // P(x,y).  If degreebound >= 0, only return those roots with
    // degree <= degreebound.  This routine only works over fields F
    // (not rings).
    vector<FX> rr_findroots(const FXY &P, int degreebound);

    // In the case of ZZ_p, store the factors of the modulus p1*p2.
    // If it's a field, p1 should be 1 and p2 should be prime.
    ZZ p1, p2;
};

// The two types of RSDecoders I've tested.  These are the types you
// would use to create a RSDecoder; i.e.:
//
// RSDecoder_ZZ_p decoder(p1, p2);
//
// or
//
// RSDecoder_GF2E decoder();
//
// Remember to also call ZZ_p::init(p1*p2) or GF2E::init(modpoly) as
// appropriate before creating the RSDecoder.
typedef RSDecoder<ZZ_p, vec_ZZ_p, ZZ_pX, ZZ_pXY, mat_ZZ_p> RSDecoder_ZZ_p;
typedef RSDecoder<GF2E, vec_GF2E, GF2EX, GF2EXY, mat_GF2E> RSDecoder_GF2E;

// Comparison functions used to let polynomials be keys for maps
bool cmp (const GF2EX& a, const GF2EX& b);
bool cmp (const ZZ_pX& a, const ZZ_pX& b);

#include "recover.h"
#include "rsdecoder_impl.h"

#endif