static void ZZ_pX_InvMod_newton_unram(struct ZZ_pX &x, const struct ZZ_pX &a, const struct ZZ_pXModulus &F, const struct ZZ_pContext &cpn, const struct ZZ_pContext &cp)
{
    //int j;
    cp.restore();
    ZZ_pX *amodp = new ZZ_pX();
    ZZ_pX *xmodp = new ZZ_pX();
    ZZ_pX *fmodp = new ZZ_pX();
    ZZ_pX_conv_modulus(*amodp, a, cp);
    ZZ_pX_conv_modulus(*fmodp, F.val(), cp);
    InvMod(*xmodp, *amodp, *fmodp);
    //cout << "xmodp: " << *xmodp << "\namodp: " << *amodp << "\nfmodp: " << *fmodp << "\n";
    cpn.restore();
    ZZ_pX *minusa = new ZZ_pX();
    ZZ_pX *xn = new ZZ_pX();
    ZZ_pX_conv_modulus(*xn, *xmodp, cpn);
    NTL::negate(*minusa, a);
    while (1 > 0)
    {
        // x_n = 2*x_{n-1} - a*x_{n-1}^2 = (2 - a*x_{n-1})*x_{n-1}
        MulMod(x, *minusa, *xn, F);
        SetCoeff(x, 0, ConstTerm(x) + 2);
        MulMod(x, x, *xn, F);
        if (x == *xn)
            break;
        *xn = x;
        //cout << "x: " << x << "\nxn: " << *xn << "\n";
        //cin >> j;
    }
    delete amodp;
    delete xmodp;
    delete fmodp;
    delete minusa;
    delete xn;
}

// The function compOrder(orders, classes,flag,m) computes the order of elements
// of the quotient group, relative to current equivalent classes. If flag==1
// then also check if the order is the same as in (Z/mZ)^* and store the order
// with negative sign if not.
static void 
compOrder(vector<long>& orders, vector<long>& classes, bool flag, long m)
{
  orders[0] = 0;
  orders[1] = 1;
  for (long i=2; i<m; i++) {
    if (classes[i] <= 1) { // ignore i not in Z_m^* and order-0 elements
      orders[i] = (classes[i]==1)? 1 : 0;
      continue;
    }

    // If not comparing order with (Z/mZ)^*, only compute the order of pivots

    if (!flag && classes[i]<i){          // not a pivot
      orders[i] = orders[classes[i]];
      continue;
    }

    // For an element i>1, the order is at least 2
    long j = MulMod(i, i, m);
    long ord = 2;
    while (classes[j] != 1) {
      j = MulMod(j, i, m); // next element in <i>
      ord++;    // count how many steps until we reach 1
    }

    // When we get here we have classes[j]==1, so if j!=1 it means that the
    // order of i in the quotient group is smaller than its order in the
    // entire group Z_m^*. If the flag is set then we store orders[i] = -ord.
    
    if (flag && j != 1) ord = -ord; // order in Z_m^* is larger than ord
    orders[i] = ord;
  }
}

static
void conjClasses(vector<unsigned long>& classes, unsigned long g, unsigned long m)
{
    for (unsigned long i=0; i<m; i++) {
        if (classes[i]==0) continue; // i \notin (Z/mZ)^*

        if (classes[i]<i) { // i is not a pivot, updated its pivot
            classes[i] = classes[classes[i]];
            continue;
        }

        // If i is a pivot, update other pivots to point to it
        unsigned long ii = i;
        unsigned long gg = g;
        unsigned long jj = MulMod(ii, gg, m);
        while (classes[jj] != i) {
            classes[classes[jj]]= i; // Merge the equivalence classes of j and i

            // Note: if classes[j]!=j then classes[j] will be updated later,
            //       when we get to i=j and use the code for "i not pivot".

            jj = MulMod(jj, g, m);
        }
    }
}

