void ECM(ZZ& q, const ZZ& N, long ncurves, long B1, long B2, long D, 
	 bool verbose) {
  // initialize ZZ_p
  ZZ_pBak bak;
  bak.save();
  ZZ_p::init(N);

  PrimeSeq seq;

  for (long i=0; i<ncurves; ++i) {
    if (verbose) {
      if (ncurves<NTL_MAX_LONG)
	std::cout<<"ECM: "<<NumBits(N)<<" bits; "
	    <<i<<"/"<<ncurves<<" curves\r"<<std::flush;
      else
	std::cout<<"ECM: "<<NumBits(N)<<" bits; "
	    <<i<<" curves\r"<<std::flush;
    }

    // try to find a factor
    ECM_one_curve(q,seq,B1,B2,D);
    if (!IsOne(q)) {
      if (verbose) {
	std::cout<<"ECM: "<<NumBits(N)<<" bits; "
             <<(i+1)<<" curves; found "<<q<<"  "<<std::endl;
      }
      return;
    }
  }
  if (verbose)
    std::cout<<"ECM: "<<NumBits(N)<<" bits; FAILED!        "<<std::endl;
  set(q);
}

void div(RR& z, const RR& a, const RR& b)
{
   if (IsZero(b))
      Error("RR: division by zero");

   if (IsZero(a)) {
      clear(z);
      return;
   }

   long la = NumBits(a.x);
   long lb = NumBits(b.x);

   long neg = (sign(a) != sign(b));

   long k = RR::prec - la + lb + 1;
   if (k < 0) k = 0;

   static RR t;
   static ZZ A, B, R;

   abs(A, a.x);
   LeftShift(A, A, k);

   abs(B, b.x);
   DivRem(t.x, R, A, B);

   t.e = a.e - b.e - k;

   normalize(z, t, !IsZero(R));

   if (neg)
      negate(z.x, z.x);
}

void CRSAFactorHintDlg::OnChangeGuessP() 
{
	UpdateData();
	if(m_choice==0){
		int base=(m_base*6+10)%20;
		CString GuessPString;
		if(m_msbLsb==0){
			GetDlgItemText(IDC_EDITGUESSP_LEFT,GuessPString);
			//	SetDlgItemText(IDC_EDITGUESSP_RIGHT,GuessPString);
		}
		else{
			GetDlgItemText(IDC_EDITGUESSP_RIGHT,GuessPString);
			//	SetDlgItemText(IDC_EDITGUESSP_LEFT,GuessPString);
		}
		ZZ tmpGuessP=setToStringValue(GuessPString,base);
		int bitsOfGuess=NumBits(tmpGuessP);
		int z=0;
		for(int i=0; i<GuessPString.GetLength()
			&&(GuessPString.GetAt(i)=='0'
			||GuessPString.GetAt(i)==' '); i++)
			if(GuessPString.GetAt(i)=='0')
				z++;
			int zeroBits=NumBits(power(to_ZZ(base),z)-1);
			bitsOfGuess+=zeroBits;
			m_GuessP=tmpGuessP;
			SetDlgItemInt(IDC_EDITB, bitsOfGuess);
	}
	UpdateData(false);
	//updateDim();
}

long GetPrimeBound(int deg_a, int deg_b, int bitsize_a, int bitsize_b, int num_of_additions){
	long res;
	int log_add;

	log_add = NumBits(num_of_additions)+1;
	res 	= 2 + NumBits(min(deg_a, deg_b)+1) + bitsize_a + bitsize_b + log_add;
	return res;
}

void exp(RR& res, const RR& x)
{
   if (x >= NTL_OVFBND || x <= -NTL_OVFBND)
      Error("RR: overflow");

   long p = RR::precision();

   // step 0: write x = n + f, n an integer and |f| <= 1/2
   //    careful -- we want to compute f to > p bits of precision


   RR f, nn;
   RR::SetPrecision(NTL_BITS_PER_LONG);
   round(nn, x);
   RR::SetPrecision(p + 10);
   sub(f, x, nn);
   long n = to_long(nn);

   // step 1: calculate t1 = e^n by repeated squaring

   RR::SetPrecision(p + NumBits(n) + 10);

   RR e;
   ComputeE(e);

   RR::SetPrecision(p + 10);

   RR t1;
   power(t1, e, n);

   // step 2: calculate t2 = e^f using Taylor series expansion

   RR::SetPrecision(p + NumBits(p) + 10);

   RR t2, s, s1, t;
   long i;

   s = 0;
   t = 1;

   for (i = 1; ; i++) {
      add(s1, s, t);
      if (s == s1) break;
      xcopy(s, s1);
      mul(t, t, f);
      div(t, t, i);
   }

   xcopy(t2, s);

   RR::SetPrecision(p);

   mul(res, t1, t2);
}

