// Create a tensor product of c1,c2. It is assumed that *this,c1,c2
// are defined relative to the same set of primes and plaintext space,
// and that *this DOES NOT point to the same object as c1,c2
void Ctxt::tensorProduct(const Ctxt& c1, const Ctxt& c2)
{
  // c1,c2 may be scaled, so multiply by the inverse scalar if needed
  long f = 1;
  if (c1.ptxtSpace>2) 
    f = rem(context.productOfPrimes(c1.primeSet),c1.ptxtSpace);
  if (f!=1) f = InvMod(f,c1.ptxtSpace);

  clear();                // clear *this, before we start adding things to it
  primeSet = c1.primeSet; // set the correct prime-set before we begin

  // The actual tensoring
  CtxtPart tmpPart(context, IndexSet::emptySet()); // a scratch CtxtPart
  for (size_t i=0; i<c1.parts.size(); i++) {
    CtxtPart thisPart = c1.parts[i];
    if (f!=1) thisPart *= f;
    for (size_t j=0; j<c2.parts.size(); j++) {
      tmpPart = c2.parts[j];
      // What secret key will the product point to?
      if (!tmpPart.skHandle.mul(thisPart.skHandle, tmpPart.skHandle))
	Error("Ctxt::tensorProduct: cannot multiply secret-key handles");

      tmpPart *= thisPart; // The element of the tensor product

      // Check if we already have a part relative to this secret-key handle
      long k = getPartIndexByHandle(tmpPart.skHandle);
      if (k >= 0) // found a matching part
	parts[k] += tmpPart;
      else
	parts.push_back(tmpPart);
    }
  }

  /* Compute the noise estimate as c1.noiseVar * c2.noiseVar * factor
   * where the factor depends on the handles of c1,c2. Specifically,
   * if the largest powerOfS in c1,c2 are n1,n2, respectively, then we
   * have factor = ((n1+n2) choose n2).
   */
  long n1=0,  n2=0;
  for (size_t i=0; i<c1.parts.size(); i++) // get largest powerOfS in c1
    if (c1.parts[i].skHandle.getPowerOfS() > n1)
      n1 = c1.parts[i].skHandle.getPowerOfS();
  for (size_t i=0; i<c2.parts.size(); i++) // get largest powerOfS in c2
    if (c2.parts[i].skHandle.getPowerOfS() > n2)
      n2 = c2.parts[i].skHandle.getPowerOfS();

  // compute ((n1+n2) choose n2)
  long factor = 1;
  for (long i=n1+1; i<=n1+n2; i++) factor *= i;
  for (long i=n2  ; i>1     ; i--) factor /= i;

  noiseVar = c1.noiseVar * c2.noiseVar * factor * context.zMStar.get_cM();
  if (f!=1) {
    // WARNING: the following line is written just so to prevent overflow
    noiseVar = (noiseVar*f)*f; // because every product was scaled by f
  }
}

// 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));
}

void InitFFTPrimeInfo(FFTPrimeInfo& info, long q, long w, long bigtab)
{
   double qinv = 1/((double) q);

   long mr = CalcMaxRoot(q);

   info.q = q;
   info.qinv = qinv;
   info.zz_p_context = 0;


   info.RootTable.SetLength(mr+1);
   info.RootInvTable.SetLength(mr+1);
   info.TwoInvTable.SetLength(mr+1);
   info.TwoInvPreconTable.SetLength(mr+1);

   long *rt = &info.RootTable[0];
   long *rit = &info.RootInvTable[0];
   long *tit = &info.TwoInvTable[0];
   mulmod_precon_t *tipt = &info.TwoInvPreconTable[0];

   long j;
   long t;

   rt[mr] = w;
   for (j = mr-1; j >= 0; j--)
      rt[j] = MulMod(rt[j+1], rt[j+1], q);

   rit[mr] = InvMod(w, q);
   for (j = mr-1; j >= 0; j--)
      rit[j] = MulMod(rit[j+1], rit[j+1], q);

   t = InvMod(2, q);
   tit[0] = 1;
   for (j = 1; j <= mr; j++)
      tit[j] = MulMod(tit[j-1], t, q);

   for (j = 0; j <= mr; j++)
      tipt[j] = PrepMulModPrecon(tit[j], q, qinv);

   info.bigtab = bigtab;
}

