#include "FHE/NTL-Subs.h"
#include "Math/Setup.h"

#include "Math/gfp.h"
#include "Math/gf2n.h"

#include "FHE/P2Data.h"
#include "FHE/QGroup.h"
#include "FHE/NoiseBounds.h"

#include "Tools/mkpath.h"

#include "FHEOffline/Proof.h"

#include <fstream>
using namespace std;

#include <NTL/ZZ.h>
#include <NTL/ZZX.h>
#include <NTL/GF2EXFactoring.h>
#include <NTL/GF2XFactoring.h>
NTL_CLIENT

#include "FHEOffline/DataSetup.h"


template <>
void generate_setup(int n_parties, int plaintext_length, int sec,
    FHE_Params& params, FFT_Data& FTD, int slack, bool round_up)
{
  Ring Rp;
  bigint p0p,p1p,p;
  SPDZ_Data_Setup_Char_p(Rp, FTD, p0p, p1p, n_parties, plaintext_length, sec,
      slack, round_up);
  params.set(Rp, {p0p, p1p});
}


template <>
void generate_setup(int n_parties, int plaintext_length, int sec,
    FHE_Params& params, P2Data& P2D, int slack, bool round_up)
{
  Ring R;
  bigint pr0,pr1;
  SPDZ_Data_Setup_Char_2(R, P2D, pr0, pr1, n_parties, plaintext_length, sec,
      slack, round_up);
  params.set(R, {pr0, pr1});
}


void generate_setup(int n, int lgp, int lg2, int sec, bool skip_2,
    int slack, bool round_up)
{
  DataSetup setup;

  // do the full setup for SHE data
  generate_setup(n, lgp, sec, setup.setup_p.params, setup.setup_p.FieldD, slack, round_up);
  if (!skip_2)
    generate_setup(n, lg2, sec, setup.setup_2.params, setup.setup_2.FieldD, slack, round_up);

  setup.write_setup(skip_2);
}


bool same_word_length(int l1, int l2)
{
  return l1 / 64 == l2 / 64;
}

template <>
int generate_semi_setup(int plaintext_length, int sec,
    FHE_Params& params, FFT_Data& FTD, bool round_up)
{
  int m = 1024;
  int lgp = plaintext_length;
  bigint p;
  generate_prime(p, lgp, m);
  int lgp0, lgp1;
  while (true)
    {
      SemiHomomorphicNoiseBounds nb(p, phi_N(m), 1, sec,
          numBits(NonInteractiveProof::slack(sec, phi_N(m))), true);
      bigint p1 = 2 * p * m, p0 = p;
      while (nb.min_p0(params.n_mults() > 0, p1) > p0)
        {
          p0 *= 2;
        }
      if (phi_N(m) < nb.min_phi_m(numBits(p0 * (params.n_mults() > 0 ? p1 : 1))))
        {
          m *= 2;
          generate_prime(p, lgp, m);
        }
      else
        {
          lgp0 = numBits(p0) + 1;
          lgp1 = numBits(p1) + 1;
          break;
        }
    }

  int extra_slack = common_semi_setup(params, m, p, lgp0, lgp1, round_up);

  FTD.init(params.get_ring(), p);
  gfp::init_field(p);
  return extra_slack;
}

template <>
int generate_semi_setup(int plaintext_length, int sec,
    FHE_Params& params, P2Data& P2D, bool round_up)
{
  if (params.n_mults() > 0)
    throw runtime_error("only implemented for 0-level BGV");
  gf2n_short::init_field(plaintext_length);
  int m;
  char_2_dimension(m, plaintext_length);
  SemiHomomorphicNoiseBounds nb(2, phi_N(m), 1, sec,
      numBits(NonInteractiveProof::slack(sec, phi_N(m))), true);
  int lgp0 = numBits(nb.min_p0(false, 0));
  int extra_slack = common_semi_setup(params, m, 2, lgp0, -1, round_up);
  load_or_generate(P2D, params.get_ring());
  return extra_slack;
}