void PaillierParty::secretShare() {
    ZZ beta = getRandomInNStar(m_n);

    std::vector<ZZ> coefficients;

    coefficients.push_back(MulMod(beta,m_m,m_n*m_m));

    for (uint32_t i=1; i < m_numOfParties; i++) {
        coefficients.push_back(getRandomInNStar(m_n*m_m));
    }

    ZZ_p::init(m_n*m_m);
    ZZ_pX polynomial;
    for (uint32_t i=0; i < m_numOfParties; i++) {
        SetCoeff(polynomial, i, conv<ZZ_p>(coefficients[i]));
    }

    for (auto &party : m_parties) {
        ZZ result = rep(eval(polynomial,ZZ_p(party.first)));
        sendZZTo(result,party.second);
    }

    ZZ_p s_i = eval(polynomial,ZZ_p(m_partyId));
    for (auto &party : m_parties) {
        ZZ value;
        receiveZZFrom(value,party.second);
        ZZ_p coefficient = conv<ZZ_p>(value);
        s_i = s_i + coefficient;
    }

    m_share = rep(s_i);

    m_pubKey = MulMod(MulMod(m_a,beta,m_n),m_m,m_n);
}

// Sets the prime defining the field for the curve and stores certain values
void Icart::setPrime(ZZ* p)
{
    //ZZ_p::init(*p);
    // Icart hash function uses 1/3 root, which is equivalent to (2p-1)/3
    exp = MulMod( SubMod( MulMod(ZZ(2), *p, *p), ZZ(1), *p), InvMod(ZZ(3),*p), *p);
    // Store inverse values to be used later
    ts = inv(ZZ_p(27));
    th = inv(ZZ_p(3));
}

/*
 * Must guarantee c+c DO NOT OVERFLOW!!!(both a, b, c are INTEGERS)
 * $a or $b may be negative, however $c must be positive
 */
template<class T> T PowMod( T a, T b, T c) {
   	T r=Mod((T)1,c);
   	a=Mod(a,c);
	while(b != 0) {
		if(b & 1) r=MulMod(r, a, c);
		a = MulMod( a, a, c);
		b >>= 1;
	}
	return r;
}

void Shares::addShares(map<string, ZZ> newShares){
    for(auto i : newShares){
        ZZ tmp = MulMod(i.second, shares[i.first], groupModulus);
        shares[i.first] = tmp;
    }
    nbrShares++;
}

void InnerProduct(zz_p& x, const vec_zz_p& a, const vec_zz_p& b,
                  long offset)
{
   if (offset < 0) LogicError("InnerProduct: negative offset");
   if (NTL_OVERFLOW(offset, 1, 0)) ResourceError("InnerProduct: offset too big");

   long n = min(a.length(), b.length()+offset);
   long i;

   long accum, t;
   long p = zz_p::modulus();
   mulmod_t pinv = zz_p::ModulusInverse();


   const zz_p *ap = a.elts();
   const zz_p *bp = b.elts();

   accum = 0;
   for (i = offset; i < n; i++) {
      t = MulMod(rep(ap[i]), rep(bp[i-offset]), p, pinv);
      accum = AddMod(accum, t, p);
   }

   x.LoopHole() = accum;
}

// Apply F(X)->F(X^k) followed by re-liearization. The automorphism is possibly
// evaluated via a sequence of steps, to ensure that we can re-linearize the
// result of every step.
void Ctxt::smartAutomorph(long k) 
{
  FHE_TIMER_START;
  // Special case: if *this is empty then do nothing
  if (this->isEmpty()) return;

  long m = context.zMStar.getM();

  k = mcMod(k, m);

  // Sanity check: verify that k \in Zm*
  assert (context.zMStar.inZmStar(k));

  long keyID=getKeyID();
  if (!inCanonicalForm(keyID)) {     // Re-linearize the input, if needed
    reLinearize(keyID);
    assert (inCanonicalForm(keyID)); // ensure that re-linearization succeeded
  }
  assert (pubKey.isReachable(k,keyID)); // reachable from 1

  while (k != 1) {
    const KeySwitch& matrix = pubKey.getNextKSWmatrix(k,keyID);
    long amt = matrix.fromKey.getPowerOfX();

    automorph(amt);
    reLinearize(keyID);

    k = MulMod(k, InvMod(amt,m), m);
  }
  FHE_TIMER_STOP;
}