NTL_CLIENT

#include <algorithm>   // defines count(...), min(...)
#include <iostream>

#include "NumbTh.h"
#include "PAlgebra.h"


// Generate the representation of Z[X]/(Phi_m(X),2) for the odd integer m
void PAlgebraModTwo::init(unsigned m)
{
  if (m == zmStar.M()) return; // nothign to do

  ((PAlgebra&)zmStar).init(m); // initialize the structure of (Z/mZ)*, if needed
  if (zmStar.M()==0 || zmStar.NSlots()==0) return; // error in zmStar
  unsigned nSlots = zmStar.NSlots();

  // Next compute the factors Ft of Phi_m(X) mod 2, for all t \in T

  //  GF2X PhimXmod = to_GF2X(zmStar.PhimX()); // Phi_m(X) mod 2
  PhimXmod = to_GF2X(zmStar.PhimX()); // Phi_m(X) mod 2

  EDF(factors, PhimXmod, zmStar.OrdTwo()); // equal-degree factorization

  // It is left to order the factors according to their representatives

  GF2XModulus F1(factors[0]); // We arbitrarily choose factors[0] as F1
  for (unsigned i=1; i<nSlots; i++) {
    unsigned t = zmStar.ith_rep(i); // Ft is minimal polynomial of x^{1/t} mod F1
    unsigned tInv = rep(inv(to_zz_p(t))); // tInv = t^{-1} mod m
    GF2X X2tInv = PowerXMod(tInv,F1);     // X2tInv = X^{1/t} mod F1
    IrredPolyMod(factors[i], X2tInv, F1);
  }
  /* Debugging sanity-check #1: we should have Ft= GCD(F1(X^t),Phi_m(X))
  GF2XModulus Pm2(PhimXmod);
  for (i=1; i<nSlots; i++) {
    unsigned t = T[i];
    GF2X X2t = PowerXMod(t,PhimXmod);  // X2t = X^t mod Phi_m(X)
    GF2X Ft = GCD(CompMod(F1,X2t,Pm2),Pm2);
    if (Ft != factors[i]) {
      cout << "Ft != F1(X^t) mod Phi_m(X), t=" << t << endl;
      exit(0);
    }
  }*******************************************************************/

  // Compute the CRT coefficients for the Ft's
  crtCoeffs.SetLength(nSlots);
  for (unsigned i=0; i<nSlots; i++) {
    GF2X te = PhimXmod / factors[i]; // \prod_{j\ne i} Fj
    te %= factors[i];              // \prod_{j\ne i} Fj mod Fi
    InvMod(crtCoeffs[i], te, factors[i]);// \prod_{j\ne i} Fj^{-1} mod Fi
  }
}

// divVec[d] = m/p_d^{e_d}, powVec[d] = p^{e_d}
// computes invVec[d] = divVec[d]^{-1} mod powVec[d]
void computeInvVec(Vec<long>& invVec,
                   const Vec<long>& divVec, const Vec<long>& powVec)
{
  long k = divVec.length();
  invVec.SetLength(k);

  for (long d = 0; d < k; d++) {
    long t1 = divVec[d] % powVec[d];
    long t2 = InvMod(t1, powVec[d]);
    invVec[d] = t2;
  }
}