int common_semi_setup(FHE_Params& params, int m, bigint p, int lgp0, int lgp1, bool round_up)
{
  cout << "Need ciphertext modulus of length " << lgp0;
  if (params.n_mults() > 0)
    cout << "+" << lgp1;
  cout << " and " << phi_N(m) << " slots" << endl;

  int extra_slack = 0;
  if (round_up)
    {
      int i;
      for (i = 0; i <= 20; i++)
        {
          if (SemiHomomorphicNoiseBounds::min_phi_m(lgp0 + i) > phi_N(m))
            break;
          if (not same_word_length(lgp0, lgp0 + i))
            break;
        }
      extra_slack = i - 1;
      lgp0 += extra_slack;
      cout << "Rounding up to " << lgp0 << ", giving extra slack of "
          << extra_slack << " bits" << endl;
    }

  Ring R;
  ::init(R, m);
  bigint p0, p1 = 1;
  if (params.n_mults() > 0)
  {
    generate_moduli(p0, p1, m, p, lgp0, lgp1);
    params.set(R, {p0, p1});
  }
  else
  {
    generate_modulus(p0, m, p, lgp0);
    params.set(R, {p0});
  }
  return extra_slack;
}

int finalize_lengths(int& lg2p0, int& lg2p1, int n, int m, int* lg2pi, bool round_up)
{
  if (n >= 2 and n <= 10)
    cout << "Difference to suggestion for p0: " << lg2p0 - lg2pi[n - 2]
        << ", for p1: " << lg2p1 - lg2pi[9 + n - 2] << endl;
  cout << "p0 needs " << int(ceil(1. * lg2p0 / 64)) << " words" << endl;
  cout << "p1 needs " << int(ceil(1. * lg2p1 / 64)) << " words" << endl;

  int extra_slack = 0;
  if (round_up)
    {
      int i = 0;
      for (i = 0; i < 10; i++)
        {
          if (phi_N(m) < NoiseBounds::min_phi_m(lg2p0 + lg2p1 + 2 * i))
            break;
          if (not same_word_length(lg2p0 + i, lg2p0))
            break;
          if (not same_word_length(lg2p1 + i, lg2p1))
            break;
        }
      i--;
      extra_slack = 2 * i;
      lg2p0 += i;
      lg2p1 += i;
      cout << "Rounding up to " << lg2p0 << "+" << lg2p1
          << ", giving extra slack of " << extra_slack << " bits" << endl;
    }

  cout << "Total length: " << lg2p0 + lg2p1 << endl;

  return extra_slack;
}



/******************************************************************************
 * Here onwards needs NTL
 ******************************************************************************/
 


/*
 * Subroutine for creating the FHE parameters
 */
int SPDZ_Data_Setup_Char_p_Sub(Ring& R, bigint& pr0, bigint& pr1, int n,
    int idx, int& m, bigint& p, int sec, int slack = 0, bool round_up = false)
{
  int lg2pi[5][2][9]
             = {  {  {130,132,132,132,132,132,132,132,132},
                     {104,104,104,106,106,108,108,110,110} },
                  {  {196,196,196,196,198,198,198,198,198},
                     {136,138,140,142,140,140,140,142,142} },
                  {  {325,325,325,325,330,330,330,330,330},
                     {205,205,205,210,205,205,205,205,210} },
                  {  {580,585,585,585,585,585,585,585,585},
                     {330,330,330,335,335,335,335,335,335} },
                  {  {1095,1095,1095,1095,1095,1095,1095,1095,1095},
                     {590,590,590,590,590,595,595,595,595} }
               };

  int lg2p0 = 0, lg2p1 = 0;
  if (n >= 2 and n <= 10)
    {
      lg2p0=lg2pi[idx][0][n-2];
      lg2p1=lg2pi[idx][1][n-2];
    }
  else if (sec == -1)
    throw runtime_error("no precomputed parameters available");

  while (sec != -1)
    {
      double phi_m_bound =
              NoiseBounds(p, phi_N(m), n, sec, slack).optimize(lg2p0, lg2p1);
      cout << "Trying primes of length " << lg2p0 << " and " << lg2p1 << endl;
      if (phi_N(m) < phi_m_bound)
        {
          int old_m = m;
          m = 2 << int(ceil(log2(phi_m_bound)));
          cout << "m = " << old_m << " too small, increasing it to " << m << endl;
          generate_prime(p, numBits(p), m);
        }
      else
        break;
    }

  init(R,m);
  int extra_slack = finalize_lengths(lg2p0, lg2p1, n, m, lg2pi[idx][0],
      round_up);
  generate_moduli(pr0, pr1, m, p, lg2p0, lg2p1);
  return extra_slack;
}