void build(zz_pXArgument& A, const zz_pX& h, const zz_pXModulus& F, long m)
{
   if (m <= 0 || deg(h) >= F.n) Error("build: bad args");

   if (m > F.n) m = F.n;

   long i;

   if (zz_pXArgBound > 0) {
      double sz = 1;
      sz = sz*F.n;
      sz = sz+6;
      sz = sz*(sizeof (long));
      sz = sz/1024;
      m = min(m, long(zz_pXArgBound/sz));
      m = max(m, 1);
   }

   zz_pXMultiplier M;

   build(M, h, F);

   A.H.SetLength(m+1);

   set(A.H[0]);
   A.H[1] = h;
   for (i = 2; i <= m; i++) 
      MulMod(A.H[i], A.H[i-1], M, F);
}

YASHE YASHE::readFromFile(std::string filename) {
  YASHE output;
  std::ifstream ifs(filename);
  boost::archive::text_iarchive ia(ifs);
  ia >> output;
  NTL::ZZ_p::init(output.cModulus);
  output.cycloMod = NTL::ZZ_pXModulus(NTL::conv<NTL::ZZ_pX>(output.cycloModX));
  {
    NTL::ZZ_pPush push(output.bigModulus); // switch to multiplication modulus
    // make another modulus for fast multiplication
    output.bigCycloMod = NTL::ZZ_pXModulus(NTL::conv<NTL::ZZ_pX>(output.cycloModX));
  }
  {
    NTL::ZZ_pPush push(output.bigPModulus); // switch to plain text modulus
    // Factor the cyclotomic polynomial modulo t
    // for batch encryption
    NTL::ZZ_pXModulus pModulusX;
    NTL::build(pModulusX, NTL::conv<NTL::ZZ_pX>(output.cycloModX));

    output.crtElements.resize(output.factors.size());
    NTL::ZZ_pX fInv, fInvInv;
    for (long i = 0; i < output.factors.size(); i++) {
      div(fInv, NTL::conv<NTL::ZZ_pX>(output.cycloModX), output.factors[i]);
      rem(fInvInv, fInv, output.factors[i]);
      InvMod(fInvInv, fInvInv, output.factors[i]);
      output.crtElements[i] = MulMod(fInv, fInvInv, pModulusX);
    }
  }
  return output;
}

void MinPolyMod(zz_pX& hh, const zz_pX& g, const zz_pXModulus& F, long m)
{
   zz_pX h, h1;
   long n = F.n;
   if (m < 1 || m > n) Error("MinPoly: bad args");

   /* probabilistically compute min-poly */

   ProbMinPolyMod(h, g, F, m);
   if (deg(h) == m) { hh = h; return; }
   CompMod(h1, h, g, F);
   if (IsZero(h1)) { hh = h; return; }

   /* not completely successful...must iterate */

   long i;

   zz_pX h2, h3;
   zz_pXMultiplier H1;
   vec_zz_p R(INIT_SIZE, n);

   for (;;) {
      R.SetLength(n);
      for (i = 0; i < n; i++) random(R[i]);
      build(H1, h1, F);
      UpdateMap(R, R, H1, F);
      DoMinPolyMod(h2, g, F, m-deg(h), R);

      mul(h, h, h2);
      if (deg(h) == m) { hh = h; return; }
      CompMod(h3, h2, g, F);
      MulMod(h1, h3, H1, F);
      if (IsZero(h1)) { hh = h; return; }
   }
}

static
void BuildMatrix(vec_GF2XVec& M, long n, const GF2EX& g, const GF2EXModulus& F,
                 long verbose)
{
   long i, j, m;
   GF2EX h;


   M.SetLength(n);
   for (i = 0; i < n; i++)
      M[i].SetSize(n, 2*GF2E::WordLength());

   set(h);
   for (j = 0; j < n; j++) {
      if (verbose && j % 10 == 0) cerr << "+";

      m = deg(h);
      for (i = 0; i < n; i++) {
         if (i <= m)
            M[i][j] = rep(h.rep[i]);
         else
            clear(M[i][j]);
      }

      if (j < n-1)
         MulMod(h, h, g, F);
   }

   for (i = 0; i < n; i++)
      add(M[i][i], M[i][i], 1);

}