ZZ ASTDiv::eval(Environment &env) { 			
	pair<VarInfo,VarInfo> p = getTypes(env);
	VarInfo leftInfo = p.first;
	VarInfo rightInfo = p.second;
	// XXX: this will cause a problem if we ever try to use
	// it during constant propagation/substitution
	const Group* lGroup = env.groups.at(leftInfo.group);
	const Group* rGroup = env.groups.at(rightInfo.group);
	const Group* retGroup = (lGroup != 0 ? lGroup : rGroup);
	// integers can just be divided
	if ((lGroup == 0 && rGroup == 0) 
		|| (leftInfo.type == VarInfo::MODULUS 
			|| rightInfo.type == VarInfo::MODULUS)) {
		return lhs->eval(env) / rhs->eval(env);
	} else if (leftInfo.type == VarInfo::EXPONENT || 
			   rightInfo.type == VarInfo::EXPONENT) {
		// want to get LHS * inv(RHS)
		// this only works if we are in a group with known order
		if (retGroup->getType() == Group::TYPE_PRIME) {
			ZZ ord = retGroup->getOrder();
			ZZ right = InvMod(rhs->eval(env), ord);
			return lhs->eval(env) * right;
			// in an RSA group, only permit computing 1/x if the caller
			// knows the order of the group
		} else if(retGroup->getType() == Group::TYPE_RSA && 
				  (retGroup->getOrder() != 0)) {
			ZZ ord = retGroup->getOrder();
			ZZ right = InvMod(rhs->eval(env), ord);
			return lhs->eval(env) * right;
		} else {
			throw CashException(CashException::CE_PARSE_ERROR,
								"That operation is not permitted in an RSA group");
		}
	} else {
		assert(retGroup);
		ZZ mod = retGroup->getModulus();
		return MulMod(lhs->eval(env), InvMod(rhs->eval(env), mod), mod);
	}
}	

ZZ BankTool::identifyDoubleSpender(const Coin& coin1, const ZZ &tPrime2, 
								   const ZZ& rValue2) const {
	// Should probably check that the R values are different
	ZZ mod = coin1.getCashGroup()->getModulus();
	ZZ order = coin1.getCashGroup()->getOrder();
	ZZ t1 = coin1.getTPrime();
    ZZ t2 = tPrime2;
	ZZ r1 = coin1.getR();
    ZZ r2 = rValue2;

    if (r2 > r1) {
		NTL::swap(r2, r1);
		NTL::swap(t2, t1);
    }

	ZZ exp = InvMod(r1 - r2, order);
	ZZ num = PowerMod(t2, r1, mod);
	ZZ denom = InvMod(PowerMod(t1, r2, mod), mod);
	ZZ base = MulMod(num, denom, mod);
	ZZ publicKeyUser = PowerMod(base, exp, mod);
	
	return publicKeyUser;
}

    long solve(int *S, int *R) { // start = S, each = R
        for (int i = 0; i < m; ++i) {
            visited[i] = false;
            inv[S[i]] = i;
        }

        long ret = 0, MOD = 1;
        for (int i = 0; i < m; ++i) {
            if (visited[i]) continue;
            // print("--- loop ---\n");
            int loop = 0;
            for (int j = i; !visited[j]; j = R[j]) {
                // print("{} ", j);
                visited[j] = true;
                vs[loop] = S[j] % 2;
                vr[loop] = j % 2;
                ++loop;
            }
            // print("\n");
            TIME += loop;

            int mod, rem;
            tie(mod, rem) = match(loop);
            // for (int i = 0; i < loop; ++i) {
            //     print("{} {}\n", vs[i], vr[i]);
            // }
            // print("mod = {}, rem = {}\n", mod, rem);

            if (mod == 0) {
                return 1e18;
            }
            long g = GCD(MOD, mod);
            if (ret % g != rem % g) {
                return 1e18;
            }
            // print("mod = {}, rem = {}\n", mod, rem);
            MOD /= g, mod /= g;
            long t = InvMod(MOD, mod) * ((rem - ret) / g % mod + mod) % mod;
            // print("t = {}\n", t);
            assert((t * MOD * g + ret) % (mod * g) == rem);
            ret = t * MOD * g + ret;
            MOD = MOD * mod * g;
            // print("MOD = {}, REM = {}\n", MOD, ret);

            // print("loop: {}\n", n);
        }
        return ret;
    }

// Division by constant
// FIXME: this is not alias friendly
SingleCRT& SingleCRT::operator/=(const ZZ &num)
{
  const IndexSet& s = map.getIndexSet();
  ZZ pi, n;

  for (long i = s.first(); i <= s.last(); i = s.next(i)) {
    pi = to_ZZ(context.ithPrime(i));
    rem(n,num,pi);
    InvMod(n,n,pi);   // n = num^{-1} mod pi

    vec_ZZ& vp = map[i].rep;
    for (long j=0; j<vp.length(); j++) MulMod(vp[j], vp[j], n, pi);
    map[i].normalize();
  }
  return *this;
}