void generate_moduli(bigint& pr0, bigint& pr1, const int m, const bigint p,
        const int lg2p0, const int lg2p1)
{
  generate_modulus(pr0, m, p, lg2p0, "0");
  generate_modulus(pr1, m, p, lg2p1, "1", pr0);
}

void generate_modulus(bigint& pr, const int m, const bigint p, const int lg2pr,
    const string& i, const bigint& pr0)
{
  int ex;

  if (lg2pr==0) { throw invalid_params(); }

  bigint step=m;
  bigint twop=1<<(numBits(m)+1);
  bigint gc=gcd(step,twop);
  step=step*twop/gc;

  ex=lg2pr-numBits(p)-numBits(step)+1;
  if (ex<0) { cout << "Something wrong in lg2p" << i << " = " << lg2pr << endl; abort();  }
  pr=1; pr=(pr<<ex)*p*step+1;
  while (!probPrime(pr) || pr==pr0 || numBits(pr)<lg2pr) { pr=pr+p*step; }
  cout << "\t pr" << i << " = " << pr << "  :   " << numBits(pr) <<  endl;

  cout << "\t\tFollowing should be both 1" << endl;
  cout << "\t\tpr" << i << " mod m = " << pr%m << endl;
  cout << "\t\tpr" << i << " mod p = " << pr%p << endl;

  cout << "Minimal MAX_MOD_SZ = " << int(ceil(1. * lg2pr / 64)) << endl;
}

/*
 * Create the char p FHE parameters
 */
void SPDZ_Data_Setup_Char_p(Ring& R, FFT_Data& FTD, bigint& pr0, bigint& pr1,
    int n, int lgp, int sec, int slack, bool round_up)
{
  bigint p;
  int idx, m;
  SPDZ_Data_Setup_Primes(p, lgp, idx, m);
  SPDZ_Data_Setup_Char_p_Sub(R, pr0, pr1, n, idx, m, p, sec, slack, round_up);

  Zp_Data Zp(p);
  gfp::init_field(p);
  FTD.init(R,Zp);
}


/* Compute Phi(N) */
int phi_N(int N)
{
  int phiN=1,p,e;
  PrimeSeq s;
  while (N!=1)
    { p=s.next();
      e=0;
      while ((N%p)==0) { N=N/p; e++; }
      if (e!=0)
        { phiN=phiN*(p-1)*power_long(p,e-1); }
    }
  return phiN;
}


/* Compute mobius function (naive method as n is small) */
int mobius(int n)
{
  int p,e,arity=0;
  PrimeSeq s;
  while (n!=1)
    { p=s.next();
      e=0;
      while ((n%p)==0) { n=n/p; e++; }
      if (e>1) { return 0; }
      if (e!=0) { arity^=1; }
    }     
  if (arity==0) { return 1; }
  return -1;
}



/* Compute cyclotomic polynomial */
ZZX Cyclotomic(int N)
{
  ZZX Num,Den,G,F;
  NTL::set(Num); NTL::set(Den);
  int m,d;

  for (d=1; d<=N; d++)
    { if ((N%d)==0)
         { clear(G);
           SetCoeff(G,N/d,1); SetCoeff(G,0,-1);
           m=mobius(d);
           if (m==1)       { Num*=G; }
           else if (m==-1) { Den*=G; }
         }
    } 
  F=Num/Den;
  return F;
}


void init(Ring& Rg,int m)
{
  Rg.mm=m;
  Rg.phim=phi_N(Rg.mm);

  Rg.pi.resize(Rg.phim);    Rg.pi_inv.resize(Rg.mm);
  for (int i=0; i<Rg.mm; i++) { Rg.pi_inv[i]=-1; }

  int k=0;
  for (int i=1; i<Rg.mm; i++)
    { if (gcd(i,Rg.mm)==1)
        { Rg.pi[k]=i;
          Rg.pi_inv[i]=k;
          k++;
        }
    }

  ZZX P=Cyclotomic(Rg.mm);
  Rg.poly.resize(Rg.phim+1);
  for (int i=0; i<Rg.phim+1; i++)
    { Rg.poly[i]=to_int(coeff(P,i)); } 
}



