static void recursiveInterpolateMod(ZZX& poly, const vec_long& x, vec_long& y,
				    const vec_zz_p& xmod, vec_zz_p& ymod,
				    long p, long p2e)
{
  if (p2e<=1) { // recursion edge condition, mod-1 poly = 0
    clear(poly);
    return;
  }

  // convert y input to zz_p
  for (long j=0; j<y.length(); j++) ymod[j] = to_zz_p(y[j] % p);

  // a polynomial p_i s.t. p_i(x[j]) = i'th p-base digit of poly(x[j])
  zz_pX polyMod;
  interpolate(polyMod, xmod, ymod);    // interpolation modulo p
  ZZX polyTmp; conv(polyTmp, polyMod); // convert to ZZX

  // update ytmp by subtracting the new digit, then dividing by p
  for (long j=0; j<y.length(); j++) {
    y[j] -= polyEvalMod(polyTmp,x[j],p2e); // mod p^e
    if (y[j]<0) y[j] += p2e;
// if (y[j] % p != 0) {
//   cerr << "@@error (p2^e="<<p2e<<"): y["<<j<<"] not divisible by "<<p<< endl;
//   exit(0);
// }
    y[j] /= p;
  } // maybe it's worth optimizing above by using multi-point evaluation

  // recursive call to get the solution of poly'(x)=y mod p^{e-1}
  recursiveInterpolateMod(poly, x, y, xmod, ymod, p, p2e/p);

  // return poly = p*poly' + polyTmp
  poly *= p;
  poly += polyTmp;
}

void Cmodulus::FFT(vec_long &y, const ZZX& x) const
{
  FHE_TIMER_START;
  zz_pBak bak; bak.save();
  context.restore();

  zz_pX& tmp = Cmodulus::getScratch_zz_pX();
  { FHE_NTIMER_START(FFT_remainder);
    conv(tmp,x);      // convert input to zpx format
  }

  if (!ALT_CRT && zMStar->getPow2()) {
    // special case when m is a power of 2

    long k = zMStar->getPow2();
    long phim = (1L << (k-1));
    long dx = deg(tmp);
    long p = zz_p::modulus();

    const zz_p *powers_p = (*powers).rep.elts();
    const mulmod_precon_t *powers_aux_p = powers_aux.elts();

    y.SetLength(phim);
    long *yp = y.elts();

    zz_p *tmp_p = tmp.rep.elts();

    for (long i = 0; i <= dx; i++)
      yp[i] = MulModPrecon(rep(tmp_p[i]), rep(powers_p[i]), p, powers_aux_p[i]);
    for (long i = dx+1; i < phim; i++)
      yp[i] = 0;

    FFTFwd(yp, yp, k-1, *zz_pInfo->p_info);

    return;
  }
    

  zz_p rt;
  conv(rt, root);  // convert root to zp format

  BluesteinFFT(tmp, getM(), rt, *powers, powers_aux, *Rb); // call the FFT routine

  // copy the result to the output vector y, keeping only the
  // entries corresponding to primitive roots of unity
  y.SetLength(zMStar->getPhiM());
  long i,j;
  long m = getM();
  for (i=j=0; i<m; i++)
    if (zMStar->inZmStar(i)) y[j++] = rep(coeff(tmp,i));
}

static
void RandomBasisElt(ZZ_pX& g, const vec_long& D, const vec_ZZVec& M)
{
   ZZ t1, t2;

   long n = D.length();

   long i, j, s;

   g.rep.SetLength(n);

   vec_ZZ_p& v = g.rep;

   for (j = n-1; j >= 0; j--) {
      if (D[j] == -1)
         random(v[j]);
      else {
         i = D[j];

         // v[j] = sum_{s=j+1}^{n-1} v[s]*M[i,s]

         clear(t1);

         for (s = j+1; s < n; s++) {
            mul(t2, rep(v[s]), M[i][s]);
            add(t1, t1, t2);
         }

         conv(v[j], t1);
      }
   }

   g.normalize();
}