void totalSums(const EncryptedArray& ea, Ctxt& ctxt)
{
  long n = ea.size();

  if (n == 1) return;

  Ctxt orig = ctxt;

  long k = NumBits(n);
  long e = 1;

  for (long i = k-2; i >= 0; i--) {
    Ctxt tmp1 = ctxt;
    ea.rotate(tmp1, e);
    ctxt += tmp1; // ctxt = ctxt + (ctxt >>> e)
    e = 2*e;

    if (bit(n, i)) {
      Ctxt tmp2 = orig;
      ea.rotate(tmp2, e);
      ctxt += tmp2; // ctxt = ctxt + (orig >>> e)
                    // NOTE: we could have also computed
                    // ctxt =  (ctxt >>> e) + orig, however,
                    // this would give us greater depth/noise
      e += 1;
    }
  }
}

void RoundToZZ(ZZ& z, const RR& a)
{
   if (a.e >= 0) {
      LeftShift(z, a.x, a.e);
      return;
   }

   long len = NumBits(a.x);

   if (-a.e > len) {
      z = 0;
      return;
   }

   if (-a.e == len) {
      if (len == 1)
         z = 0;
      else
         z = sign(a.x);

      return;
   }

   NTL_TLS_LOCAL(RR, t);

   ConvPrec(t, a, len+a.e);

   LeftShift(z, t.x, t.e);
}

static
void normalize1(RR& z, const ZZ& y_x, long y_e, long prec, long residual)
{
   long len = NumBits(y_x);

   if (len > prec) {
      long correction = ZZ_RoundCorrection(y_x, len - prec, residual);

      RightShift(z.x, y_x, len - prec);

      if (correction) 
         add(z.x, z.x, correction);

      z.e = y_e + len - prec;
   }
   else if (len == 0) {
      clear(z.x);
      z.e = 0;
   }
   else {
      z.x = y_x;
      z.e = y_e;
   }

   if (!IsOdd(z.x))
      z.e += MakeOdd(z.x);

   if (z.e >= NTL_OVFBND)
      ResourceError("RR: overflow");

   if (z.e <= -NTL_OVFBND)
      ResourceError("RR: underflow");
}

void SqrRoot(RR& z, const RR& a)
{
   if (sign(a) < 0)
      ArithmeticError("RR: attempt to take square root of negative number");

   if (IsZero(a)) {
      clear(z);
      return;
   }

   RR t;
   ZZ T1, T2;
   long k;

   k = 2*RR::prec - NumBits(a.x) + 1;

   if (k < 0) k = 0;

   if ((a.e - k) & 1) k++;

   LeftShift(T1, a.x, k);
   // since k >= 2*prec - bits(a) + 1, T1 has at least 2*prec+1 bits,           
   // thus T1 >= 2^(2*prec)                                                     

   SqrRoot(t.x, T1); // t.x >= 2^prec thus t.x contains the round bit           
   t.e = (a.e - k)/2;
   sqr(T2, t.x);  

   // T1-T2 is the (lower part of the) sticky bit                               
   normalize(z, t, T2 < T1);
}

NTL_CLIENT
#include "FHE.h"
#include "timing.h"
#include "EncryptedArray.h"

#include <cassert>
#include <cstdio>


// computes ctxt^{2^d-1} using a method that takes
// O(log d) automorphisms and multiplications
void fastPower(Ctxt& ctxt, long d) 
{
  if (d == 1) return;

  Ctxt orig = ctxt;

  long k = NumBits(d);
  long e = 1;

  for (long i = k-2; i >= 0; i--) {
    Ctxt tmp1 = ctxt;
    tmp1.smartAutomorph(1L << e);
    ctxt.multiplyBy(tmp1);
    e = 2*e;

    if (bit(d, i)) {
      ctxt.smartAutomorph(2);
      ctxt.multiplyBy(orig);
      e += 1;
    }
  }
}

NTL_START_IMPL

static
long CharPolyBound(const mat_ZZ& a)
// This bound is computed via interpolation
// through complex roots of unity.

{
   long n = a.NumRows();
   long i;
   ZZ res, t1, t2;

   set(res);

   for (i = 0; i < n; i++) {
      InnerProduct(t1, a[i], a[i]);
      abs(t2, a[i][i]);
      mul(t2, t2, 2);
      add(t2, t2, 1);
      add(t1, t1, t2);
      if (t1 > 1) {
         SqrRoot(t1, t1);
         add(t1, t1, 1);
      }
      mul(res, res, t1);
   }

   return NumBits(res);
}