// Compute the mapping between linear array and a hypercube corresponding
/// to a single generator tree
void ComputeOneGenMapping(Permut& genMap, const OneGeneratorTree& T)
{
  Vec<long> dims(INIT_SIZE, T.getNleaves());
  Vec<long> coefs(INIT_SIZE,T.getNleaves());
  for (long i=T.getNleaves()-1, leaf=T.lastLeaf(); i>=0;
                                i--, leaf=T.prevLeaf(leaf)) {
    dims[i] = T[leaf].getData().size;
    coefs[i] = T[leaf].getData().e;
  }

  // A representation of an integer with digits from dims
  Vec<long> rep(INIT_SIZE, T.getNleaves());
  for (long i=0; i<rep.length(); i++) rep[i]=0; // initialize to zero

  // initialize to all zero
  long sz = T[0].getData().size;
  genMap.SetLength(sz);
  for (long i=0; i<sz; i++) genMap[i]=0;

  // compute the permutation
  for (long i=1; i<sz; i++) {
    addOne(rep, dims); // representation of i in base dims
    for (long j=0; j<coefs.length(); j++) {
      long tmp = MulMod(rep[j], coefs[j], sz);
      genMap[i] = AddMod(genMap[i], tmp, sz);
    }
  }
}

void build(ZZ_pXArgument& A, const ZZ_pX& h, const ZZ_pXModulus& F, long m)
{
   if (m <= 0 || deg(h) >= F.n) LogicError("build: bad args");

   if (m > F.n) m = F.n;

   long i;

   if (ZZ_pXArgBound > 0) {
      double sz = ZZ_p::storage();
      sz = sz*F.n;
      sz = sz + NTL_VECTOR_HEADER_SIZE + sizeof(vec_ZZ_p);
      sz = sz/1024;
      m = min(m, long(ZZ_pXArgBound/sz));
      m = max(m, 1);
   }

   ZZ_pXMultiplier M;

   build(M, h, F);

   A.H.SetLength(m+1);

   set(A.H[0]);
   A.H[1] = h;
   for (i = 2; i <= m; i++) 
      MulMod(A.H[i], A.H[i-1], M, F);
}

void CompMod(ZZ_pX& x, const ZZ_pX& g, const ZZ_pXArgument& A, 
             const ZZ_pXModulus& F)
{
   if (deg(g) <= 0) {
      x = g;
      return;
   }


   ZZ_pX s, t;
   ZZVec scratch(F.n, ZZ_p::ExtendedModulusSize());

   long m = A.H.length() - 1;
   long l = ((g.rep.length()+m-1)/m) - 1;

   ZZ_pXMultiplier M;
   build(M, A.H[m], F);

   InnerProduct(t, g.rep, l*m, l*m + m - 1, A.H, F.n, scratch);
   for (long i = l-1; i >= 0; i--) {
      InnerProduct(s, g.rep, i*m, i*m + m - 1, A.H, F.n, scratch);
      MulMod(t, t, M, F);
      add(t, t, s);
   }

   x = t;
}

void CompMod(zz_pX& x, const zz_pX& g, const zz_pXArgument& A, 
             const zz_pXModulus& F)
{
   if (deg(g) <= 0) {
      x = g;
      return;
   }


   zz_pX s, t;
   vec_zz_p scratch(INIT_SIZE, F.n);

   long m = A.H.length() - 1;
   long l = ((g.rep.length()+m-1)/m) - 1;

   zz_pXMultiplier M;
   build(M, A.H[m], F);

   InnerProduct(t, g.rep, l*m, l*m + m - 1, A.H, F.n, scratch);
   for (long i = l-1; i >= 0; i--) {
      InnerProduct(s, g.rep, i*m, i*m + m - 1, A.H, F.n, scratch);
      MulMod(t, t, M, F);
      add(t, t, s);
   }

   x = t;
}