void Cmodulus::FFT_aux(vec_long &y, zz_pX& tmp) const
{

  if (zMStar->getPow2()) {
    // special case when m is a power of 2

    long k = zMStar->getPow2();
    long phim = (1L << (k-1));
    long dx = deg(tmp);
    long p = zz_p::modulus();

    const zz_p *powers_p = (*powers).rep.elts();
    const mulmod_precon_t *powers_aux_p = powers_aux.elts();

    y.SetLength(phim);
    long *yp = y.elts();

    zz_p *tmp_p = tmp.rep.elts();

    for (long i = 0; i <= dx; i++)
      yp[i] = MulModPrecon(rep(tmp_p[i]), rep(powers_p[i]), p, powers_aux_p[i]);
    for (long i = dx+1; i < phim; i++)
      yp[i] = 0;

#ifdef FHE_OPENCL
    AltFFTFwd(yp, yp, k-1, *altFFTInfo);
#else
    FFTFwd(yp, yp, k-1, *zz_pInfo->p_info);
#endif

    return;
  }
    

  zz_p rt;
  conv(rt, root);  // convert root to zp format

  BluesteinFFT(tmp, getM(), rt, *powers, powers_aux, *Rb); // call the FFT routine

  // copy the result to the output vector y, keeping only the
  // entries corresponding to primitive roots of unity
  y.SetLength(zMStar->getPhiM());
  long i,j;
  long m = getM();
  for (i=j=0; i<m; i++)
    if (zMStar->inZmStar(i)) y[j++] = rep(coeff(tmp,i));
}

// Interpolate the integer polynomial such that poly(x[i] mod p)=y[i] (mod p^e)
// It is assumed that the points x[i] are all distinct modulo p
void interpolateMod(ZZX& poly, const vec_long& x, const vec_long& y,
		    long p, long e)
{
  poly = ZZX::zero();       // initialize to zero
  long p2e = power_long(p,e); // p^e

  vec_long ytmp(INIT_SIZE, y.length()); // A temporary writable copy
  for (long j=0; j<y.length(); j++) {
    ytmp[j] = y[j] % p2e;
    if (ytmp[j] < 0) ytmp[j] += p2e;
  }

  zz_pBak bak; bak.save();    // Set the current modulus to p
  zz_p::init(p);

  vec_zz_p xmod(INIT_SIZE, x.length()); // convert to zz_p
  for (long j=0; j<x.length(); j++) xmod[j] = to_zz_p(x[j] % p);

  vec_zz_p ymod(INIT_SIZE, y.length()); // scratch space
  recursiveInterpolateMod(poly, x, ytmp, xmod, ymod, p, p2e);
}

void ResHalfGCD(zz_pX& U, zz_pX& V, vec_zz_p& cvec, vec_long& dvec)
{
   long d_red = (deg(U)+1)/2;

   if (IsZero(V) || deg(V) <= deg(U) - d_red) {
      return;
   }

   long du = deg(U);


   long d1 = (d_red + 1)/2;
   if (d1 < 1) d1 = 1;
   if (d1 >= d_red) d1 = d_red - 1;

   zz_pXMatrix M1;

   ResHalfGCD(M1, U, V, d1, cvec, dvec);
   mul(U, V, M1);

   long d2 = deg(V) - du + d_red;

   if (IsZero(V) || d2 <= 0) {
      return;
   }

   M1(0,0).kill();
   M1(0,1).kill();
   M1(1,0).kill();
   M1(1,1).kill();


   zz_pX Q;

   append(cvec, LeadCoeff(V));
   append(dvec, dvec[dvec.length()-1]-deg(U)+deg(V));
   DivRem(Q, U, U, V);
   swap(U, V);

   ResHalfGCD(M1, U, V, d2, cvec, dvec);

   mul(U, V, M1); 
}