void power(mat_ZZ& X, const mat_ZZ& A, const ZZ& e)
{
    if (A.NumRows() != A.NumCols()) Error("power: non-square matrix");

    if (e == 0) {
        ident(X, A.NumRows());
        return;
    }

    mat_ZZ T1, T2;
    long i, k;

    k = NumBits(e);
    T1 = A;

    for (i = k-2; i >= 0; i--) {
        sqr(T2, T1);
        if (bit(e, i))
            mul(T1, T2, A);
        else
            T1 = T2;
    }

    if (e < 0)
        inv(X, T1);
    else
        X = T1;
}

void replicate0(const EncryptedArray& ea, Ctxt& ctxt, long pos)
{
  long dim = ea.dimension();

  for (long d = 0; d < dim; d++) {
    if (!ea.nativeDimension(d)) {
      long shamt = -ea.coordinate(d, pos);
      ea.rotate1D(ctxt, d, shamt, true); // "don't care"
    }

    Ctxt ctxt_orig = ctxt; 

    long sz = ea.sizeOfDimension(d);
    long k = NumBits(sz);
    long e = 1;

    // now process bits k-2 down to 0
    for (long j = k-2; j >= 0; j--) {
      // e -> 2*e
      Ctxt tmp = ctxt;
      ea.rotate1D(tmp, d, e, true); // "don't care"
      ctxt += tmp;
      e = 2*e;
      
      long b = bit(sz, j); // bit j of sz
      // e -> e+b
      if (b) {
        ea.rotate1D(ctxt, d, 1, true); // "don't care"
        ctxt += ctxt_orig;
        e++;
      }
    }
  }
}

void power(RR& z, const RR& a, long e)
{
   RR b, res;

   long n = NumBits(e);

   RRPush push;
   long p = RR::precision();
   RR::SetPrecision(p + n + 10);

   xcopy(b, a);

   set(res);
   long i;

   for (i = n-1; i >= 0; i--) {
      sqr(res, res);
      if (bit(e, i))
         mul(res, res, b);
   }

   RR::SetPrecision(p);

   if (e < 0) 
      inv(z, res);
   else
      xcopy(z, res);
}

void round(RR& z, const RR& a)
{
   if (a.e >= 0) {
      xcopy(z, a);
      return;
   }

   long len = NumBits(a.x);

   if (-a.e > len) {
      z = 0;
      return;
   }

   if (-a.e == len) {
      if (len == 1)
         z = 0;
      else
         z = sign(a.x);

      return;
   }

   NTL_TLS_LOCAL(RR, t);
   ConvPrec(t, a, len+a.e);
   xcopy(z, t);
}

void log(RR& res, const RR& x)
{
   if (x <= 0) Error("argument to log must be positive");

   long p = RR::precision();

   RR::SetPrecision(p + NumBits(p) + 10);

   RR y;
   long n;

   // re-write x = 2^n * (1 - y), where -1/2 < y < 1/4  (so 3/4 < 1-y < 3/2)

   if (x > 0.75 && x < 1.5) {
      n = 0;
      sub(y, 1, x);
   }
   else {
      n = Lg2(x) - 1;
      RR t;
      power2(t, -n);
      mul(t, t, x);
      while (t > 1.5) {
         mul(t, t, 0.5);
         n++;
      }

      sub(y, 1, t);
   }

   // compute s = - ln(1-y) by power series expansion

   RR s, s1, t, t1;

   s = 0;
   xcopy(t, y);
   xcopy(t1, y);

   long i;

   for (i = 2; ; i++) {
      add(s1, s, t);
      if (s == s1) break;
      xcopy(s, s1);
      mul(t1, t1, y);
      div(t, t1, i);
   }

   if (n == 0)
      t = 0;
   else {
      ComputeLn2(t);
      mul(t, t, n);
   }

   RR::SetPrecision(p);

   sub(res, t, s);
}