// Процедура проверки подписи
void verifysign (ZZ &R, ZZ &r, ZZ &s, ZZ &e, Qxy &Q, Qxy &P, ZZ &q  )
{
    ZZ v, z1, z2;
    Qxy C;
    ZZ e1;
    if (InvModStatus( v, e, q ))
    {
        cout << "\nError in signature\n";
        goto end;
    }
    v = InvMod(e,q);
    z1 = (s*v)%q;
    z2 = (-r*v)%q;
    C = P*z1 + Q*z2;
    R = (conv<ZZ>(C.putx()))%q;

    cout << "\nv (dec) = \n" << v << endl;
    cout << "\nv (hex) = \n";
    show_dec_in_hex (v, L);
    cout << endl;

    cout << "\nz1 (dec) = \n" << z1 << endl;
    cout << "\nz1 (hex) = \n";
    show_dec_in_hex (z1, L);
    cout << endl;

    cout << "\nz2 (dec) = \n" << z2 << endl;
    cout << "\nz2 (hex) = \n";
    show_dec_in_hex (z2, L);
    cout << endl;

    cout << "\nPoint C:\n";
    C.putQxy();

    cout << "\nR (dec) = \n" << R << endl;
    cout << "\nR (hex) = \n";
    show_dec_in_hex (R, L);
    cout << endl;

end:
    if ( r == R )
        cout << "\nr = R\nSignature is OK";
    else
        cout << "\nSignature is FAILD";
    cout << endl;
}

// mod-switch down to primeSet \intersect s, after this call we have
// primeSet<=s. s must contain either all special primes or none of them.
void Ctxt::modDownToSet(const IndexSet &s)
{
  IndexSet intersection = primeSet & s;
  //  assert(!empty(intersection));       // some primes must be left
  if (empty(intersection)) {
    cerr << "modDownToSet called from "<<primeSet<<" to "<<s<<endl;
    exit(1);
  }
  if (intersection==primeSet) return; // nothing to do, removing no primes
  FHE_TIMER_START;

  IndexSet setDiff = primeSet / intersection; // set-minus

  // Scale down all the parts: use either a simple "drop down" (just removing
  // primes, i.e., reducing the ctxt modulo the samaller modulus), or a "real
  // modulus switching" with rounding, basically whichever yeilds smaller
  // noise. Recall that we keep the invariant that a ciphertext mod Q is
  // decrypted to Q*m (mod p), so if we just "drop down" we still need to
  // multiply by (Q^{-1} mod p).

  // Get an estimate for the added noise term for modulus switching
  xdouble addedNoiseVar = modSwitchAddedNoiseVar();
  if (noiseVar*ptxtSpace*ptxtSpace < addedNoiseVar) {     // just "drop down"
    long prodInv = InvMod(rem(context.productOfPrimes(setDiff),ptxtSpace), ptxtSpace);
    for (size_t i=0; i<parts.size(); i++) {
      parts[i].removePrimes(setDiff);         // remove the primes not in s
      parts[i] *= prodInv;
      // WARNING: the following line is written just so to prevent overflow
      noiseVar = noiseVar*prodInv*prodInv;
    }
    //    cerr << "DEGENERATE DROP\n";
  } 
  else {                                      // do real mod switching
    for (size_t i=0; i<parts.size(); i++) 
      parts[i].scaleDownToSet(intersection, ptxtSpace);

    // update the noise estimate
    double f = context.logOfProduct(setDiff);
    noiseVar /= xexp(2*f);
    noiseVar += addedNoiseVar;
  }
  primeSet.remove(setDiff); // remove the primes not in s
  assert(verifyPrimeSet()); // sanity-check: ensure primeSet is still valid
  FHE_TIMER_STOP;
}

// Divide a cipehrtext by p, for plaintext space p^r, r>1. It is assumed
// that the ciphertext encrypts a polynomial which is zero mod p. If this
// is not the case then the result will not be a valid ciphertext anymore.
// As a side-effect, the plaintext space is reduced from p^r to p^{r-1}.
void Ctxt::divideByP()
{
  // Special case: if *this is empty then do nothing
  if (this->isEmpty()) return;

  long p = getContext().zMStar.getP();
  assert (ptxtSpace>p);

  // multiply all the parts by p^{-1} mod Q (Q=productOfPrimes)
  ZZ pInverse, Q;
  getContext().productOfPrimes(Q, getPrimeSet());
  InvMod(pInverse, conv<ZZ>(p), Q);
  for (size_t i=0; i<parts.size(); i++)
    parts[i] *= pInverse;

  noiseVar /= (p * (double)p);  // noise is reduced by a p factor
  ptxtSpace /= p;               // and so is the plaintext space
}