void ResIterHalfGCD(ZZ_pXMatrix& M_out, ZZ_pX& U, ZZ_pX& V, long d_red,
                    vec_ZZ_p& cvec, vec_long& dvec)
{
   M_out(0,0).SetMaxLength(d_red);
   M_out(0,1).SetMaxLength(d_red);
   M_out(1,0).SetMaxLength(d_red);
   M_out(1,1).SetMaxLength(d_red);

   set(M_out(0,0));   clear(M_out(0,1));
   clear(M_out(1,0)); set(M_out(1,1));

   long goal = deg(U) - d_red;

   if (deg(V) <= goal)
      return;

   ZZVec tmp(deg(U)+1, ZZ_p::ExtendedModulusSize());
   ZZ_pX Q, t(INIT_SIZE, d_red);


   while (deg(V) > goal) {
      append(cvec, LeadCoeff(V));
      append(dvec, dvec[dvec.length()-1]-deg(U)+deg(V));
      PlainDivRem(Q, U, U, V, tmp);
      swap(U, V);

      mul(t, Q, M_out(1,0));
      sub(t, M_out(0,0), t);
      M_out(0,0) = M_out(1,0);
      M_out(1,0) = t;

      mul(t, Q, M_out(1,1));
      sub(t, M_out(0,1), t);
      M_out(0,1) = M_out(1,1);
      M_out(1,1) = t;
   }
}

void Cmodulus::iFFT(zz_pX &x, const vec_long& y)const
{
  FHE_TIMER_START;
  zz_pBak bak; bak.save();
  context.restore();

  if (zMStar->getPow2()) {
    // special case when m is a power of 2

    long k = zMStar->getPow2();
    long phim = (1L << (k-1));
    long p = zz_p::modulus();

    const zz_p *ipowers_p = (*ipowers).rep.elts();
    const mulmod_precon_t *ipowers_aux_p = ipowers_aux.elts();

    const long *yp = y.elts();

    vec_long& tmp = Cmodulus::getScratch_vec_long();
    tmp.SetLength(phim);
    long *tmp_p = tmp.elts();

#ifdef FHE_OPENCL
    AltFFTRev1(tmp_p, yp, k-1, *altFFTInfo);
#else
    FFTRev1(tmp_p, yp, k-1, *zz_pInfo->p_info);
#endif

    x.rep.SetLength(phim);
    zz_p *xp = x.rep.elts();

    for (long i = 0; i < phim; i++)
      xp[i].LoopHole() = MulModPrecon(tmp_p[i], rep(ipowers_p[i]), p, ipowers_aux_p[i]);


    x.normalize();

    return;
  }



  zz_p rt;
  long m = getM();

  // convert input to zpx format, initializing only the coeffs i s.t. (i,m)=1
  x.rep.SetLength(m);
  long i,j;
  for (i=j=0; i<m; i++)
    if (zMStar->inZmStar(i)) x.rep[i].LoopHole() = y[j++]; // DIRT: y[j] already reduced
  x.normalize();
  conv(rt, rInv);  // convert rInv to zp format

  BluesteinFFT(x, m, rt, *ipowers, ipowers_aux, *iRb); // call the FFT routine

  // reduce the result mod (Phi_m(X),q) and copy to the output polynomial x
  { FHE_NTIMER_START(iFFT_division);
    rem(x, x, *phimx); // out %= (Phi_m(X),q)
  }

  // normalize
  zz_p mm_inv;
  conv(mm_inv, m_inv);
  x *= mm_inv; 
}

inline void clear(vec_long& v) {
  memset(v.elts(),0,sizeof(long)*v.length());
}