zz_pInfoT::zz_pInfoT(long NewP, long maxroot)
{
   if (maxroot < 0) LogicError("zz_pContext: maxroot may not be negative");

   if (NewP <= 1) LogicError("zz_pContext: p must be > 1");
   if (NumBits(NewP) > NTL_SP_NBITS) ResourceError("zz_pContext: modulus too big");

   ZZ P, B, M, M1, MinusM;
   long n, i;
   long q, t;

   p = NewP;

   pinv = 1/double(p);

   p_info = 0;

   conv(P, p);

   sqr(B, P);
   LeftShift(B, B, maxroot+NTL_FFTFudge);

   set(M);
   n = 0;
   while (M <= B) {
      UseFFTPrime(n);
      q = GetFFTPrime(n);
      n++;
      mul(M, M, q);
   }

   if (n > 4) LogicError("zz_pInit: too many primes");

   NumPrimes = n;
   PrimeCnt = n;
   MaxRoot = CalcMaxRoot(q);

   if (maxroot < MaxRoot)
      MaxRoot = maxroot;

   negate(MinusM, M);
   MinusMModP = rem(MinusM, p);

   CoeffModP.SetLength(n);
   x.SetLength(n);
   u.SetLength(n);

   for (i = 0; i < n; i++) {
      q = GetFFTPrime(i);

      div(M1, M, q);
      t = rem(M1, q);
      t = InvMod(t, q);
      if (NTL_zz_p_QUICK_CRT) mul(M1, M1, t);
      CoeffModP[i] = rem(M1, p);
      x[i] = ((double) t)/((double) q);
      u[i] = t;
   }
}

void KarSqr(ZZX& c, const ZZX& a)
{
   if (IsZero(a)) {
      clear(c);
      return;
   }

   vec_ZZ mem;

   const ZZ *ap;
   ZZ *cp;

   long sa = a.rep.length();

   if (&a == &c) {
      mem = a.rep;
      ap = mem.elts();
   }
   else
      ap = a.rep.elts();

   c.rep.SetLength(sa+sa-1);
   cp = c.rep.elts();

   long maxa, xover;

   maxa = MaxBits(a);

   xover = 2;

   if (sa < xover)
      PlainSqr(cp, ap, sa);
   else {
      /* karatsuba */

      long n, hn, sp, depth;

      n = sa;
      sp = 0;
      depth = 0;
      do {
         hn = (n+1) >> 1;
         sp += hn+hn+hn - 1;
         n = hn;
         depth++;
      } while (n >= xover);

      ZZVec stk;
      stk.SetSize(sp,
         ((2*maxa + NumBits(sa) + 2*depth + 10)
          + NTL_ZZ_NBITS-1)/NTL_ZZ_NBITS);

      KarSqr(cp, ap, sa, stk.elts());
   }

   c.normalize();
}

void pow(RR& res, const RR& x, const RR& y)
{
   if (y == 0) {
      res = 1;
      return;
   }

   if (x == 0) {
      res = 0;
      return;
   }

   if (x == 1) {
      res = 1;
      return;
   }

   if (x < 0) {
      Error("pow: sorry...first argument to pow must be nonnegative");
   }

   long p = RR::precision();

   // calculate working precison...one could use p + NTL_BITS_PER_LONG + 10,
   // for example, but we want the behaviour to be machine independent.
   // so we calculate instead a rough approximation to log |y log(x)|

   RR t1, t2;
   long k;

   if (x > 0.5 && x < 1.5) {
      xcopy(t1, x - 1);
      k = Lg2(t1);
   }
   else {
      k = NumBits(Lg2(x));
   }

   k += Lg2(y);

   if (k > NTL_BITS_PER_LONG+10) Error("RR: overflow");

   if (k < 0) k = 0;


   RR::SetPrecision(p + k + 10);

   t1 = y*log(x);

   RR::SetPrecision(p);

   t2 = exp(t1);

   res = t2;
}

void sub(RR& z, const RR& a, const RR& b)
{
   NTL_TLS_LOCAL(RR, t);

   if (IsZero(a.x)) {
      negate(z, b);
      return;
   }

   if (IsZero(b.x)) {
      xcopy(z, a);
      return;
   }

   if (a.e > b.e) {
      if (a.e-b.e - max(RR::prec-NumBits(a.x),0) >= NumBits(b.x) + 2)
         normalize(z, a, -sign(b));
      else {
         LeftShift(t.x, a.x, a.e-b.e);
         sub(t.x, t.x, b.x);
         t.e = b.e;
         xcopy(z, t);
      }
   }
   else if (a.e < b.e) {
      if (b.e-a.e - max(RR::prec-NumBits(b.x),0) >= NumBits(a.x) + 2) {
         normalize(z, b, -sign(a));
         negate(z.x, z.x);
      }
      else {
         LeftShift(t.x, b.x, b.e-a.e);
         sub(t.x, a.x, t.x);
         t.e = a.e;
         xcopy(z, t);
      }
   }
   else {
      sub(t.x, a.x, b.x);
      t.e = a.e;
      normalize(z, t);
   }
}