// Computes a(b) mod c
GF2X Subs_Mod(const GF2X& a,const GF2X& b,const GF2X& c)
{
  GF2X ans,pb,bb=b%c;
  ans=to_GF2X(coeff(a,0));
  pb=bb;
  for (int i=1; i<=deg(a); i++)
    { ans=ans+pb*coeff(a,i);
      pb=MulMod(pb,bb,c);
    }
  return ans;
}


// Computes a(x^pow) mod c  where x^m=1
GF2X Subs_PowX_Mod(const GF2X& a,int pow,int m,const GF2X& c)
{
  GF2X ans; ans.SetMaxLength(m);
  for (int i=0; i<=deg(a); i++)
    { int j=MulMod(i,pow,m);
      if (IsOne(coeff(a,i))) { SetCoeff(ans,j,1); }
    }
  ans=ans%c;
  return ans;
}




void init(P2Data& P2D,const Ring& Rg)
{
  GF2X G,F;
  SetCoeff(G,gf2n_short::degree(),1);
  SetCoeff(G,0,1);
  for (int i=0; i<gf2n_short::get_nterms(); i++)
    { SetCoeff(G,gf2n_short::get_t(i),1); }
  //cout << "G = " << G << endl;

  for (int i=0; i<=Rg.phi_m(); i++)
    { if (((Rg.Phi()[i])%2)!=0)
        { SetCoeff(F,i,1); }
    }
  //cout << "F = " << F << endl;

  // seed randomness to achieve same result for all players
  // randomness is used in SFCanZass and FindRoot
  SetSeed(ZZ(0));
  
  // Now factor F modulo 2
  vec_GF2X facts=SFCanZass(F);

  // Check all is compatible
  int d=deg(facts[0]);
  if (d%deg(G)!=0)
    { throw invalid_params(); }

  // Compute the quotient group
  QGroup QGrp;
  int Gord=-1,e=Rg.phi_m()/d; // e = # of plaintext slots, phi(m)/degree
  int seed=1;
  while (Gord!=e)
    { QGrp.assign(Rg.m(),seed);       // QGrp encodes the the quotient group Z_m^*/<2>
      Gord=QGrp.order();
      if (Gord!=e) { cout << "Group order wrong, need to repeat the Haf-Mc algorithm" << endl; seed++; }
    }
  //cout << " l = " << Gord << " , d = " << d << endl;
  if ((Gord*gf2n_short::degree())!=Rg.phi_m())
    { cout << "Plaintext type requires Gord*gf2n_short::degree ==  phi_m" << endl;
      throw not_implemented();
    }

  vector<GF2X> Fi(Gord);
  vector<GF2X> Rts(Gord);
  vector<GF2X> u(Gord);
  /*
     Find map from Type 0 (mod G) -> Type 1 (mod Fi)
     for the first of the Fi's only
  */
  Fi[0]=facts[0];
  GF2E::init(facts[0]);     // work with the extension field GF_2[X]/Fi[0]
  GF2EX Ga=to_GF2EX(G);     // represent G as a polynomial over the extension field
  Rts[0]=rep(FindRoot(Ga)); // Find a roof of G in this field

  cout << "Fixing field ordering and the maps (Need to count to " << Gord << " here)\n\t";
  GF2E::init(G);
  GF2X g;
  vector<int> used(facts.length());
  for (int i=0; i<facts.length(); i++) { used[i]=0; }
  used[0]=1;
  for (int i=0; i<Gord; i++)
    { cout << i << " " << flush;
      if (i!=0)
        { int hpow=QGrp.nth_element(i);
          Rts[i]=Subs_PowX_Mod(Rts[0],hpow,Rg.m(),F);
          bool flag=false;
          for (int j=0; j<facts.length(); j++)
            { if (used[j]==0)
                { g=Subs_PowX_Mod(facts[0],hpow,Rg.m(),facts[j]);
                  if (IsZero(g))
                    { g=Subs_Mod(G,Rts[i],facts[j]);
                      if (!IsZero(g))
                        { cout << "Something wrong - G" << endl;
                          throw invalid_params();
                        }
                      Fi[i]=facts[j];
                      used[j]=1;
                      flag=true;
                      break;
                    }
                }
            }
          if (flag==false)
           { cout << "Something gone wrong" << endl;
             throw invalid_params();
           }
          Rts[i]=Rts[i]%Fi[i];
        }

     // Now sort out the projection map (for CRT reconstruction)
     GF2X te=(F/Fi[i]);
     GF2X tei=InvMod(te%Fi[i],Fi[i]);
     u[i]=MulMod(te,tei,F); // u[i] = \prod_{j!=i} F[j]*(F[j]^{-1} mod F[i])
   }
  cout << endl;

  // Make the forward matrix
  //   This is a deg(F) x (deg(G)*Gord)  matrix which maps elements
  //   vectors in the SIMD representation into plaintext vectors
  
  imatrix A;
  A.resize(Rg.phi_m(), imatrix::value_type(Gord*gf2n_short::degree()));
  P2D.A.resize(A[0].size(), imatrix::value_type(A.size()));
  for (int slot=0; slot<Gord; slot++)
    { for (int co=0; co<gf2n_short::degree(); co++)
        { // Work out how x^co in given slot maps to plaintext vector
          GF2X av;
          SetCoeff(av,co,1);
          // av is mod G, now map to mod Fi
          av=Subs_Mod(av,Rts[slot],Fi[slot]);
          // Now need to map using CRT to the plaintext vector
          av=MulMod(av,u[slot],F);
          //cout << slot << " " << co << " : " << av << endl;
          for (int k=0; k<Rg.phi_m(); k++)
	    {
              int i = slot*gf2n_short::degree()+co;
              if (IsOne(coeff(av,k)))
                { A[k][i]=1; }
              else
                { A[k][i]=0; }
	      P2D.A[i][k] = A[k][i];
	    }
       }
    }
  //cout << "Forward Matrix : " << endl; print(P2D.A);

  // Find pseudo inverse modulo 2
  pinv(P2D.Ai, A);
  P2D.Ai.resize(Gord*gf2n_short::degree());

  //cout << "Inverse Matrix : " << endl; print(P2D.Ai);

  P2D.slots=Gord;

}