long CRT(vec_ZZ& gg, ZZ& a, const vec_zz_p& G)
{
   long n = gg.length();
   if (G.length() != n) Error("CRT: vector length mismatch");

   long p = zz_p::modulus();

   ZZ new_a;
   mul(new_a, a, p);

   long a_inv;
   a_inv = rem(a, p);
   a_inv = InvMod(a_inv, p);

   long p1;
   p1 = p >> 1;

   ZZ a1;
   RightShift(a1, a, 1);

   long p_odd = (p & 1);

   long modified = 0;

   long h;

   ZZ g;
   long i;
   for (i = 0; i < n; i++) {
      if (!CRTInRange(gg[i], a)) {
         modified = 1;
         rem(g, gg[i], a);
         if (g > a1) sub(g, g, a);
      }
      else
         g = gg[i];

      h = rem(g, p);
      h = SubMod(rep(G[i]), h, p);
      h = MulMod(h, a_inv, p);
      if (h > p1)
         h = h - p;

      if (h != 0) {
         modified = 1;

         if (!p_odd && g > 0 && (h == p1))
            MulSubFrom(g, a, h);
         else
            MulAddTo(g, a, h);
      }

      gg[i] = g;
   }

   a = new_a;

   return modified;
}

NTL_START_IMPL


long IterIrredTest(const GF2X& f)
{
   long df = deg(f);

   if (df <= 0) return 0;
   if (df == 1) return 1;

   GF2XModulus F;

   build(F, f);
   
   GF2X h;
   SetX(h);
   SqrMod(h, h, F);

   long i, d, limit, limit_sqr;
   GF2X g, X, t, prod;


   SetX(X);

   i = 0;
   g = h;
   d = 1;
   limit = 2;
   limit_sqr = limit*limit;

   set(prod);

   while (2*d <= df) {
      add(t, g, X);
      MulMod(prod, prod, t, F);
      i++;
      if (i == limit_sqr) {
         GCD(t, f, prod);
         if (!IsOne(t)) return 0;

         set(prod);
         limit++;
         limit_sqr = limit*limit;
         i = 0;
      }

      d = d + 1;
      if (2*d <= deg(f)) {
         SqrMod(g, g, F);
      }
   }

   if (i > 0) {
      GCD(t, f, prod);
      if (!IsOne(t)) return 0;
   }

   return 1;
}

NTL::ZZ_pX YASHE::keyGen() {
  /**
   * The secret key is computed as
   *
   *      f' <- X_key
   *      f = (t*f' + 1) mod q
   *      secretKey = f
   * 
   * Secret keys are randomly generated until
   * an invertible key f^-1 is found
   */
  NTL::ZZ_pX secretKey, secretKeyInv;
  long inverseStatus;
  do {
    secretKey = pModulus * randomKeyPoly() + 1;
    inverseStatus = InvModStatus(secretKeyInv, secretKey, cycloMod);
  } while (inverseStatus == 1);

  /**
   * The public key is computed by
   *
   *      g <- X_key
   *      h = (t*g*f^-1) mod q
   *      publicKey = h
   */
  publicKey = MulMod(randomKeyPoly(), secretKeyInv, cycloMod);
  publicKey *= pModulus;

  /** 
   * The evaluation key is computed by
   *
   *        e, s <- X_err
   *        gamma = (powersOfRadix(f) + e + h*s) mod q
   *        evaluationKey = gamma
   */
  std::vector<NTL::ZZ_pX> evalKey;
  powersOfRadix(evalKey, secretKey);
  evalKeyMult.resize(decompSize);
  for (long i = 0; i < decompSize; i++) {
    evalKey[i] += randomErrPoly();
    evalKey[i] += MulMod(publicKey, randomErrPoly(), cycloMod);
    NTL::build(evalKeyMult[i], evalKey[i], cycloMod);
  }

  return secretKey;
}