static
void NullSpace(long& r, vec_long& D, vec_ZZVec& M, long verbose)
{
   long k, l, n;
   long i, j;
   long pos;
   ZZ t1, t2;
   ZZ *x, *y;

   const ZZ& p = ZZ_p::modulus();

   n = M.length();

   D.SetLength(n);
   for (j = 0; j < n; j++) D[j] = -1;

   r = 0;

   l = 0;
   for (k = 0; k < n; k++) {

      if (verbose && k % 10 == 0) cerr << "+";

      pos = -1;
      for (i = l; i < n; i++) {
         rem(t1, M[i][k], p);
         M[i][k] = t1;
         if (pos == -1 && !IsZero(t1))
            pos = i;
      }

      if (pos != -1) {
         swap(M[pos], M[l]);

         // make M[l, k] == -1 mod p, and make row l reduced

         InvMod(t1, M[l][k], p);
         NegateMod(t1, t1, p);
         for (j = k+1; j < n; j++) {
            rem(t2, M[l][j], p);
            MulMod(M[l][j], t2, t1, p);
         }

         for (i = l+1; i < n; i++) {
            // M[i] = M[i] + M[l]*M[i,k]

            t1 = M[i][k];   // this is already reduced

            x = M[i].elts() + (k+1);
            y = M[l].elts() + (k+1);

            for (j = k+1; j < n; j++, x++, y++) {
               // *x = *x + (*y)*t1

               mul(t2, *y, t1);
               add(*x, *x, t2);
            }
         }

         D[k] = l;   // variable k is defined by row l
         l++;

      }
      else {
         r++;
      }
   }
}