/*
 * Create the FHE parameters
 */
void char_2_dimension(int& m, int& lg2)
{
  switch (lg2)
    { case -1:
        m=17;
        lg2=8;
        break;
      case 40:
        m=13175;
        break;
      case -40:
        m=5797;
        lg2=40;
        break;
      default:
        throw runtime_error("field size not supported");
        break;
    }
}

void SPDZ_Data_Setup_Char_2(Ring& R, P2Data& P2D, bigint& pr0, bigint& pr1,
    int n, int lg2, int sec, int slack, bool round_up)
{
  int lg2pi[2][9]
             = {  {70,70,70,70,70,70,70,70,70},
                  {70,75,75,75,75,80,80,80,80}
               };

  cout << "Setting up parameters\n";
  if ((n<2 || n>10) and sec == -1) { throw invalid_params(); }

  int m,lg2p0,lg2p1,ex;

  char_2_dimension(m, lg2);

  if (sec == -1)
  {
    lg2p0=lg2pi[0][n-2];
    lg2p1=lg2pi[1][n-2];
  }
  else
  {
    NoiseBounds(2, phi_N(m), n, sec, slack).optimize(lg2p0, lg2p1);
    finalize_lengths(lg2p0, lg2p1, n, m, lg2pi[0], round_up);
  }

  if (NoiseBounds::min_phi_m(lg2p0 + lg2p1) > phi_N(m))
    throw runtime_error("number of slots too small");

  cout << "m = " << m << endl;
  init(R,m);

  if (lg2p0==0 || lg2p1==0) { throw invalid_params(); }

  // We want pr0=pr1=1 mod (m*twop) where twop is the smallest 
  // power of two bigger than 2*m. This means we have m'th roots of
  // unity and twop'th roots of unity. This means FFT's are easier
  // to implement
  int lg2m=numBits(m);
  bigint step=m<<(lg2m+1);
  
  ex=lg2p0-2*lg2m;
  pr0=1; pr0=(pr0<<ex)*step+1;
  while (!probPrime(pr0)) { pr0=pr0+step; }
  cout << "\t pr0 = " << pr0 << "  :   " << numBits(pr0) << endl;

  ex=lg2p1-2*lg2m;
  pr1=1; pr1=(pr1<<ex)*step+1;
  while (!probPrime(pr1) || pr1==pr0) { pr1=pr1+step; }
  cout << "\t pr1 = " << pr1 << "  :   " << numBits(pr1) <<  endl;

  cout << "\t\tFollowing should be both 1" << endl;
  cout << "\t\tpr1 mod m = " << pr1%m << endl;
  cout << "\t\tpr1 mod 2^lg2m = " << pr1%(1<<lg2m) << endl;

  gf2n_short::init_field(lg2);
  load_or_generate(P2D, R);
}