// generate all matrices of the form s(X^{g^i})->s(X) for generators g of
// Zm* /<2> and i<ord(g). If g has different orders in Zm* and Zm* /<2>
// then generate also matrices of the form s(X^{g^{-i}})->s(X)
void add1DMatrices(FHESecKey& sKey, long keyID)
{
  const FHEcontext &context = sKey.getContext();
  long m = context.zMStar.getM();

  // key-switching matrices for the automorphisms
  for (long i = 0; i < (long)context.zMStar.numOfGens(); i++) {
    for (long j = 1; j < (long)context.zMStar.OrderOf(i); j++) {
      long val = PowerMod(context.zMStar.ZmStarGen(i), j, m); // val = g^j
      // From s(X^val) to s(X)
      sKey.GenKeySWmatrix(1, val, keyID, keyID);
      if (!context.zMStar.SameOrd(i))
	// also from s(X^{1/val}) to s(X)
	sKey.GenKeySWmatrix(1, InvMod(val,m), keyID, keyID);
    }
  }
  sKey.setKeySwitchMap(); // re-compute the key-switching map
}

// 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;

  // A hack: record this automorphism rather than actually performing it
  if (isSetAutomorphVals()) { // defined in NumbTh.h
    recordAutomorphVal(k);
    return;
  }
  // Special case: if *this is empty then do nothing
  if (this->isEmpty()) return;

  // Sanity check: verify that k \in Zm*
  long m = context.zMStar.getM();
  k = mcMod(k, m);
  assert (context.zMStar.inZmStar(k));

  long keyID=getKeyID();
  if (!pubKey.isReachable(k,keyID)) {// must have key-switching matrices for it
    throw std::logic_error("no key-switching matrices for k="+std::to_string(k)
                           + ", keyID="+std::to_string(keyID));
  }

  if (!inCanonicalForm(keyID)) {     // Re-linearize the input, if needed
    reLinearize(keyID);
    assert (inCanonicalForm(keyID)); // ensure that re-linearization succeeded
  }

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

    // A hack: record this automorphism rather than actually performing it
    if (isSetAutomorphVals2()) { // defined in NumbTh.h
      recordAutomorphVal2(amt);
      return;
    }
    //cerr << "********* automorph " << amt << "\n";
    automorph(amt);
    reLinearize(keyID);
    k = MulMod(k, InvMod(amt,m), m);
  }
  FHE_TIMER_STOP;
}

// Assumes current zz_p modulus is p^r
// computes S = F^{-1} mod G via Hensel lifting
void InvModpr(zz_pX& S, const zz_pX& F, const zz_pX& G, long p, long r)
{
  ZZX ff, gg, ss, tt;

  ff = to_ZZX(F); 
  gg = to_ZZX(G);

  zz_pBak bak;
  bak.save();
  zz_p::init(p);

  zz_pX f, g, s, t;
  f = to_zz_pX(ff);
  g = to_zz_pX(gg);
  s = InvMod(f, g);
  t = (1-s*f)/g;
  assert(s*f + t*g == 1);
  ss = to_ZZX(s);
  tt = to_ZZX(t);

  ZZ pk = to_ZZ(1);

  for (long k = 1; k < r; k++) {
    // lift from p^k to p^{k+1}
    pk = pk * p;

    assert(divide(ss*ff + tt*gg - 1, pk));

    zz_pX d = to_zz_pX( (1 - (ss*ff + tt*gg))/pk );
    zz_pX s1, t1;
    s1 = (s * d) % g;
    t1 = (d-s1*f)/g;
    ss = ss + pk*to_ZZX(s1);
    tt = tt + pk*to_ZZX(t1);
  }

  bak.restore();

  S = to_zz_pX(ss);

  assert((S*F) % G == 1);
}