{
  pk.clear();
  //  cerr << "FHEPubKey[";
  seekPastChar(str, '['); // defined in NumbTh.cpp

  // sanity check, verify that basic context parameters are correct
  unsigned long m, p, r;
  vector<long> gens, ords;
  readContextBase(str, m, p, r, gens, ords);
  assert(comparePAlgebra(pk.getContext().zMStar, m, p, r, gens, ords));

  // Get the public encryption key itself
  str >> pk.pubEncrKey;

  // Get the vector of secret-key Hamming-weights
  vec_long vl;
  str >> vl;
  pk.skHwts.resize(vl.length());
  for (long i=0; i<(long)pk.skHwts.size(); i++) pk.skHwts[i] = vl[i];

  // Get the key-switching matrices
  long nMatrices;
  str >> nMatrices;
  pk.keySwitching.resize(nMatrices);
  for (long i=0; i<nMatrices; i++)  // read the matrix from input str
    pk.keySwitching[i].readMatrix(str, pk.getContext());

  // Get the key-switching map
  Vec< Vec<long> > vvl;
  str >> vvl;
  pk.keySwitchMap.resize(vvl.length());

  DynamicPtxtPowers(long _x, long _p, long nPowers, long _d=1) : p(_p)
  {
    assert (_x>=0 && _p>1 && nPowers>0); // Sanity check
    v.SetLength(nPowers);
    dpth.SetLength(nPowers);
    for (long i=1; i<nPowers; i++) // Initializes nPowers empty slots
      dpth[i]=v[i]=-1;

    v[0] = _x % p;         // store X itself in v[0]
    dpth[0] = _d;
    //    cerr << " initial power="<<v[0]<<", initial depth="<<dpth[0];
  }

// Returns the e'th power, computing it as needed
long DynamicPtxtPowers::getPower(long e)
{
  // FIXME: Do we want to allow the vector to grow? If so then begin by
  // checking e<v.length() and resizing if not. Currently throws an exception.

  if (v.at(e-1)<0) { // Not computed yet, compute it now
    
    long k = 1L<<(NextPowerOfTwo(e)-1); // largest power of two smaller than e

    v[e-1] = getPower(e-k);             // compute X^e = X^{e-k} * X^k
    v[e-1] = MulMod(v[e-1], getPower(k), p);
    dpth[e-1] = max(getDepth(k),getDepth(e-k)) +1;
    nMults++;
  }
  return v[e-1];
}

 long size() const { return v.length(); }

static 
void mul_aux(mat_zz_p& X, const mat_zz_p& A, const mat_zz_p& B)  
{  
   long n = A.NumRows();  
   long l = A.NumCols();  
   long m = B.NumCols();  
  
   if (l != B.NumRows())  
      Error("matrix mul: dimension mismatch");  
  
   X.SetDims(n, m); 

   if (m > 1) {  // new preconditioning code

      long p = zz_p::modulus();
      double pinv = zz_p::ModulusInverse();

      mul_aux_vec.SetLength(m);
      long *acc = mul_aux_vec.elts();

      long i, j, k;

      for (i = 0; i < n; i++) {
         const zz_p* ap = A[i].elts();

         for (j = 0; j < m; j++) acc[j] = 0;

         for (k = 0;  k < l; k++) {   
            long aa = rep(ap[k]);
            if (aa != 0) {
               const zz_p* bp = B[k].elts();
               long T1;
               mulmod_precon_t aapinv = PrepMulModPrecon(aa, p, pinv);

               for (j = 0; j < m; j++) {
        	  T1 = MulModPrecon(rep(bp[j]), aa, p, aapinv);
        	  acc[j] = AddMod(acc[j], T1, p);
               } 
            }
         }

         zz_p *xp = X[i].elts();
         for (j = 0; j < m; j++)
            xp[j].LoopHole() = acc[j];    
      }

   }
   else {  // just use the old code, w/o preconditioning

      long p = zz_p::modulus();
      double pinv = zz_p::ModulusInverse();

      long i, j, k;  
      long acc, tmp;  

      for (i = 1; i <= n; i++) {  
	 for (j = 1; j <= m; j++) {  
            acc = 0;  
            for(k = 1; k <= l; k++) {  
               tmp = MulMod(rep(A(i,k)), rep(B(k,j)), p, pinv);  
               acc = AddMod(acc, tmp, p);  
            }  
            X(i,j).LoopHole() = acc;  
	 } 
      }
  
   }
}  

 bool wasComputed(long i) { return (i>=0 && i<v.length() && v[i]>=0); }

 long getDepth(long e) { return (e>0)? dpth.at(e-1) : 0; }

void mul(vec_zz_p& x, const vec_zz_p& a, const mat_zz_p& B)
{
   long l = a.length();
   long m = B.NumCols();
  
   if (l != B.NumRows())  
      Error("matrix mul: dimension mismatch");  

   if (m == 0) { 

      x.SetLength(0);
      
   }
   else if (m == 1) {

      long p = zz_p::modulus();
      double pinv = zz_p::ModulusInverse();

      long acc, tmp;
      long k;

      acc = 0;  
      for(k = 1; k <= l; k++) {  
         tmp = MulMod(rep(a(k)), rep(B(k,1)), p, pinv);  
         acc = AddMod(acc, tmp, p);  
      } 

      x.SetLength(1);
      x(1).LoopHole()  = acc;
          
   }
   else {  // m > 1.  precondition


      long p = zz_p::modulus();
      double pinv = zz_p::ModulusInverse();

      mul_aux_vec.SetLength(m);
      long *acc = mul_aux_vec.elts();

      long j, k;


      const zz_p* ap = a.elts();

      for (j = 0; j < m; j++) acc[j] = 0;

      for (k = 0;  k < l; k++) {
         long aa = rep(ap[k]);
         if (aa != 0) {
            const zz_p* bp = B[k].elts();
            long T1;
            mulmod_precon_t aapinv = PrepMulModPrecon(aa, p, pinv);

            for (j = 0; j < m; j++) {
               T1 = MulModPrecon(rep(bp[j]), aa, p, aapinv);
               acc[j] = AddMod(acc[j], T1, p);
            }
         } 
      }

      x.SetLength(m);
      zz_p *xp = x.elts();
      for (j = 0; j < m; j++)
         xp[j].LoopHole() = acc[j];    
      
   }
}