void load_or_generate(P2Data& P2D, const Ring& R)
{
  try
  {
      P2D.load(R);
  }
  catch (...)
  {
      cout << "Loading failed" << endl;
      init(P2D,R);
      P2D.store(R);
  }
}



/*
 * Create FHE parameters for a general plaintext modulus p
 *   Basically this is for general large primes only
 */
void SPDZ_Data_Setup_Char_p_General(Ring& R, PPData& PPD, bigint& pr0,
        bigint& pr1, int n, int sec, bigint& p)
{
  cout << "Setting up parameters" << endl;

  int lgp=numBits(p);
  int mm,idx;

  // mm is the minimum value of m we will accept
  if (lgp<48)
    { mm=100;  // Test case
      idx=0;
    }
  else if (lgp <96)
    { mm=8192;
      idx=1;
    }
  else if (lgp<192)
    { mm=16384;
      idx=2;
    }
  else if (lgp<384)
    { mm=16384;
      idx=3;
    }
  else if (lgp<768)
    { mm=32768;
      idx=4;
    }
  else
    { throw invalid_params(); }

  // Now find the small factors of p-1 and their exponents
  bigint t=p-1;
  vector<long> primes(100),exp(100);

  PrimeSeq s;
  long pr;
  pr=s.next();
  int len=0;
  while (pr<2*mm)
    { int e=0;
      while ((t%pr)==0)
	{ e++;
	  t=t/pr;
        }
      if (e!=0)
	{ primes[len]=pr;
          exp[len]=e;
	  if (len!=0) { cout << " * "; }
          cout << pr << "^" << e << flush;
          len++;
        }
      pr=s.next();
    }
  cout << endl;

  // We want to find the best m which divides pr-1, such that
  //    - 2*m > phi(m) > mm
  //    - m has the smallest number of factors
  vector<int> ee;
  ee.resize(len);
  for (int i=0; i<len; i++) { ee[i]=0; }
  int min_hwt=-1,m=-1,bphi_m=-1,bmx=-1;
  bool flag=true;
  while (flag)
    { int cand_m=1,hwt=0,mx=0;
      for (int i=0; i<len; i++)
        { //cout << ee[i] << " ";
          if (ee[i]!=0)
            { hwt++;
              for (int j=0; j<ee[i]; j++)
                { cand_m*=primes[i]; }
              if (ee[i]>mx) { mx=ee[i]; }
            }
        }
      // Put "if" here to stop searching for things which will never work
      if (cand_m>1 && cand_m<4*mm)
        { //cout << " : " << cand_m << " : " << hwt << flush;
          int phim=phi_N(cand_m);
          //cout << " : " << phim << " : " << mm << endl;
          if (phim>mm && phim<3*mm)
	    { if (m==-1 || hwt<min_hwt || (hwt==min_hwt && mx<bmx) || (hwt==min_hwt && mx==bmx && phim<bphi_m))
	        { m=cand_m;
                  min_hwt=hwt;
                  bphi_m=phim;
                  bmx=mx;
                  //cout << "\t OK" << endl;
                }
            }
        }
      else
        { //cout << endl; 
        }
      int i=0;
      ee[i]=ee[i]+1;
      while (ee[i]>exp[i] && flag) 
	{ ee[i]=0; 
          i++;
          if (i==len) { flag=false; i=0; }
          else        { ee[i]=ee[i]+1;   }
	}
    }
  if (m==-1)
    { throw bad_value(); }
  cout << "Chosen value of m=" << m << "\t\t phi(m)=" << bphi_m << " : " << min_hwt << " : " << bmx << endl;

  SPDZ_Data_Setup_Char_p_Sub(R,pr0,pr1,n,idx,m,p,sec);
  int mx=0;
  for (int i=0; i<R.phi_m(); i++)
    { if (mx<R.Phi()[i]) { mx=R.Phi()[i]; } }
  cout << "Max Coeff = " << mx << endl;

  Zp_Data Zp(p);
  gfp::init_field(p);
  PPD.init(R,Zp);
}