ZZ PaillierParty::decrypt(const std::map<uint32_t,ZZ> &partialCiphers,
                          const std::vector<ZZ> &pubKeys) {
    ZZ sum(0);
    for (auto &pubKey : pubKeys) {
        sum = AddMod(sum,pubKey,m_n);
    }

    ZZ phi = MulMod(m_a,sum,m_n);

    ZZ prod(1);
    for (auto &partialCipher: partialCiphers) {
        uint32_t partyId = partialCipher.first;
        prod = MulMod(prod,PowerMod(partialCipher.second, 2*m_lagrangeBasis[partyId],m_field),m_field);
    }

    ZZ LResult = L_function(prod);

    ZZ combinedShare = InvMod(MulMod(4*m_delta*m_delta,phi,m_n),m_n);
    return MulMod(LResult, combinedShare, m_n);
}

static void
init_representatives(Vec<long>& representatives, long dim, 
                     const Vec<long>& mvec, const PAlgebra& zMStar)
{
  assert(dim >= 0 && dim < mvec.length());

  // special case
  if (dim >= LONG(zMStar.numOfGens())) {
    representatives.SetLength(1);
    representatives[0] = 1;
    return;
  }
  
  long m = mvec[dim];
  long D = zMStar.OrderOf(dim);
  long g = InvMod(zMStar.ZmStarGen(dim) % m, m);

  representatives.SetLength(D);
  for (long i = 0; i < D; i++)
    representatives[i] = PowerMod(g, i, m);
}

static void
init_representatives(Vec<long>& representatives, long dim, 
                     const Vec<long>& mvec, const PAlgebra& zMStar)
{
  //OLD: assert(dim >= 0 && dim < mvec.length());
  helib::assertInRange(dim, 0l, mvec.length(), "Invalid argument: dim must be between 0 and mvec.length()");

  // special case
  if (dim >= LONG(zMStar.numOfGens())) {
    representatives.SetLength(1);
    representatives[0] = 1;
    return;
  }
  
  long m = mvec[dim];
  long D = zMStar.OrderOf(dim);
  long g = InvMod(zMStar.ZmStarGen(dim) % m, m);

  representatives.SetLength(D);
  for (long i = 0; i < D; i++)
    representatives[i] = PowerMod(g, i, m);
}

ZZ get_PowerMod(const ZZ& numar, const ZZ& exponent, const ZZ& baza)
{
	if (exponent == 0) 
		return to_ZZ(1);

	long k = NumBits(exponent);

	ZZ res;
	res = 1;

	for (long i = k-1; i >= 0; i--)
	{
		res = (res*res) % baza;
		if (bit(exponent, i) == 1)
			res = (res*numar) % n;
	}

	if (exponent < 0)
		return InvMod(res, baza);
	else
		return res;
}

static
void compute_a_vals(Vec<ZZ>& a, long p, long e)
// computes a[m] = a(m)/m! for m = p..(e-1)(p-1)+1,
// as defined by Chen and Han.
// a.length() is set to (e-1)(p-1)+2

{
   ZZ p_to_e = power_ZZ(p, e);
   ZZ p_to_2e = power_ZZ(p, 2*e);

   long len = (e-1)*(p-1)+2;

   ZZ_pPush push(p_to_2e);

   ZZ_pX x_plus_1_to_p = power(ZZ_pX(INIT_MONO, 1) + 1, p);
   ZZ_pX denom = InvTrunc(x_plus_1_to_p - ZZ_pX(INIT_MONO, p), len);
   ZZ_pX poly = MulTrunc(x_plus_1_to_p, denom, len);
   poly *= p;

   a.SetLength(len);

   ZZ m_fac(1);
   for (long m = 2; m < p; m++) {
      m_fac = MulMod(m_fac, m, p_to_2e);
   }

   for (long m = p; m < len; m++) {
      m_fac = MulMod(m_fac, m, p_to_2e);
      ZZ c = rep(coeff(poly, m));
      ZZ d = GCD(m_fac, p_to_2e);
      if (d == 0 || d > p_to_e || c % d != 0) Error("cannot divide");
      ZZ m_fac_deflated = (m_fac / d) % p_to_e;
      ZZ c_deflated = (c / d) % p_to_e;
      a[m] = MulMod(c_deflated, InvMod(m_fac_deflated, p_to_e), p_to_e);
   }

}