text:
#include <NTL/ZZ.h>
#include <NTL/vec_ZZ.h>
#include <NTL/Lazy.h>
#include <NTL/fileio.h>

NTL_START_IMPL

const ZZ& ZZ::zero()
{
   NTL_THREAD_LOCAL static ZZ z;
   return z;
}

const ZZ& ZZ_expo(long e)
{
   NTL_THREAD_LOCAL static ZZ expo_helper;
   conv(expo_helper, e);
   return expo_helper;
}

void AddMod(ZZ& x, const ZZ& a, long b, const ZZ& n)
{
   NTL_ZZRegister(B);
   conv(B, b);
   AddMod(x, a, B, n);
}

void SubMod(ZZ& x, const ZZ& a, long b, const ZZ& n)
{
   NTL_ZZRegister(B);
   conv(B, b);
   SubMod(x, a, B, n);
}

void SubMod(ZZ& x, long a, const ZZ& b, const ZZ& n)
{
   NTL_ZZRegister(A);
   conv(A, a);
   SubMod(x, A, b, n);
}

// ****** input and output

NTL_THREAD_LOCAL static long iodigits = 0;
NTL_THREAD_LOCAL static long ioradix = 0;

// iodigits is the greatest integer such that 10^{iodigits} < NTL_WSP_BOUND
// ioradix = 10^{iodigits}

static void InitZZIO()
{
   long x;

   x = (NTL_WSP_BOUND-1)/10;
   iodigits = 0;
   ioradix = 1;

   while (x) {
      x = x / 10;
      iodigits++;
      ioradix = ioradix * 10;
   }

   if (iodigits <= 0) TerminalError("problem with I/O");
}

istream& operator>>(istream& s, ZZ& x)
{
   long c;
   long cval;
   long sign;
   long ndigits;
   long acc;
   NTL_ZZRegister(a);

   if (!s) NTL_INPUT_ERROR(s, "bad ZZ input");

   if (!iodigits) InitZZIO();

   a = 0;

   SkipWhiteSpace(s);
   c = s.peek();

   if (c == '-') {
      sign = -1;
      s.get();
      c = s.peek();
   }
   else
      sign = 1;

   cval = CharToIntVal(c);

   if (cval < 0 || cval > 9) NTL_INPUT_ERROR(s, "bad ZZ input");

   ndigits = 0;
   acc = 0;
   while (cval >= 0 && cval <= 9) {
      acc = acc*10 + cval;
      ndigits++;

      if (ndigits == iodigits) {
         mul(a, a, ioradix);
         add(a, a, acc);
         ndigits = 0;
         acc = 0;
      }

      s.get();
      c = s.peek();
      cval = CharToIntVal(c);
   }

   if (ndigits != 0) {
      long mpy = 1;
      while (ndigits > 0) {
         mpy = mpy * 10;
         ndigits--;
      }

      mul(a, a, mpy);
      add(a, a, acc);
   }

   if (sign == -1)
      negate(a, a);

   x = a;
   return s;
}

// The class _ZZ_local_stack should be defined in an empty namespace,
// but since I don't want to rely on namespaces, we just give it a funny
// name to avoid accidental name clashes.

struct _ZZ_local_stack {
   long top;
   Vec<long> data;

   _ZZ_local_stack() { top = -1; }

   long pop() { return data[top--]; }
   long empty() { return (top == -1); }
   void push(long x);
};

void _ZZ_local_stack::push(long x)
{
   if (top+1 >= data.length())
      data.SetLength(max(32, long(1.414*data.length())));

   top++;
   data[top] = x;
}

static
void PrintDigits(ostream& s, long d, long justify)
{
   NTL_THREAD_LOCAL static Vec<char> buf(INIT_SIZE, iodigits);

   long i = 0;

   while (d) {
      buf[i] = IntValToChar(d % 10);
      d = d / 10;
      i++;
   }

   if (justify) {
      long j = iodigits - i;
      while (j > 0) {
         s << "0";
         j--;
      }
   }

   while (i > 0) {
      i--;
      s << buf[i];
   }
}

ostream& operator<<(ostream& s, const ZZ& a)
{
   ZZ b;
   _ZZ_local_stack S;
   long r;
   long k;

   if (!iodigits) InitZZIO();

   b = a;

   k = sign(b);

   if (k == 0) {
      s << "0";
      return s;
   }

   if (k < 0) {
      s << "-";
      negate(b, b);
   }

   do {
      r = DivRem(b, b, ioradix);
      S.push(r);
   } while (!IsZero(b));

   r = S.pop();
   PrintDigits(s, r, 0);

   while (!S.empty()) {
      r = S.pop();
      PrintDigits(s, r, 1);
   }

   return s;
}

long GCD(long a, long b)
{
   long u, v, t, x;

   if (a < 0) {
      if (a < -NTL_MAX_LONG) ResourceError("GCD: integer overflow");
      a = -a;
   }

   if (b < 0) {
      if (b < -NTL_MAX_LONG) ResourceError("GCD: integer overflow");
      b = -b;
   }

   if (b==0)
      x = a;
   else {
      u = a;
      v = b;
      do {
         t = u % v;
         u = v;
         v = t;
      } while (v != 0);

      x = u;
   }

   return x;
}

void XGCD(long& d, long& s, long& t, long a, long b)
{
   long  u, v, u0, v0, u1, v1, u2, v2, q, r;

   long aneg = 0, bneg = 0;

   if (a < 0) {
      if (a < -NTL_MAX_LONG) ResourceError("XGCD: integer overflow");
      a = -a;
      aneg = 1;
   }

   if (b < 0) {
      if (b < -NTL_MAX_LONG) ResourceError("XGCD: integer overflow");
      b = -b;
      bneg = 1;
   }

   u1=1; v1=0;
   u2=0; v2=1;
   u = a; v = b;

   while (v != 0) {
      q = u / v;
      r = u % v;
      u = v;
      v = r;
      u0 = u2;
      v0 = v2;
      u2 =  u1 - q*u2;
      v2 = v1- q*v2;
      u1 = u0;
      v1 = v0;
   }

   if (aneg)
      u1 = -u1;

   if (bneg)
      v1 = -v1;

   d = u;
   s = u1;
   t = v1;
}

long InvMod(long a, long n)
{
   long d, s, t;

   XGCD(d, s, t, a, n);
   if (d != 1) InvModError("InvMod: inverse undefined");
   if (s < 0)
      return s + n;
   else
      return s;
}

long PowerMod(long a, long ee, long n)
{
   long x, y;

   unsigned long e;

   if (ee < 0)
      e = - ((unsigned long) ee);
   else
      e = ee;

   x = 1;
   y = a;
   while (e) {
      if (e & 1) x = MulMod(x, y, n);
      y = MulMod(y, y, n);
      e = e >> 1;
   }

   if (ee < 0) x = InvMod(x, n);

   return x;
}

long ProbPrime(long n, long NumTests)
{
   long m, x, y, z;
   long i, j, k;

   if (n <= 1) return 0;

   if (n == 2) return 1;
   if (n % 2 == 0) return 0;

   if (n == 3) return 1;
   if (n % 3 == 0) return 0;

   if (n == 5) return 1;
   if (n % 5 == 0) return 0;

   if (n == 7) return 1;
   if (n % 7 == 0) return 0;

   if (n >= NTL_SP_BOUND) {
      return ProbPrime(to_ZZ(n), NumTests);
   }

   m = n - 1;
   k = 0;
   while((m & 1) == 0) {
      m = m >> 1;
      k++;
   }

   // n - 1 == 2^k * m, m odd

   for (i = 0; i < NumTests; i++) {
      do {
         x = RandomBnd(n);
      } while (x == 0);
      // x == 0 is not a useful candidtae for a witness!

      if (x == 0) continue;
      z = PowerMod(x, m, n);
      if (z == 1) continue;

      j = 0;
      do {
         y = z;
         z = MulMod(y, y, n);
         j++;
      } while (j != k && z != 1);

      if (z != 1 || y !=  n-1) return 0;
   }

   return 1;
}

long MillerWitness(const ZZ& n, const ZZ& x)
{
   ZZ m, y, z;
   long j, k;

   if (x == 0) return 0;

   add(m, n, -1);
   k = MakeOdd(m);
   // n - 1 == 2^k * m, m odd

   PowerMod(z, x, m, n);
   if (z == 1) return 0;

   j = 0;
   do {
      y = z;
      SqrMod(z, y, n);
      j++;
   } while (j != k && z != 1);

   if (z != 1) return 1;
   add(y, y, 1);
   if (y != n) return 1;
   return 0;
}

// ComputePrimeBound computes a reasonable bound for trial
// division in the Miller-Rabin test.
// It is computed a bit on the "low" side, since being a bit
// low doesn't hurt much, but being too high can hurt a lot.

static
long ComputePrimeBound(long bn)
{
   long wn = (bn+NTL_ZZ_NBITS-1)/NTL_ZZ_NBITS;

   long fn;

   if (wn <= 36)
      fn = wn/4 + 1;
   else
      fn = long(1.67*sqrt(double(wn)));

   long prime_bnd;

   if (NumBits(bn) + NumBits(fn) > NTL_SP_NBITS)
      prime_bnd = NTL_SP_BOUND;
   else
      prime_bnd = bn*fn;

   return prime_bnd;
}

long ProbPrime(const ZZ& n, long NumTrials)
{
   if (n <= 1) return 0;

   if (n.SinglePrecision()) {
      return ProbPrime(to_long(n), NumTrials);
   }

   long prime_bnd = ComputePrimeBound(NumBits(n));

   PrimeSeq s;
   long p;

   p = s.next();
   while (p && p < prime_bnd) {
      if (rem(n, p) == 0)
         return 0;

      p = s.next();
   }

   ZZ W;
   W = 2;

   // first try W == 2....the exponentiation
   // algorithm runs slightly faster in this case

   if (MillerWitness(n, W))
      return 0;

   long i;

   for (i = 0; i < NumTrials; i++) {
      do {
         RandomBnd(W, n);
      } while (W == 0);
      // W == 0 is not a useful candidate for a witness!

      if (MillerWitness(n, W))
         return 0;
   }

   return 1;
}

void RandomPrime(ZZ& n, long l, long NumTrials)
{
   if (l <= 1)
      LogicError("RandomPrime: l out of range");

   if (l == 2) {
      if (RandomBnd(2))
         n = 3;
      else
         n = 2;

      return;
   }

   do {
      RandomLen(n, l);
      if (!IsOdd(n)) add(n, n, 1);
   } while (!ProbPrime(n, NumTrials));
}

void NextPrime(ZZ& n, const ZZ& m, long NumTrials)
{
   ZZ x;

   if (m <= 2) {
      n = 2;
      return;
   }

   x = m;

   while (!ProbPrime(x, NumTrials))
      add(x, x, 1);

   n = x;
}

long NextPrime(long m, long NumTrials)
{
   long x;

   if (m <= 2)
      return 2;

   x = m;

   while (x < NTL_SP_BOUND && !ProbPrime(x, NumTrials))
      x++;

   if (x >= NTL_SP_BOUND)
      ResourceError("NextPrime: no more primes");

   return x;
}

long NextPowerOfTwo(long m)
{
   long k;
   unsigned long n, um;

   if (m < 0) return 0;

   um = m;
   n = 1;
   k = 0;

   while (n < um) {
      n = n << 1;
      k++;
   }

   if (k >= NTL_BITS_PER_LONG-1)
      ResourceError("NextPowerOfTwo: overflow");

   return k;
}

long NumBits(long a)
{
   unsigned long aa;
   if (a < 0)
      aa = - ((unsigned long) a);
   else
      aa = a;

   long k = 0;
   while (aa) {
      k++;
      aa = aa >> 1;
   }

   return k;
}

long bit(long a, long k)
{
   unsigned long aa;
   if (a < 0)
      aa = - ((unsigned long) a);
   else
      aa = a;

   if (k < 0 || k >= NTL_BITS_PER_LONG)
      return 0;
   else
      return long((aa >> k) & 1);
}

long divide(ZZ& q, const ZZ& a, const ZZ& b)
{
   NTL_ZZRegister(qq);
   NTL_ZZRegister(r);

   if (IsZero(b)) {
      if (IsZero(a)) {
         clear(q);
         return 1;
      }
      else
         return 0;
   }

   if (IsOne(b)) {
      q = a;
      return 1;
   }

   DivRem(qq, r, a, b);
   if (!IsZero(r)) return 0;
   q = qq;
   return 1;
}

long divide(const ZZ& a, const ZZ& b)
{
   NTL_ZZRegister(r);

   if (IsZero(b)) return IsZero(a);
   if (IsOne(b)) return 1;

   rem(r, a, b);
   return IsZero(r);
}

long divide(ZZ& q, const ZZ& a, long b)
{
   NTL_ZZRegister(qq);

   if (!b) {
      if (IsZero(a)) {
         clear(q);
         return 1;
      }
      else
         return 0;
   }

   if (b == 1) {
      q = a;
      return 1;
   }

   long r = DivRem(qq, a, b);
   if (r) return 0;
   q = qq;
   return 1;
}

long divide(const ZZ& a, long b)
{
   if (!b) return IsZero(a);
   if (b == 1) {
      return 1;
   }

   long r = rem(a,  b);
   return (r == 0);
}

void InvMod(ZZ& x, const ZZ& a, const ZZ& n)
{
   // NOTE: the underlying LIP routines write to the first argument,
   // even if inverse is undefined

   NTL_ZZRegister(xx);
   if (InvModStatus(xx, a, n))
      InvModError("InvMod: inverse undefined", a, n);
   x = xx;
}

void PowerMod(ZZ& x, const ZZ& a, const ZZ& e, const ZZ& n)
{
   // NOTE: this ensures that all modular inverses are computed
   // in the routine InvMod above, rather than the LIP-internal
   // modular inverse routine
   if (e < 0) {
      ZZ a_inv;
      ZZ e_neg;

      InvMod(a_inv, a, n);
      negate(e_neg, e);
      LowLevelPowerMod(x, a_inv, e_neg, n);
   }
   else
      LowLevelPowerMod(x, a, e, n);
}

#ifdef NTL_EXCEPTIONS

void InvModError(const char *s, const ZZ& a, const ZZ& n)
{
   throw InvModErrorObject(s, a, n);
}

#else

void InvModError(const char *s, const ZZ& a, const ZZ& n)
{
   TerminalError(s);
}

#endif

long RandomPrime_long(long l, long NumTrials)
{
   if (l <= 1 || l >= NTL_BITS_PER_LONG)
      ResourceError("RandomPrime: length out of range");

   long n;
   do {
      n = RandomLen_long(l);
   } while (!ProbPrime(n, NumTrials));

   return n;
}

static Lazy< Vec<char> > lowsieve_storage;
// This is a GLOBAL VARIABLE

PrimeSeq::PrimeSeq()
{
   movesieve = 0;
   pshift = -1;
   pindex = -1;
   exhausted = 0;
}

long PrimeSeq::next()
{
   if (exhausted) {
      return 0;
   }

   if (pshift < 0) {
      shift(0);
      return 2;
   }

   for (;;) {
      const char *p = movesieve;
      long i = pindex;

      while ((++i) < NTL_PRIME_BND) {
         if (p[i]) {
            pindex = i;
            return pshift + 2 * i + 3;
         }
      }

      long newshift = pshift + 2*NTL_PRIME_BND;

      if (newshift > 2 * NTL_PRIME_BND * (2 * NTL_PRIME_BND + 1)) {
         /* end of the road */
         exhausted = 1;
         return 0;
      }

      shift(newshift);
   }
}

void PrimeSeq::shift(long newshift)
{
   long i;
   long j;
   long jstep;
   long jstart;
   long ibound;
   char *p;

   if (!lowsieve_storage.built())
      start();

   const char *lowsieve = lowsieve_storage->elts();

   if (newshift < 0) {
      pshift = -1;
   }
   else if (newshift == 0) {
      pshift = 0;
      movesieve = lowsieve;
   }
   else if (newshift != pshift) {
      if (movesieve_mem.length() == 0) {
         movesieve_mem.SetLength(NTL_PRIME_BND);
      }

      pshift = newshift;
      movesieve = p = movesieve_mem.elts();
      for (i = 0; i < NTL_PRIME_BND; i++)
         p[i] = 1;

      jstep = 3;
      ibound = pshift + 2 * NTL_PRIME_BND + 1;
      for (i = 0; jstep * jstep <= ibound; i++) {
         if (lowsieve[i]) {
            if (!((jstart = (pshift + 2) / jstep + 1) & 1))
               jstart++;
            if (jstart <= jstep)
               jstart = jstep;
            jstart = (jstart * jstep - pshift - 3) / 2;
            for (j = jstart; j < NTL_PRIME_BND; j += jstep)
               p[j] = 0;
         }
         jstep += 2;
      }
   }

   pindex = -1;
   exhausted = 0;
}

void PrimeSeq::start()
{
   long i;
   long j;
   long jstep;
   long jstart;
   long ibnd;
   char *p;

   do {
      Lazy< Vec<char> >::Builder builder(lowsieve_storage);
      if (!builder()) break;

      UniquePtr< Vec<char> > ptr;
      ptr.make();
      ptr->SetLength(NTL_PRIME_BND);

      p = ptr->elts();

      for (i = 0; i < NTL_PRIME_BND; i++)
         p[i] = 1;

      jstep = 1;
      jstart = -1;
      ibnd = (SqrRoot(2 * NTL_PRIME_BND + 1) - 3) / 2;
      for (i = 0; i <= ibnd; i++) {
         jstart += 2 * ((jstep += 2) - 1);
         if (p[i])
            for (j = jstart; j < NTL_PRIME_BND; j += jstep)
               p[j] = 0;
      }

      builder.move(ptr);
   } while (0);

}

void PrimeSeq::reset(long b)
{
   if (b > (2*NTL_PRIME_BND+1)*(2*NTL_PRIME_BND+1)) {
      exhausted = 1;
      return;
   }

   if (b <= 2) {
      shift(-1);
      return;
   }

   if ((b & 1) == 0) b++;

   shift(((b-3) / (2*NTL_PRIME_BND))* (2*NTL_PRIME_BND));
   pindex = (b - pshift - 3)/2 - 1;
}

long Jacobi(const ZZ& aa, const ZZ& nn)
{
   ZZ a, n;
   long t, k;
   long d;

   a = aa;
   n = nn;
   t = 1;

   while (a != 0) {
      k = MakeOdd(a);
      d = trunc_long(n, 3);
      if ((k & 1) && (d == 3 || d == 5)) t = -t;

      if (trunc_long(a, 2) == 3 && (d & 3) == 3) t = -t;
      swap(a, n);
      rem(a, a, n);
   }

   if (n == 1)
      return t;
   else
      return 0;
}

void SqrRootMod(ZZ& x, const ZZ& aa, const ZZ& nn)
{
   if (aa == 0 || aa == 1) {
      x = aa;
      return;
   }

   // at this point, we must have nn >= 5

   if (trunc_long(nn, 2) == 3) {  // special case, n = 3 (mod 4)
      ZZ n, a, e, z;

      n = nn;
      a  = aa;

      add(e, n, 1);
      RightShift(e, e, 2);

      PowerMod(z, a, e, n);
      x = z;

      return;
   }

   ZZ n, m;
   int h, nlen;

   n = nn;
   nlen = NumBits(n);

   sub(m, n, 1);
   h = MakeOdd(m);  // h >= 2

   if (nlen > 50 && h < SqrRoot(nlen)) {
      long i, j;
      ZZ a, b, a_inv, c, r, m1, d;

      a = aa;
      InvMod(a_inv, a, n);

      if (h == 2)
         b = 2;
      else {
         do {
            RandomBnd(b, n);
         } while (Jacobi(b, n) != -1);
      }

      PowerMod(c, b, m, n);

      add(m1, m, 1);
      RightShift(m1, m1, 1);
      PowerMod(r, a, m1, n);

      for (i = h-2; i >= 0; i--) {
         SqrMod(d, r, n);
         MulMod(d, d, a_inv, n);
         for (j = 0; j < i; j++)
            SqrMod(d, d, n);
         if (!IsOne(d))
            MulMod(r, r, c, n);
         SqrMod(c, c, n);
      }

      x = r;
      return;
   }

   long i, k;
   ZZ ma, t, u, v, e;
   ZZ t1, t2, t3, t4;

   n = nn;
   NegateMod(ma, aa, n);

   // find t such that t^2 - 4*a is not a square

   MulMod(t1, ma, 4, n);
   do {
      RandomBnd(t, n);
      SqrMod(t2, t, n);
      AddMod(t2, t2, t1, n);
   } while (Jacobi(t2, n) != -1);

   // compute u*X + v = X^{(n+1)/2} mod f, where f = X^2 - t*X + a

   add(e, n, 1);
   RightShift(e, e, 1);

   u = 0;
   v = 1;

   k = NumBits(e);

   for (i = k - 1; i >= 0; i--) {
      add(t2, u, v);
      sqr(t3, t2);  // t3 = (u+v)^2
      sqr(t1, u);
      sqr(t2, v);
      sub(t3, t3, t1);
      sub(t3, t3, t2); // t1 = u^2, t2 = v^2, t3 = 2*u*v
      rem(t1, t1, n);
      mul(t4, t1, t);
      add(t4, t4, t3);
      rem(u, t4, n);

      mul(t4, t1, ma);
      add(t4, t4, t2);
      rem(v, t4, n);

      if (bit(e, i)) {
         MulMod(t1, u, t, n);
         AddMod(t1, t1, v, n);
         MulMod(v, u, ma, n);
         u = t1;
      }

   }

   x = v;
}

// Chinese Remaindering.
//
// This version in new to v3.7, and is significantly
// simpler and faster than the previous version.
//
// This function takes as input g, a, G, p,
// such that a > 0, 0 <= G < p, and gcd(a, p) = 1.
// It computes a' = a*p and g' such that
//   * g' = g (mod a);
//   * g' = G (mod p);
//   * -a'/2 < g' <= a'/2.
// It then sets g := g' and a := a', and returns 1 iff g has changed.
//
// Under normal use, the input value g satisfies -a/2 < g <= a/2;
// however, this was not documented or enforced in earlier versions,
// so to maintain backward compatability, no restrictions are placed
// on g.  This routine runs faster, though, if -a/2 < g <= a/2,
// and the first thing the routine does is to make this condition
// hold.
//
// Also, under normal use, both a and p are odd;  however, the routine
// will still work even if this is not so.
//
// The routine is based on the following simple fact.
//
// Let -a/2 < g <= a/2, and let h satisfy
//   * g + a h = G (mod p);
//   * -p/2 < h <= p/2.
// Further, if p = 2*h and g > 0, set
//   g' := g - a h;
// otherwise, set
//   g' := g + a h.
// Then g' so defined satisfies the above requirements.
//
// It is trivial to see that g's satisfies the congruence conditions.
// The only thing is to check that the "balancing" condition
// -a'/2 < g' <= a'/2 also holds.

long CRT(ZZ& gg, ZZ& a, long G, long p)
{
   if (p >= NTL_SP_BOUND) {
      ZZ GG, pp;
      conv(GG, G);
      conv(pp, p);
      return CRT(gg, a, GG, pp);
   }

   long modified = 0;

   NTL_ZZRegister(g);

   if (!CRTInRange(gg, a)) {
      modified = 1;
      ZZ a1;
      rem(g, gg, a);
      RightShift(a1, a, 1);
      if (g > a1) sub(g, g, a);
   }
   else
      g = gg;

   long p1;
   p1 = p >> 1;

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

   long h;
   h = rem(g, p);
   h = SubMod(G, h, p);
   h = MulMod(h, a_inv, p);
   if (h > p1)
      h = h - p;

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

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

   mul(a, a, p);
   gg = g;

   return modified;
}

long CRT(ZZ& gg, ZZ& a, const ZZ& G, const ZZ& p)
{
   long modified = 0;

   ZZ g;

   if (!CRTInRange(gg, a)) {
      modified = 1;
      ZZ a1;
      rem(g, gg, a);
      RightShift(a1, a, 1);
      if (g > a1) sub(g, g, a);
   }
   else
      g = gg;

   ZZ p1;
   RightShift(p1, p, 1);

   ZZ a_inv;
   rem(a_inv, a, p);
   InvMod(a_inv, a_inv, p);

   ZZ h;
   rem(h, g, p);
   SubMod(h, G, h, p);
   MulMod(h, h, a_inv, p);
   if (h > p1)
      sub(h, h, p);

   if (h != 0) {
      modified = 1;
      ZZ ah;
      mul(ah, a, h);

      if (!IsOdd(p) && g > 0 &&  (h == p1))
         sub(g, g, ah);
      else
         add(g, g, ah);
   }

   mul(a, a, p);
   gg = g;

   return modified;
}

void sub(ZZ& x, const ZZ& a, long b)
{
   NTL_ZZRegister(B);
   conv(B, b);
   sub(x, a, B);
}

void sub(ZZ& x, long a, const ZZ& b)
{
   NTL_ZZRegister(A);
   conv(A, a);
   sub(x, A, b);
}

void power2(ZZ& x, long e)
{
   if (e < 0) ArithmeticError("power2: negative exponent");
   set(x);
   LeftShift(x, x, e);
}

void conv(ZZ& x, const char *s)
{
   long c;
   long cval;
   long sign;
   long ndigits;
   long acc;
   long i = 0;

   NTL_ZZRegister(a);

   if (!s) InputError("bad ZZ input");

   if (!iodigits) InitZZIO();

   a = 0;

   c = s[i];
   while (IsWhiteSpace(c)) {
      i++;
      c = s[i];
   }

   if (c == '-') {
      sign = -1;
      i++;
      c = s[i];
   }
   else
      sign = 1;

   cval = CharToIntVal(c);
   if (cval < 0 || cval > 9) InputError("bad ZZ input");

   ndigits = 0;
   acc = 0;
   while (cval >= 0 && cval <= 9) {
      acc = acc*10 + cval;
      ndigits++;

      if (ndigits == iodigits) {
         mul(a, a, ioradix);
         add(a, a, acc);
         ndigits = 0;
         acc = 0;
      }

      i++;
      c = s[i];
      cval = CharToIntVal(c);
   }

   if (ndigits != 0) {
      long mpy = 1;
      while (ndigits > 0) {
         mpy = mpy * 10;
         ndigits--;
      }

      mul(a, a, mpy);
      add(a, a, acc);
   }

   if (sign == -1)
      negate(a, a);

   x = a;
}

void bit_and(ZZ& x, const ZZ& a, long b)
{
   NTL_ZZRegister(B);
   conv(B, b);
   bit_and(x, a, B);
}

void bit_or(ZZ& x, const ZZ& a, long b)
{
   NTL_ZZRegister(B);
   conv(B, b);
   bit_or(x, a, B);
}

void bit_xor(ZZ& x, const ZZ& a, long b)
{
   NTL_ZZRegister(B);
   conv(B, b);
   bit_xor(x, a, B);
}

long power_long(long a, long e)
{
   if (e < 0) ArithmeticError("power_long: negative exponent");

   if (e == 0) return 1;

   if (a == 1) return 1;
   if (a == -1) {
      if (e & 1)
         return -1;
      else
         return 1;
   }

   // no overflow check --- result is computed correctly
   // modulo word size

   unsigned long res = 1;
   unsigned long aa = a;
   long i;

   for (i = 0; i < e; i++)
      res *= aa;

   return to_long(res);
}

//  RANDOM NUMBER GENERATION

// Idea for this PRNG.  Iteratively hash seed using md5
// to get 256 bytes to initialize arc4.
// Then use arc4 to get a pseudo-random byte stream.

// I've taken care that the pseudo-random numbers generated by
// the routines RandomBnd, RandomBits, and RandomLen
// are completely platform independent.

// I make use of the md5 compression function,
// which I've modified to work on 64-bit machines

/*
 *  BEGIN RSA's md5 stuff
 *
 */

/*
 **********************************************************************
 ** md5.c                                                            **
 ** RSA Data Security, Inc. MD5 Message Digest Algorithm             **
 ** Created: 2/17/90 RLR                                             **
 ** Revised: 1/91 SRD,AJ,BSK,JT Reference C Version                  **
 **********************************************************************
 */

/*
 **********************************************************************
 ** Copyright (C) 1990, RSA Data Security, Inc. All rights reserved. **
 **                                                                  **
 ** License to copy and use this software is granted provided that   **
 ** it is identified as the "RSA Data Security, Inc. MD5 Message     **
 ** Digest Algorithm" in all material mentioning or referencing this **
 ** software or this function.                                       **
 **                                                                  **
 ** License is also granted to make and use derivative works         **
 ** provided that such works are identified as "derived from the RSA **
 ** Data Security, Inc. MD5 Message Digest Algorithm" in all         **
 ** material mentioning or referencing the derived work.             **
 **                                                                  **
 ** RSA Data Security, Inc. makes no representations concerning      **
 ** either the merchantability of this software or the suitability   **
 ** of this software for any particular purpose.  It is provided "as **
 ** is" without express or implied warranty of any kind.             **
 **                                                                  **
 ** These notices must be retained in any copies of any part of this **
 ** documentation and/or software.                                   **
 **********************************************************************
 */

#if (NTL_BITS_PER_LONG <= 32)
#define TRUNC32(x) (x)
#else
#define TRUNC32(x) ((x) & ((1UL << 32)-1UL))
#endif

/* F, G and H are basic MD5 functions: selection, majority, parity */
#define F(x, y, z) (((x) & (y)) | ((~x) & (z)))
#define G(x, y, z) (((x) & (z)) | ((y) & (~z)))
#define H(x, y, z) ((x) ^ (y) ^ (z))
#define I(x, y, z) (TRUNC32((y) ^ ((x) | (~z))))

/* ROTATE_LEFT rotates x left n bits */
#define ROTATE_LEFT(x, n) (TRUNC32(((x) << (n)) | ((x) >> (32-(n)))))

/* FF, GG, HH, and II transformations for rounds 1, 2, 3, and 4 */
/* Rotation is separate from addition to prevent recomputation */
#define FF(a, b, c, d, x, s, ac) \
  {(a) = TRUNC32((a) + F((b), (c), (d)) + (x) + (ac)); \
   (a) = ROTATE_LEFT((a), (s)); \
   (a) = TRUNC32((a) + (b)); \
  }
#define GG(a, b, c, d, x, s, ac) \
  {(a) = TRUNC32((a) + G((b), (c), (d)) + (x) + (ac)); \
   (a) = ROTATE_LEFT((a), (s)); \
   (a) = TRUNC32((a) + (b)); \
  }
#define HH(a, b, c, d, x, s, ac) \
  {(a) = TRUNC32((a) + H((b), (c), (d)) + (x) + (ac)); \
   (a) = ROTATE_LEFT((a), (s)); \
   (a) = TRUNC32((a) + (b)); \
  }
#define II(a, b, c, d, x, s, ac) \
  {(a) = TRUNC32((a) + I((b), (c), (d)) + (x) + (ac)); \
   (a) = ROTATE_LEFT((a), (s)); \
   (a) = TRUNC32((a) + (b)); \
  }

static
void MD5_default_IV(unsigned long *buf)
{
   buf[0] = 0x67452301UL;
   buf[1] = 0xefcdab89UL;
   buf[2] = 0x98badcfeUL;
   buf[3] = 0x10325476UL;
}

/* Basic MD5 step. Transform buf based on in.
 */

static
void MD5_compress(unsigned long *buf, unsigned long *in)
{
  unsigned long a = buf[0], b = buf[1], c = buf[2], d = buf[3];

  /* Round 1 */
#define S11 7
#define S12 12
#define S13 17
#define S14 22
  FF ( a, b, c, d, in[ 0], S11, 3614090360UL); /* 1 */
  FF ( d, a, b, c, in[ 1], S12, 3905402710UL); /* 2 */
  FF ( c, d, a, b, in[ 2], S13,  606105819UL); /* 3 */
  FF ( b, c, d, a, in[ 3], S14, 3250441966UL); /* 4 */
  FF ( a, b, c, d, in[ 4], S11, 4118548399UL); /* 5 */
  FF ( d, a, b, c, in[ 5], S12, 1200080426UL); /* 6 */
  FF ( c, d, a, b, in[ 6], S13, 2821735955UL); /* 7 */
  FF ( b, c, d, a, in[ 7], S14, 4249261313UL); /* 8 */
  FF ( a, b, c, d, in[ 8], S11, 1770035416UL); /* 9 */
  FF ( d, a, b, c, in[ 9], S12, 2336552879UL); /* 10 */
  FF ( c, d, a, b, in[10], S13, 4294925233UL); /* 11 */
  FF ( b, c, d, a, in[11], S14, 2304563134UL); /* 12 */
  FF ( a, b, c, d, in[12], S11, 1804603682UL); /* 13 */
  FF ( d, a, b, c, in[13], S12, 4254626195UL); /* 14 */
  FF ( c, d, a, b, in[14], S13, 2792965006UL); /* 15 */
  FF ( b, c, d, a, in[15], S14, 1236535329UL); /* 16 */

  /* Round 2 */
#define S21 5
#define S22 9
#define S23 14
#define S24 20
  GG ( a, b, c, d, in[ 1], S21, 4129170786UL); /* 17 */
  GG ( d, a, b, c, in[ 6], S22, 3225465664UL); /* 18 */
  GG ( c, d, a, b, in[11], S23,  643717713UL); /* 19 */
  GG ( b, c, d, a, in[ 0], S24, 3921069994UL); /* 20 */
  GG ( a, b, c, d, in[ 5], S21, 3593408605UL); /* 21 */
  GG ( d, a, b, c, in[10], S22,   38016083UL); /* 22 */
  GG ( c, d, a, b, in[15], S23, 3634488961UL); /* 23 */
  GG ( b, c, d, a, in[ 4], S24, 3889429448UL); /* 24 */
  GG ( a, b, c, d, in[ 9], S21,  568446438UL); /* 25 */
  GG ( d, a, b, c, in[14], S22, 3275163606UL); /* 26 */
  GG ( c, d, a, b, in[ 3], S23, 4107603335UL); /* 27 */
  GG ( b, c, d, a, in[ 8], S24, 1163531501UL); /* 28 */
  GG ( a, b, c, d, in[13], S21, 2850285829UL); /* 29 */
  GG ( d, a, b, c, in[ 2], S22, 4243563512UL); /* 30 */
  GG ( c, d, a, b, in[ 7], S23, 1735328473UL); /* 31 */
  GG ( b, c, d, a, in[12], S24, 2368359562UL); /* 32 */

  /* Round 3 */
#define S31 4
#define S32 11
#define S33 16
#define S34 23
  HH ( a, b, c, d, in[ 5], S31, 4294588738UL); /* 33 */
  HH ( d, a, b, c, in[ 8], S32, 2272392833UL); /* 34 */
  HH ( c, d, a, b, in[11], S33, 1839030562UL); /* 35 */
  HH ( b, c, d, a, in[14], S34, 4259657740UL); /* 36 */
  HH ( a, b, c, d, in[ 1], S31, 2763975236UL); /* 37 */
  HH ( d, a, b, c, in[ 4], S32, 1272893353UL); /* 38 */
  HH ( c, d, a, b, in[ 7], S33, 4139469664UL); /* 39 */
  HH ( b, c, d, a, in[10], S34, 3200236656UL); /* 40 */
  HH ( a, b, c, d, in[13], S31,  681279174UL); /* 41 */
  HH ( d, a, b, c, in[ 0], S32, 3936430074UL); /* 42 */
  HH ( c, d, a, b, in[ 3], S33, 3572445317UL); /* 43 */
  HH ( b, c, d, a, in[ 6], S34,   76029189UL); /* 44 */
  HH ( a, b, c, d, in[ 9], S31, 3654602809UL); /* 45 */
  HH ( d, a, b, c, in[12], S32, 3873151461UL); /* 46 */
  HH ( c, d, a, b, in[15], S33,  530742520UL); /* 47 */
  HH ( b, c, d, a, in[ 2], S34, 3299628645UL); /* 48 */

  /* Round 4 */
#define S41 6
#define S42 10
#define S43 15
#define S44 21
  II ( a, b, c, d, in[ 0], S41, 4096336452UL); /* 49 */
  II ( d, a, b, c, in[ 7], S42, 1126891415UL); /* 50 */
  II ( c, d, a, b, in[14], S43, 2878612391UL); /* 51 */
  II ( b, c, d, a, in[ 5], S44, 4237533241UL); /* 52 */
  II ( a, b, c, d, in[12], S41, 1700485571UL); /* 53 */
  II ( d, a, b, c, in[ 3], S42, 2399980690UL); /* 54 */
  II ( c, d, a, b, in[10], S43, 4293915773UL); /* 55 */
  II ( b, c, d, a, in[ 1], S44, 2240044497UL); /* 56 */
  II ( a, b, c, d, in[ 8], S41, 1873313359UL); /* 57 */
  II ( d, a, b, c, in[15], S42, 4264355552UL); /* 58 */
  II ( c, d, a, b, in[ 6], S43, 2734768916UL); /* 59 */
  II ( b, c, d, a, in[13], S44, 1309151649UL); /* 60 */
  II ( a, b, c, d, in[ 4], S41, 4149444226UL); /* 61 */
  II ( d, a, b, c, in[11], S42, 3174756917UL); /* 62 */
  II ( c, d, a, b, in[ 2], S43,  718787259UL); /* 63 */
  II ( b, c, d, a, in[ 9], S44, 3951481745UL); /* 64 */

  buf[0] = TRUNC32(buf[0] + a);
  buf[1] = TRUNC32(buf[1] + b);
  buf[2] = TRUNC32(buf[2] + c);
  buf[3] = TRUNC32(buf[3] + d);
}

/*
 *  END RSA's md5 stuff
 *
 */

static
void words_from_bytes(unsigned long *txtl, const unsigned char *txtc, long n)
{
   long i;
   unsigned long v;

   for (i = 0; i < n; i++) {
      v = txtc[4*i];
      v += ((unsigned long) (txtc[4*i+1])) << 8;
      v += ((unsigned long) (txtc[4*i+2])) << 16;
      v += ((unsigned long) (txtc[4*i+3])) << 24;
      txtl[i] = v;
   }
}

static
void bytes_from_words(unsigned char *txtc, const unsigned long *txtl, long n)
{
   long i;
   unsigned long v;

   for (i = 0; i < n; i++) {
      v = txtl[i];
      txtc[4*i] = v & 255;
      v = v >> 8;
      txtc[4*i+1] = v & 255;
      v = v >> 8;
      txtc[4*i+2] = v & 255;
      v = v >> 8;
      txtc[4*i+3] = v & 255;
   }
}

static
void MD5_compress1(unsigned long *buf, unsigned char *in, long n)
{
   unsigned long txtl[16];
   unsigned char txtc[64];
   long i, j, k;

   if (n < 0) n = 0;

   i = 0;
   while (i < n) {
      k = n-i;
      if (k > 64) k = 64;
      for (j = 0; j < k; j++)
         txtc[j] = in[i+j];
      for (; j < 64; j++)
         txtc[j] = 0;
      words_from_bytes(txtl, txtc, 16);
      MD5_compress(buf, txtl);
      i += k;
   }
}

// the "cipherpunk" version of arc4

struct _ZZ_arc4_key
{
    unsigned char state[256];
    unsigned char x;
    unsigned char y;
};

static inline
void swap_byte(unsigned char *a, unsigned char *b)
{
    unsigned char swapByte;

    swapByte = *a;
    *a = *b;
    *b = swapByte;
}

static
void prepare_key(unsigned char *key_data_ptr,
                 long key_data_len, _ZZ_arc4_key *key)
{
    unsigned char index1;
    unsigned char index2;
    unsigned char* state;
    long counter;

    state = &key->state[0];
    for(counter = 0; counter < 256; counter++)
       state[counter] = counter;
    key->x = 0;
    key->y = 0;
    index1 = 0;
    index2 = 0;
    for(counter = 0; counter < 256; counter++)
    {
         index2 = (key_data_ptr[index1] + state[counter] + index2) & 255;
         swap_byte(&state[counter], &state[index2]);

         index1 = (index1 + 1) % key_data_len;
    }
}

static
void arc4(unsigned char *buffer_ptr, long buffer_len, _ZZ_arc4_key *key)
{
    unsigned char x;
    unsigned char y;
    unsigned char* state;
    unsigned char xorIndex;
    long counter;

    x = key->x;
    y = key->y;

    state = &key->state[0];
    for(counter = 0; counter < buffer_len; counter ++)
    {
         x = (x + 1) & 255;
         y = (state[x] + y) & 255;
         swap_byte(&state[x], &state[y]);

         xorIndex = (state[x] + state[y]) & 255;

         buffer_ptr[counter] = state[xorIndex];
     }
     key->x = x;
     key->y = y;
}

// global state information for PRNG

NTL_THREAD_LOCAL static long ran_initialized = 0;
NTL_THREAD_LOCAL static _ZZ_arc4_key ran_key;

static const unsigned long default_md5_tab[16] = {
744663023UL, 1011602954UL, 3163087192UL, 3383838527UL,
3305324122UL, 3197458079UL, 2266495600UL, 2760303563UL,
346234297UL, 1919920720UL, 1896169861UL, 2192176675UL,
2027150322UL, 2090160759UL, 2134858730UL, 1131796244UL
};

static
void build_arc4_tab(unsigned char *seed_bytes, const ZZ& s)
{
   long nb = NumBytes(s);

   unsigned char *txt;

   Vec<unsigned char> txt_storage;
   txt_storage.SetLength(nb + 68);
   txt = txt_storage.elts();

   BytesFromZZ(txt + 4, s, nb);

   bytes_from_words(txt + nb + 4, default_md5_tab, 16);

   unsigned long buf[4];

   unsigned long i;
   for (i = 0; i < 16; i++) {
      MD5_default_IV(buf);
      bytes_from_words(txt, &i, 1);

      MD5_compress1(buf, txt, nb + 68);

      bytes_from_words(seed_bytes + 16*i, buf, 4);
   }
}

void SetSeed(const ZZ& s)
{
   unsigned char seed_bytes[256];

   build_arc4_tab(seed_bytes, s);
   prepare_key(seed_bytes, 256, &ran_key);

   // we discard the first 1024 bytes of the arc4 stream, as this is
   // recommended practice.

   arc4(seed_bytes, 256, &ran_key);
   arc4(seed_bytes, 256, &ran_key);
   arc4(seed_bytes, 256, &ran_key);
   arc4(seed_bytes, 256, &ran_key);

   ran_initialized = 1;
}

static
void ran_bytes(unsigned char *bytes, long n)
{
   if (!ran_initialized) {
      ZZ x;
      const string& id = UniqueID();
      ZZFromBytes(x, (const unsigned char *) id.c_str(), id.length());
      SetSeed(x);
   }
   arc4(bytes, n, &ran_key);
}

unsigned long RandomWord()
{
   unsigned char buf[NTL_BITS_PER_LONG/8];
   long i;
   unsigned long res;

   ran_bytes(buf, NTL_BITS_PER_LONG/8);

   res = 0;
   for (i = NTL_BITS_PER_LONG/8 - 1; i >= 0; i--) {
      res = res << 8;
      res = res | buf[i];
   }

   return res;
}

long RandomBits_long(long l)
{
   if (l <= 0) return 0;
   if (l >= NTL_BITS_PER_LONG)
      ResourceError("RandomBits: length too big");

   unsigned char buf[NTL_BITS_PER_LONG/8];
   unsigned long res;
   long i;

   long nb = (l+7)/8;
   ran_bytes(buf, nb);

   res = 0;
   for (i = nb - 1; i >= 0; i--) {
      res = res << 8;
      res = res | buf[i];
   }

   return long(res & ((1UL << l)-1UL));
}

unsigned long RandomBits_ulong(long l)
{
   if (l <= 0) return 0;
   if (l > NTL_BITS_PER_LONG)
      ResourceError("RandomBits: length too big");

   unsigned char buf[NTL_BITS_PER_LONG/8];
   unsigned long res;
   long i;

   long nb = (l+7)/8;
   ran_bytes(buf, nb);

   res = 0;
   for (i = nb - 1; i >= 0; i--) {
      res = res << 8;
      res = res | buf[i];
   }

   if (l < NTL_BITS_PER_LONG)
      res = res & ((1UL << l)-1UL);

   return res;
}

long RandomLen_long(long l)
{
   if (l <= 0) return 0;
   if (l == 1) return 1;
   if (l >= NTL_BITS_PER_LONG)
      ResourceError("RandomLen: length too big");

   return RandomBits_long(l-1) + (1L << (l-1));
}

void RandomBits(ZZ& x, long l)
{
   if (l <= 0) {
      x = 0;
      return;
   }

   if (NTL_OVERFLOW(l, 1, 0))
      ResourceError("RandomBits: length too big");

   long nb = (l+7)/8;

   NTL_THREAD_LOCAL static Vec<unsigned char> buf_mem;
   Vec<unsigned char>::Watcher watch_buf_mem(buf_mem);

   buf_mem.SetLength(nb);
   unsigned char *buf = buf_mem.elts();

   ran_bytes(buf, nb);

   NTL_ZZRegister(res);

   ZZFromBytes(res, buf, nb);
   trunc(res, res, l);

   x = res;
}

void RandomLen(ZZ& x, long l)
{
   if (l <= 0) {
      x = 0;
      return;
   }

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

   if (NTL_OVERFLOW(l, 1, 0))
      ResourceError("RandomLen: length too big");

   // pre-allocate space to avoid two allocations
   long nw = (l + NTL_ZZ_NBITS - 1)/NTL_ZZ_NBITS;
   x.SetSize(nw);

   RandomBits(x, l-1);
   SetBit(x, l-1);
}

const long RandomBndExcess = 8;

void RandomBnd(ZZ& x, const ZZ& bnd)
{
   if (bnd <= 1) {
      x = 0;
      return;
   }

   long k = NumBits(bnd);

   if (weight(bnd) == 1) {
      RandomBits(x, k-1);
      return;
   }

   long l = k + RandomBndExcess;

   NTL_ZZRegister(t);
   NTL_ZZRegister(r);
   NTL_ZZRegister(t1);

   do {
      RandomBits(t, l);
      rem(r, t, bnd);
      sub(t1, bnd, r);
      add(t, t, t1);
   } while (NumBits(t) > l);

   x = r;
}

long RandomBnd(long bnd)
{
   if (bnd <= 1) return 0;

   long k = NumBits(bnd);

   if (((bnd - 1) & bnd) == 0)
      return RandomBits_long(k-1);

   long l = k + RandomBndExcess;

   if (l > NTL_BITS_PER_LONG-2) {
      NTL_ZZRegister(Bnd);
      NTL_ZZRegister(res);

      Bnd = bnd;
      RandomBnd(res, Bnd);
      return to_long(res);
   }

   long t, r;

   do {
      t = RandomBits_long(l);
      r = t % bnd;
   } while (t + bnd - r > (1L << l));

   return r;
}

// More prime generation stuff...

static
double Log2(double x)
{
   NTL_THREAD_LOCAL static double log2 = log(2.0);
   return log(x)/log2;
}

// Define p(k,t) to be the conditional probability that a random, odd, k-bit
// number is composite, given that it passes t iterations of the
// Miller-Rabin test.
// This routine returns 0 or 1, and if it returns 1 then
// p(k,t) <= 2^{-n}.
// This basically encodes the estimates of Damgard, Landrock, and Pomerance;
// it uses floating point arithmetic, but is coded in such a way
// that its results should be correct, assuming that the log function
// is computed with reasonable precision.
//
// It is assumed that k >= 3 and t >= 1; if this does not hold,
// then 0 is returned.

static
long ErrBoundTest(long kk, long tt, long nn)

{
   const double fudge = (1.0 + 1024.0/NTL_FDOUBLE_PRECISION);
   const double log2_3 = Log2(3.0);
   const double log2_7 = Log2(7.0);
   const double log2_20 = Log2(20.0);

   double k = kk;
   double t = tt;
   double n = nn;

   if (k < 3 || t < 1) return 0;
   if (n < 1) return 1;

   // the following test is largely academic
   if (9*t > NTL_FDOUBLE_PRECISION) LogicError("ErrBoundTest: t too big");

   double log2_k = Log2(k);

   if ((n + log2_k)*fudge <= 2*t)
      return 1;

   if ((2*log2_k + 4.0 + n)*fudge <= 2*sqrt(k))
      return 2;

   if ((t == 2 && k >= 88) || (3 <= t && 9*t <= k && k >= 21)) {
      if ((1.5*log2_k + t + 4.0 + n)*fudge <= 0.5*Log2(t) + 2*(sqrt(t*k)))
         return 3;
   }

   if (k <= 9*t && 4*t <= k && k >= 21) {
      if ( ((log2_3 + log2_7 + log2_k + n)*fudge <= log2_20 + 5*t)  &&
           ((log2_3 + (15.0/4.0)*log2_k + n)*fudge <= log2_7 + k/2 + 2*t) &&
           ((2*log2_3 + 2 + log2_k + n)*fudge <= k/4 + 3*t) )
         return 4;
   }

   if (4*t >= k && k >= 21) {
      if (((15.0/4.0)*log2_k + n)*fudge <= log2_7 + k/2 + 2*t)
         return 5;
   }

   return 0;
}

void GenPrime(ZZ& n, long k, long err)
{
   if (k <= 1) LogicError("GenPrime: bad length");

   if (k > (1L << 20)) ResourceError("GenPrime: length too large");

   if (err < 1) err = 1;
   if (err > 512) err = 512;

   if (k == 2) {
      if (RandomBnd(2))
         n = 3;
      else
         n = 2;

      return;
   }

   long t;

   t = 1;
   while (!ErrBoundTest(k, t, err))
      t++;

   RandomPrime(n, k, t);
}

long GenPrime_long(long k, long err)
{
   if (k <= 1) LogicError("GenPrime: bad length");

   if (k >= NTL_BITS_PER_LONG) ResourceError("GenPrime: length too large");

   if (err < 1) err = 1;
   if (err > 512) err = 512;

   if (k == 2) {
      if (RandomBnd(2))
         return 3;
      else
         return 2;
   }

   long t;

   t = 1;
   while (!ErrBoundTest(k, t, err))
      t++;

   return RandomPrime_long(k, t);
}

void GenGermainPrime(ZZ& n, long k, long err)
{
   if (k <= 1) LogicError("GenGermainPrime: bad length");

   if (k > (1L << 20)) ResourceError("GenGermainPrime: length too large");

   if (err < 1) err = 1;
   if (err > 512) err = 512;

   if (k == 2) {
      if (RandomBnd(2))
         n = 3;
      else
         n = 2;

      return;
   }

   long prime_bnd = ComputePrimeBound(k);

   if (NumBits(prime_bnd) >= k/2)
      prime_bnd = (1L << (k/2-1));

   ZZ two;
   two = 2;

   ZZ n1;

   PrimeSeq s;

   ZZ iter;
   iter = 0;

   for (;;) {
      iter++;

      RandomLen(n, k);
      if (!IsOdd(n)) add(n, n, 1);

      s.reset(3);
      long p;

      long sieve_passed = 1;

      p = s.next();
      while (p && p < prime_bnd) {
         long r = rem(n, p);

         if (r == 0) {
            sieve_passed = 0;
            break;
         }

         // test if 2*r + 1 = 0 (mod p)
         if (r == p-r-1) {
            sieve_passed = 0;
            break;
         }

         p = s.next();
      }

      if (!sieve_passed) continue;

      if (MillerWitness(n, two)) continue;

      // n1 = 2*n+1
      mul(n1, n, 2);
      add(n1, n1, 1);

      if (MillerWitness(n1, two)) continue;

      // now do t M-R iterations...just to make sure

      // First compute the appropriate number of M-R iterations, t
      // The following computes t such that
      //       p(k,t)*8/k <= 2^{-err}/(5*iter^{1.25})
      // which suffices to get an overall error probability of 2^{-err}.
      // Note that this method has the advantage of not requiring
      // any assumptions on the density of Germain primes.

      long err1 = max(1, err + 7 + (5*NumBits(iter) + 3)/4 - NumBits(k));
      long t;
      t = 1;
      while (!ErrBoundTest(k, t, err1))
         t++;

      ZZ W;
      long MR_passed = 1;

      long i;
      for (i = 1; i <= t; i++) {
         do {
            RandomBnd(W, n);
         } while (W == 0);
         // W == 0 is not a useful candidate witness!

         if (MillerWitness(n, W)) {
            MR_passed = 0;
            break;
         }
      }

      if (MR_passed) break;
   }
}

long GenGermainPrime_long(long k, long err)
{
   if (k >= NTL_BITS_PER_LONG-1)
      ResourceError("GenGermainPrime_long: length too long");

   ZZ n;
   GenGermainPrime(n, k, err);
   return to_long(n);
}

NTL_END_IMPL

max_stars_repo_name: Brainloop-Security/secret-sharing

max_issues_repo_name: Brainloop-Security/secret-sharing

max_forks_repo_name: Brainloop-Security/secret-sharing


text:
/* Copyright (C) 2017 IBM Corp.
 *  Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *     http://www.apache.org/licenses/LICENSE-2.0
 * Unless required by applicable law or agreed to in writing,
 * software distributed under the License is distributed on an
 * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,
 * either express or implied. See the License for the specific
 * language governing permissions and limitations under the License.
 */
/*****************************************************************************
 * TDMatrix.cpp - Matrixes with trapdoors,
 *   implementing the Efficient algorithm SampleD - Algorithm 3 from MP12.
 *****************************************************************************/
#include <stdexcept>
#include <ctime>
#include <NTL/mat_ZZ.h>
#include <NTL/mat_lzz_p.h>
#include <iostream>
#include <fstream>

NTL_CLIENT
#include "utils/timing.h"
#include "TDMatrix.h"
#include "mat_l.h"

//#define DEBUGPRINT
//#define DEBUG
//#define PRINTDOT //show how program progresses

// A static variable for debugging purposes
std::atomic<int> TDMatrix::maxSample(0); // keep the largest sample ever drawn

/***** Some local functions *****/

// Compute the covariance matrix sigmaP = sigmaX*I - sigmaG*[R/I]*(R^t|I)
static void
computeCovarianceMatrix(mat_l& covMat,const mat_l& R,long sigmaG,long sigmaX);

// Multiply a vector by the G matrix modulo the current NTL modulus, y=G*x.
// G is and n-by-m matrix with m=n*e*numOfFactors, and is defind by means
// of the vector of factors and the exponent e, each row of G is mostly zero,
// except for a progression
//    ( 1, f1,...,f1^e,  (f1^e f2), (f1^e f2^2), ..., (f1^e...fk^{e-1}) )
template<class T>
void multByG(vec_zz_p &y, const Vec<T> &x,
             long n, const vec_l& factors, long e)
{
    FHE_TIMER_START;
    long k = factors.length();
    long m = n*k*e;

    assert(x.length()==m);
    y.SetLength(n);

    for (long row=0; row<n; row++)   // Multiply x by next row of G
    {
        long index = row*k*e;      // first non-zero entry in this row
        zz_p val = to_zz_p(1L);  // value in current entry of G
        y[row] = to_zz_p(0L);    // accumulated sum
        for (long f=0; f<k; f++)     // go over all the factors
        {
            zz_p fact = to_zz_p(factors[f]); // convert to zz_p
            for (long i=0; i<e; i++)    // use each factor e times
            {
                zz_p addedTerm = val;
                addedTerm *= x[index++];  // update sum, advance index
                y[row] += addedTerm;
                val *= fact;              // compute next entry in this row of G
            }
        } // end of current factor
    } // end of current row
}
// Explicit instantiations for Vec<long> and Vec<zz_p>
template void multByG<zz_p>(vec_zz_p &y, const Vec<zz_p> &x,
                            long n, const vec_l& factors, long e);
template void multByG<long>(vec_zz_p &y, const Vec<long> &x,
                            long n, const vec_l& factors, long e);

// Choose a (pseudo)random A with a trapdoor R s.t. A*[R/I]=G
// A=(aBar| aPrime ) with aPrime=G-aBar*R isn't computed explicitly here
void TDMatrix::initTDmatrix(TDMatrixParams& prms, const CRTmatrix* abar)
{
    FHE_TIMER_START;
    params = &prms;

    if (abar)   // aBar is specified by the caller
        aBarCRT = *abar;
    else        // Choose aBar at random
        randomFill(aBarCRT, prms.n, prms.mBar, prms);

    mat_l sigmaP;
    do
    {
        // Choose R as a random small matrix
        setSmall(trapDoorMatR, params->mBar,params->wLen, params->r);

        // Compute the covariance matrix sigmaP = sigmaX*I - sigmaG*[R/I]*(R^t|I)

        computeCovarianceMatrix(sigmaP, trapDoorMatR,
                                /*sigmaG=*/prms.r *prms.maxFactor, prms.sigmaX);

#ifdef DEBUGPRINT
        long minVal =0;
        long maxVal=0;
        for (long i = 0; i < sigmaP.NumCols(); i++)
            for (long j = 0; j < sigmaP.NumRows(); j++)
            {
                if (sigmaP[i][j] < minVal)
                    minVal = sigmaP[i][j];
                if (sigmaP[i][j] > maxVal)
                    maxVal = sigmaP[i][j];
             }

        ZZ q = prms.getQ();
        cout << "minVal=" << minVal << endl;
        cout << "maxVal=" << maxVal << endl;
        cout << "q=" << q << endl;
        cout << "after computeCovarianceMatrix"  << endl;
#endif // DEBUGPRINT
    }
    while (!(gaussSamp.InitSampler(sigmaP)));  // error, try again

#ifdef DEBUGPRINT
        cout << "after gaussSamp.InitSampler"  << endl;
#endif // DEBUGPRINT
}

// Compute the covariance matrix sigmaP = sigmaX*I - sigmaG*[R/I]*(R^t|I)
// Note that [R/I]*(R^t|I) is a block matrix with this form:
//    RR = ( R^2 | R^t )
//         (  R  |  I  )
void computeCovarianceMatrix(mat_l& covMat, const mat_l& R, long sigmaG,long sigmaX)
{

    FHE_TIMER_START;

    long mBar = R.NumRows();
    long wLen = R.NumCols();
    long m = mBar + wLen;

    covMat.SetDims(m,m); // allocate space

    mat_l tmp(INIT_SIZE, mBar, mBar);
    square(tmp, R);   // Set the top-left mBar-by-mBar as R*R^t

    // Reset the top-left to sigmaX*I - sigmaG*R*R^t
    for (long i=0; i<mBar; i++) for (long j=0; j<mBar; j++)
        {
            if (i==j)
                covMat[i][j] = sigmaX - sigmaG*tmp[i][j];
            else
                covMat[i][j] = -sigmaG*tmp[i][j];
        }

    // Set the top-right ro -sigmaG*R^t and the bootom-left to -sigmaG*R
    for (long i=0; i<mBar; i++) for (long j=0; j<wLen; j++)
        {
            covMat[i][j+mBar] = covMat[j+mBar][i] = -sigmaG * R[i][j];
        }

    // Finally set the bottom-right to (sigmaX-sigmaG)*I;
    for (long i=mBar; i<m; i++) for (long j=mBar; j<m; j++)
        {
            if (i==j) covMat[i][j] = sigmaX -sigmaG;
            else      covMat[i][j] = 0;
        }
}

// Compute A explicitly, A = ( aBar | G - aBar*R )
void TDMatrix::getA(CRTmatrix& A) const
{
    FHE_TIMER_START;
    A.params = params;

    long n = params->n;
    long m = params->m;
    long kFactors = params->kFactors;
    long mBar = params->mBar;

    // Compute A modulo each factor separately

    zz_pPush ppush; // backup NTL's current modulus
    A.SetLength(kFactors);
    for (long f=0; f<kFactors; f++)
    {
        params->zzp_context[f].restore();
        mat_zz_p& aMod = A[f];
        const mat_zz_p& aBarMod = aBarCRT[f];

        aMod.SetDims(n,m); // allocate space

        // Copy aBar
        for (long i=0; i<n; i++) for (long j=0; j<mBar; j++)
            {
                aMod[i][j] = aBarMod[i][j];
            }

        // Compute aBar * R
        mat_zz_p tmp;
        mat_zz_p rMod;
        conv(rMod, trapDoorMatR); // conver to zz_p format
        mul(tmp, aBarMod, rMod); // tmp = aBar * R (mod P_f)

        // Copy -aBar*R to left part of A
        for (long i=0; i<n; i++) for (long j=mBar; j<m; j++)
            {
                aMod[i][j] = -tmp[i][j-mBar];
            }

        // Add G to left part of A
        for (long i=0; i<n; i++)    // one row at a time
        {
            // Add (1, f1, f1^2, ..., f1^e...fk^{e-1}) to each row

            long index = (i * kFactors * params->e) + mBar; // 1st index to add to

            zz_p val = to_zz_p(1L);      // The value to add to the next entry
            for (long j=0; j < kFactors; j++)
                for (long ei=0; ei<params->e; ei++, index++)
                {
                    aMod[i][index] += val;
                    val *= params->factors[j];// Change value to add before next entry
                }
        }
    }

}

/** sampleWithTrapdoor - Implementation of Efficient algorithm SampleD -
 * Algorithm 3 from MP12. Choose p at random and z to match the needed
 * solution. This debugging function also returns p and z.
 **/
int TDMatrix::sampleWithTrapdoor(vec_l &xOut, vec_l& p, vec_l& z,
                                 Vec<vec_zz_p>& U) const
{
    FHE_TIMER_START;
    long mBar = params->mBar;
    long m = params->m;
    long n = params->n;
    long e = params->e;
    long kFactors = params->kFactors;

    //choose perturbation P from (0,sigma)

    Vec<double> zeroVector(INIT_SIZE, m);
    clear(zeroVector);

    gaussSamp.SampleDiscreteGaussian(p, zeroVector);
#ifdef DEBUGPRINT
    cout << "\np = " << p << endl;
#endif
    // Split the pertubation vector in two
    vec_l p1(INIT_SIZE, mBar);
    for (long i = 0; i < mBar; i++)
        p1[i] = p[i];

    vec_l p2(INIT_SIZE, params->wLen);
    for (long i = mBar; i < m; i++)
        p2[i - mBar] = p[i];

    //compute wBar and w

    /* We need to compute v = u-A*p = u-wBar-w (in CRT representation), where
     *
     *    A*p = (Abar| Aprime ) * [p1/p2] = Abar*p1 + Aprime*p2
     *        = Abar*p1 + (G-Abar*R)*p2   = Abar*(p1-R*p2) +  G*p2
     *                                      \___wBar____/    \_w_/
     */

    // p1-R*p2 is computed over the integers, stored in p1
    vec_l tmp(INIT_SIZE, mBar);
    mul(tmp, trapDoorMatR, p2); // R*p2
    p1 -= tmp;                  // p1-R*p2

    // All other vectors are computed modulo each of the factors

    // FIXME: There's room for parallelization here, but in our case we have
    // thousands of calls to sampleWithTrapdoor that can be run in parallel,
    // so there is no reason to parallelize also at this lower level.

    zz_pPush push; // backup NTL's current modulus

    Vec<vec_zz_p> vVec = U; // allocate space and initialize v=u
    vec_zz_p wMod, wBarMod;
    for (long iE = 0; iE < kFactors; iE++)
    {
        params->zzp_context[iE].restore();

        // Compute Abar*(p1-R*p2) modulo the current factor
        const mat_zz_p& aBarMod = aBarCRT[iE];
        mul(wBarMod, aBarMod, conv<vec_zz_p>(p1));  // wBar = Abar*(p1-R*p2)
        // convert p1-R*p2 to vec_zz_p and multiply by aBar modulo each factor

        vVec[iE] -= wBarMod; // u - wBar

        multByG(wMod, p2, n, params->factors, e); // w = G*p2
        vVec[iE] -= wMod;    // u - wBar - w
    }

    // Now that we computed v modulo all the primes, sample z s.t. G*z=v

    sampleG(z, vVec, params);
#ifdef DEBUGPRINT
    cout << "z = " << z << endl;
#endif
    // finally, return x = p + [R/I]*z = p + [R*z / z]

    xOut = p;
    mul(tmp, trapDoorMatR, z); // add R*z to top mBar entries
#ifdef DEBUGPRINT
    cout << "R*z = " << tmp << endl;
#endif
    for (long i=0; i<mBar; i++)
    {
      xOut[i] += tmp[i];
      long absX = abs(xOut[i]);
      if (absX > TDMatrix::maxSample) TDMatrix::maxSample.store(absX);
         // not quite thread-safe, the stored value may not be the maximum
    }

    // add z to bottom m-mBar entries
    for (long i = mBar; i< params->m; i++)
    {
      xOut[i] += z[i - mBar];
      long absX = abs(xOut[i]);
      if (absX > TDMatrix::maxSample) TDMatrix::maxSample.store(absX);
         // not quite thread-safe, the stored value may not be the maximum
    }
    return 0;
}

//Write TDMatrix different parameters to file. The handle to the open
//file is sent as input to the function. returns number of items read by
//the function
long TDMatrix::writeToFile(FILE* handle)
{
    FHE_TIMER_START;
    long count = 0;

    count += params->writeToFile(handle); // write the params
    count += aBarCRT.writeToFile(handle); // write the Abar matrix

    long numRows = trapDoorMatR.NumRows();
    long numCols = trapDoorMatR.NumCols();
    count += fwrite(&numRows,sizeof(long),1,handle);
    count += fwrite(&numCols,sizeof(long),1,handle);
    for (long i = 0; i < numRows; i++)  // write the rows of R, one at a time
    {
        count+= fwrite(trapDoorMatR[i].elts(), sizeof(long), numCols, handle);
        // each row is implemented as a C vector
    }

    count += gaussSamp.writeToFile(handle); // write the DGaussSampler

    return count;
}

//Read the TDMatrix parameters from a flie. The handle to the open file is
//received as input to this function, as well as a pointer to the TDMatrixParams
//if the pointer is not provided and params as not been initialized before, an
//error occurs, as a pointer to this structure is needed
//returns number of items read by the function

long TDMatrix::readFromFile(FILE* handle, TDMatrixParams* prmBuf)
{
    FHE_TIMER_START;
    assert(params != NULL || prmBuf != NULL); // some pointer must be provided
    long count=0;

    // If buffer is given, make params point to it and don't overwrite from iput
    if (prmBuf != NULL)
    {
        TDMatrixParams p;
        count += p.readFromFile(handle);
        assert(p == *prmBuf); // sanity check
        params = prmBuf;      // point to given params
    }
    else
        count += params->readFromFile(handle); // overwrite params from input

    count += aBarCRT.readFromFile(handle,prmBuf); // read the Abar matrix

    long numRows,numCols;
    count += fread(&numRows,sizeof(long),1,handle);
    count += fread(&numCols,sizeof(long),1,handle);
    trapDoorMatR.SetDims(numRows,numCols);
    for (long i = 0; i < numRows; i++)  // write the rows of R, one at a time
    {
        count+= fread(trapDoorMatR[i].elts(), sizeof(long), numCols, handle);
        // each row is implemented as a C vector
    }

    count += gaussSamp.readFromFile(handle); // read the DGaussSampler

    return count;
}

// operator==, //compares all variables in TDMatrix
bool operator==(const TDMatrix& A, const TDMatrix& B)
{
    FHE_TIMER_START;
    if (A.getParams() != B.getParams())
        return false;

    if (A.getABar() != B.getABar())
        return false;

    if (A.getR() != B.getR())
        return false;

    return true;
}

#ifdef DEBUG
static void printTrapdoor(const mat_l& td, const TDMatrixParams& params)
{
    FHE_TIMER_START;
    double maxR = 0;
    double AvgR=1;
    //  long TotalItems;

    for (long i=0; i < params.mBar; i++)
        for (long j=0; j < params.wLen; j++)
        {
            if (abs(td[i][j]) > maxR)
                maxR = abs(td[i][j]);
            AvgR *= (double)(i*params.wLen+j)/(i*params.wLen+j+1);
            AvgR+=(double)abs(td[i][j]) /(i*params.wLen+j+1);
        }
    cout << "maxR = " << maxR << ", AvgR = " << AvgR << "\n";

    mat_l rSquare;
    //how large is R^2?
    square(rSquare,td);
    //  cout << "trapDoorMatR=" << td << "rSquare= " << rSquare << "\n";

    maxR = 0;
    AvgR=1;
    for (long i=0; i < params.mBar; i++)
        for (long j=0; j < params.mBar; j++)
        {
            if (abs(rSquare[i][j]) > maxR)
                maxR = abs(rSquare[i][j]);
            AvgR *= (double)(i*params.wLen+j)/(i*params.wLen+j+1);
            AvgR+=(double)abs(rSquare[i][j]) /(i*params.wLen+j+1);
        }
    cout << "maxRSquare = " << maxR << ", AvgRSquare = " << AvgR << "\n";
}
#endif //DEBUG

ext: cpp

max_stars_repo_name: shaih/BPobfus

max_issues_repo_name: shaih/BPobfus

max_forks_repo_name: shaih/BPobfus


text: #include <boost/spirit/repository/home/karma/directive.hpp>

ext: hpp

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
/**
    Copyright 2001-2019 Red Anchor Trading Co. Ltd.

    Distributed under the Boost Software License, Version 1.0.
    See <http://www.boost.org/LICENSE_1_0.txt>
 */

#ifndef FOST_CONFIG_HPP
#define FOST_CONFIG_HPP
#pragma once

// Some control #defines
#ifdef FOST_OS_WINDOWS
#ifndef WINVER
#define WINVER 0x0500 // Support only Windows 2000 and above
#define _WIN32_WINNT WINVER
#define _WIN32_WINDOWS WINVER
#endif

#define _ATL_APARTMENT_THREADED
#define _ATL_NO_AUTOMATIC_NAMESPACE
#define _ATL_CSTRING_EXPLICIT_CONSTRUCTORS
#define _ATL_ALL_WARNINGS

#define BOOST_ALL_DYN_LINK
#define BOOST_THREAD_USE_DLL

#ifdef _MSC_VER
#include "config.msvc.hpp"
#endif

// Might need to be 1500 for MSVC 9.0
#if (_MSC_VER < 1600)
// In stdint.h on Linux
typedef __int8 int8_t;
typedef unsigned __int8 uint8_t;
typedef __int16 int16_t;
typedef unsigned __int16 uint16_t;
typedef __int32 int32_t;
typedef unsigned __int32 uint32_t;
typedef __int64 int64_t;
typedef unsigned __int64 uint64_t;
#else
#include <stdint.h>
#endif

#define FSL_ABSTRACT __declspec(novtable)
#define FSL_EXPORT __declspec(dllexport)
#else
#include <stdint.h>

typedef long LONG;

#define FSL_ABSTRACT
#define FSL_EXPORT

#undef Null
#endif

#include <limits>
#include <vector>
#include <list>
#include <map>
#include <set>
#include <typeinfo>
#include <string>
#include <sstream>
#include <fstream>

#if INT_MAX != LONG_MAX && LONG_MAX != LLONG_MAX
#define FOST_USE_LONG
#endif

#include <boost/utility.hpp>
#include <boost/smart_ptr.hpp>

#ifdef FOST_OS_WINDOWS
#ifdef FOST_HAVE_MFC
#include <afxwin.h>
#include <afxdisp.h>
#include <atlbase.h>
#include <atlcom.h>
#include <comdef.h>
#else
#define WINDOWS_LEAN_AND_MEAN 1
#define _WINSOCKAPI_
#include <windows.h>
#undef _WINSOCKAPI_
#endif
/// Microsoft are nuts -- the stuff they `#define` It turns out we can't
/// undefine it because loads of their headers rely on this
// #undef interface
#endif

#endif // FOST_CONFIG_HPP

max_stars_repo_name: KayEss/fost-base

max_issues_repo_name: KayEss/fost-base

max_forks_repo_name: KayEss/fost-base


text:
// Boost.Geometry (aka GGL, Generic Geometry Library)
//
// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

#include <geometry_test_common.hpp>

#include <boost/geometry/io/wkt/wkt.hpp>

#include <boost/geometry/algorithms/correct.hpp>
#include <boost/geometry/algorithms/within.hpp>

#include <boost/geometry/geometries/box.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/multi_polygon.hpp>

#include "test_within.hpp"

template <typename P>
void test_all()
{
    typedef bg::model::multi_polygon<bg::model::polygon<P> > mp;

    // test multi-with-one-polygon (trivial case)
    test_geometry<P, mp>("POINT(1 1)", "MULTIPOLYGON(((0 0,0 2,2 2,2 0,0 0)))", true);
    test_geometry<P, mp>("POINT(3 3)", "MULTIPOLYGON(((0 0,0 2,2 2,2 0,0 0)))", false);
    test_geometry<P, mp>("POINT(0 1)", "MULTIPOLYGON(((0 0,0 2,2 2,2 0,0 0)))", false);
    test_geometry<P, mp>("POINT(4 4)", "MULTIPOLYGON(((0 0,0 2,2 2,2 0,0 0)))", false);

    // test if it is in one of them
    std::string multi("MULTIPOLYGON("
        "((0 0,0 2,2 2,2 0,0 0))"
        "((3 3,3 6,6 6,6 3,3 3))"
        ")");
    test_geometry<P, mp>("POINT(4 4)", multi, true);
    test_geometry<P, mp>("POINT(1 1)", multi, true);
    test_geometry<P, mp>("POINT(0 1)", multi, false);
}

int test_main( int , char* [] )
{
    //test_all<bg::model::d2::point_xy<int> >();
    test_all<bg::model::d2::point_xy<double> >();

#if defined(HAVE_TTMATH)
    test_all<bg::model::d2::point_xy<ttmath_big> >();
#endif

    return 0;
}


text:
// Copyright 2012 Cloudera Inc.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#include "exec/aggregation-node.h"

#include <math.h>
#include <sstream>
#include <boost/functional/hash.hpp>
#include <thrift/protocol/TDebugProtocol.h>

#include <x86intrin.h>

#include "codegen/codegen-anyval.h"
#include "codegen/llvm-codegen.h"
#include "exec/hash-table.inline.h"
#include "exprs/agg-fn-evaluator.h"
#include "exprs/expr.h"
#include "runtime/descriptors.h"
#include "runtime/mem-pool.h"
#include "runtime/raw-value.h"
#include "runtime/row-batch.h"
#include "runtime/runtime-state.h"
#include "runtime/string-value.inline.h"
#include "runtime/tuple.h"
#include "runtime/tuple-row.h"
#include "udf/udf-internal.h"
#include "util/debug-util.h"
#include "util/runtime-profile.h"

#include "gen-cpp/Exprs_types.h"
#include "gen-cpp/PlanNodes_types.h"

using namespace impala;
using namespace std;
using namespace boost;
using namespace llvm;

namespace impala {

const char* AggregationNode::LLVM_CLASS_NAME = "class.impala::AggregationNode";

// TODO: pass in maximum size; enforce by setting limit in mempool
AggregationNode::AggregationNode(ObjectPool* pool, const TPlanNode& tnode,
                                 const DescriptorTbl& descs)
  : ExecNode(pool, tnode, descs),
    agg_tuple_id_(tnode.agg_node.agg_tuple_id),
    agg_tuple_desc_(NULL),
    singleton_output_tuple_(NULL),
    codegen_process_row_batch_fn_(NULL),
    process_row_batch_fn_(NULL),
    is_merge_(tnode.agg_node.__isset.is_merge ? tnode.agg_node.is_merge : false),
    needs_finalize_(tnode.agg_node.need_finalize),
    build_timer_(NULL),
    get_results_timer_(NULL),
    hash_table_buckets_counter_(NULL) {
}

Status AggregationNode::Init(const TPlanNode& tnode) {
  RETURN_IF_ERROR(ExecNode::Init(tnode));
  RETURN_IF_ERROR(
      Expr::CreateExprTrees(pool_, tnode.agg_node.grouping_exprs, &probe_exprs_));
  for (int i = 0; i < tnode.agg_node.aggregate_functions.size(); ++i) {
    AggFnEvaluator* evaluator;
    RETURN_IF_ERROR(AggFnEvaluator::Create(
        pool_, tnode.agg_node.aggregate_functions[i], &evaluator));
    aggregate_evaluators_.push_back(evaluator);
  }
  return Status::OK;
}

Status AggregationNode::Prepare(RuntimeState* state) {
  SCOPED_TIMER(runtime_profile_->total_time_counter());
  RETURN_IF_ERROR(ExecNode::Prepare(state));

  tuple_pool_.reset(new MemPool(mem_tracker()));
  build_timer_ = ADD_TIMER(runtime_profile(), "BuildTime");
  get_results_timer_ = ADD_TIMER(runtime_profile(), "GetResultsTime");
  hash_table_buckets_counter_ =
      ADD_COUNTER(runtime_profile(), "BuildBuckets", TCounterType::UNIT);
  hash_table_load_factor_counter_ =
      ADD_COUNTER(runtime_profile(), "LoadFactor", TCounterType::DOUBLE_VALUE);

  agg_tuple_desc_ = state->desc_tbl().GetTupleDescriptor(agg_tuple_id_);
  RETURN_IF_ERROR(Expr::Prepare(probe_exprs_, state, child(0)->row_desc(), false));

  // Construct build exprs from agg_tuple_desc_
  for (int i = 0; i < probe_exprs_.size(); ++i) {
    SlotDescriptor* desc = agg_tuple_desc_->slots()[i];
    Expr* expr = new SlotRef(desc);
    state->obj_pool()->Add(expr);
    build_exprs_.push_back(expr);
  }
  RETURN_IF_ERROR(Expr::Prepare(build_exprs_, state, row_desc(), false));

  int j = probe_exprs_.size();
  for (int i = 0; i < aggregate_evaluators_.size(); ++i, ++j) {
    // skip non-materialized slots; we don't have evaluators instantiated for those
    while (!agg_tuple_desc_->slots()[j]->is_materialized()) {
      DCHECK_LT(j, agg_tuple_desc_->slots().size() - 1)
          << "#eval= " << aggregate_evaluators_.size()
          << " #probe=" << probe_exprs_.size();
      ++j;
    }
    SlotDescriptor* desc = agg_tuple_desc_->slots()[j];
    RETURN_IF_ERROR(aggregate_evaluators_[i]->Prepare(state, child(0)->row_desc(), desc));
  }

  // TODO: how many buckets?
  hash_tbl_.reset(new HashTable(state, build_exprs_, probe_exprs_, 1, true, true,
      id(), mem_tracker()));

  if (probe_exprs_.empty()) {
    // create single output tuple now; we need to output something
    // even if our input is empty
    singleton_output_tuple_ = ConstructAggTuple();
    hash_tbl_->Insert(reinterpret_cast<TupleRow*>(&singleton_output_tuple_));
    output_iterator_ = hash_tbl_->Begin();
  }

  if (state->codegen_enabled()) {
    DCHECK(state->codegen() != NULL);
    Function* update_tuple_fn = CodegenUpdateAggTuple(state->codegen());
    if (update_tuple_fn != NULL) {
      codegen_process_row_batch_fn_ =
          CodegenProcessRowBatch(state->codegen(), update_tuple_fn);
      if (codegen_process_row_batch_fn_ != NULL) {
        // Update to using codegen'd process row batch.
        state->codegen()->AddFunctionToJit(
            codegen_process_row_batch_fn_,
            reinterpret_cast<void**>(&process_row_batch_fn_));
        AddRuntimeExecOption("Codegen Enabled");
      }
    }
  }
  return Status::OK;
}

Status AggregationNode::Open(RuntimeState* state) {
  SCOPED_TIMER(runtime_profile_->total_time_counter());
  RETURN_IF_ERROR(ExecNode::Open(state));

  RETURN_IF_ERROR(Expr::Open(probe_exprs_, state));
  RETURN_IF_ERROR(Expr::Open(build_exprs_, state));

  for (int i = 0; i < aggregate_evaluators_.size(); ++i) {
    RETURN_IF_ERROR(aggregate_evaluators_[i]->Open(state));
  }

  RETURN_IF_ERROR(children_[0]->Open(state));

  RowBatch batch(children_[0]->row_desc(), state->batch_size(), mem_tracker());
  int64_t num_input_rows = 0;
  while (true) {
    bool eos;
    RETURN_IF_CANCELLED(state);
    RETURN_IF_ERROR(state->CheckQueryState());
    RETURN_IF_ERROR(children_[0]->GetNext(state, &batch, &eos));
    SCOPED_TIMER(build_timer_);

    if (VLOG_ROW_IS_ON) {
      for (int i = 0; i < batch.num_rows(); ++i) {
        TupleRow* row = batch.GetRow(i);
        VLOG_ROW << "input row: " << PrintRow(row, children_[0]->row_desc());
      }
    }
    if (process_row_batch_fn_ != NULL) {
      process_row_batch_fn_(this, &batch);
    } else if (probe_exprs_.empty()) {
      ProcessRowBatchNoGrouping(&batch);
    } else {
      ProcessRowBatchWithGrouping(&batch);
    }
    COUNTER_SET(hash_table_buckets_counter_, hash_tbl_->num_buckets());
    COUNTER_SET(hash_table_load_factor_counter_, hash_tbl_->load_factor());
    num_input_rows += batch.num_rows();
    // We must set output_iterator_ here, rather than outside the loop, because
    // output_iterator_ must be set if the function returns within the loop
    output_iterator_ = hash_tbl_->Begin();

    batch.Reset();
    RETURN_IF_ERROR(state->CheckQueryState());
    if (eos) break;
  }

  // We have consumed all of the input from the child and transfered ownership of the
  // resources we need, so the child can be closed safely to release its resources.
  child(0)->Close(state);
  VLOG_FILE << "aggregated " << num_input_rows << " input rows into "
            << hash_tbl_->size() << " output rows";
  return Status::OK;
}

Status AggregationNode::GetNext(RuntimeState* state, RowBatch* row_batch, bool* eos) {
  SCOPED_TIMER(runtime_profile_->total_time_counter());
  RETURN_IF_ERROR(ExecDebugAction(TExecNodePhase::GETNEXT, state));
  RETURN_IF_CANCELLED(state);
  RETURN_IF_ERROR(state->CheckQueryState());
  SCOPED_TIMER(get_results_timer_);

  if (ReachedLimit()) {
    *eos = true;
    return Status::OK;
  }
  Expr** conjuncts = &conjuncts_[0];
  int num_conjuncts = conjuncts_.size();

  while (!output_iterator_.AtEnd() && !row_batch->AtCapacity()) {
    int row_idx = row_batch->AddRow();
    TupleRow* row = row_batch->GetRow(row_idx);
    Tuple* agg_tuple = output_iterator_.GetRow()->GetTuple(0);
    FinalizeAggTuple(agg_tuple);
    output_iterator_.Next<false>();
    row->SetTuple(0, agg_tuple);
    if (ExecNode::EvalConjuncts(conjuncts, num_conjuncts, row)) {
      VLOG_ROW << "output row: " << PrintRow(row, row_desc());
      row_batch->CommitLastRow();
      ++num_rows_returned_;
      if (ReachedLimit()) break;
    }
  }
  *eos = output_iterator_.AtEnd() || ReachedLimit();
  COUNTER_SET(rows_returned_counter_, num_rows_returned_);
  return Status::OK;
}

void AggregationNode::Close(RuntimeState* state) {
  if (is_closed()) return;

  // Iterate through the remaining rows in the hash table and call Serialize/Finalize on
  // them in order to free any memory allocated by UDAs
  while (!output_iterator_.AtEnd()) {
    Tuple* agg_tuple = output_iterator_.GetRow()->GetTuple(0);
    FinalizeAggTuple(agg_tuple);
    output_iterator_.Next<false>();
  }

  if (tuple_pool_.get() != NULL) tuple_pool_->FreeAll();
  if (hash_tbl_.get() != NULL) hash_tbl_->Close();
  for (int i = 0; i < aggregate_evaluators_.size(); ++i) {
    aggregate_evaluators_[i]->Close(state);
  }
  Expr::Close(probe_exprs_, state);
  Expr::Close(build_exprs_, state);
  ExecNode::Close(state);
}

Tuple* AggregationNode::ConstructAggTuple() {
  Tuple* agg_tuple = Tuple::Create(agg_tuple_desc_->byte_size(), tuple_pool_.get());
  vector<SlotDescriptor*>::const_iterator slot_desc = agg_tuple_desc_->slots().begin();

  // copy grouping values
  for (int i = 0; i < probe_exprs_.size(); ++i, ++slot_desc) {
    if (hash_tbl_->last_expr_value_null(i)) {
      agg_tuple->SetNull((*slot_desc)->null_indicator_offset());
    } else {
      void* src = hash_tbl_->last_expr_value(i);
      void* dst = agg_tuple->GetSlot((*slot_desc)->tuple_offset());
      RawValue::Write(src, dst, (*slot_desc)->type(), tuple_pool_.get());
    }
  }

  // Initialize aggregate output.
  for (int i = 0; i < aggregate_evaluators_.size(); ++i, ++slot_desc) {
    while (!(*slot_desc)->is_materialized()) ++slot_desc;
    AggFnEvaluator* evaluator = aggregate_evaluators_[i];
    evaluator->Init(agg_tuple);
    // Codegen specific path.
    // To minimize branching on the UpdateAggTuple path, initialize the result value
    // so that UpdateAggTuple doesn't have to check if the aggregation
    // dst slot is null.
    //  - sum/count: 0
    //  - min: max_value
    //  - max: min_value
    // TODO: remove when we don't use the irbuilder for codegen here.
    // This optimization no longer applies with AnyVal
    if ((*slot_desc)->type().type != TYPE_STRING &&
        (*slot_desc)->type().type != TYPE_TIMESTAMP &&
        (*slot_desc)->type().type != TYPE_CHAR &&
        (*slot_desc)->type().type != TYPE_DECIMAL) {
      ExprValue default_value;
      void* default_value_ptr = NULL;
      switch (evaluator->agg_op()) {
        case AggFnEvaluator::MIN:
          default_value_ptr = default_value.SetToMax((*slot_desc)->type());
          RawValue::Write(default_value_ptr, agg_tuple, *slot_desc, NULL);
          break;
        case AggFnEvaluator::MAX:
          default_value_ptr = default_value.SetToMin((*slot_desc)->type());
          RawValue::Write(default_value_ptr, agg_tuple, *slot_desc, NULL);
          break;
        default:
          break;
      }
    }
  }
  return agg_tuple;
}

void AggregationNode::UpdateAggTuple(Tuple* tuple, TupleRow* row) {
  DCHECK(tuple != NULL || aggregate_evaluators_.empty());
  for (vector<AggFnEvaluator*>::const_iterator evaluator = aggregate_evaluators_.begin();
      evaluator != aggregate_evaluators_.end(); ++evaluator) {
    if (is_merge_) {
      (*evaluator)->Merge(row, tuple);
    } else {
      (*evaluator)->Update(row, tuple);
    }
  }
}

void AggregationNode::FinalizeAggTuple(Tuple* tuple) {
  DCHECK(tuple != NULL || aggregate_evaluators_.empty());
  for (vector<AggFnEvaluator*>::const_iterator evaluator = aggregate_evaluators_.begin();
      evaluator != aggregate_evaluators_.end(); ++evaluator) {
    if (needs_finalize_) {
      (*evaluator)->Finalize(tuple);
    } else {
      (*evaluator)->Serialize(tuple);
    }
  }
}

void AggregationNode::DebugString(int indentation_level, stringstream* out) const {
  *out << string(indentation_level * 2, ' ');
  *out << "AggregationNode(tuple_id=" << agg_tuple_id_
       << " is_merge=" << is_merge_ << " needs_finalize=" << needs_finalize_
       << " probe_exprs=" << Expr::DebugString(probe_exprs_)
       << " agg_exprs=" << AggFnEvaluator::DebugString(aggregate_evaluators_);
  ExecNode::DebugString(indentation_level, out);
  *out << ")";
}

IRFunction::Type GetHllUpdateFunction(const ColumnType& type) {
  switch (type.type) {
    case TYPE_BOOLEAN: return IRFunction::HLL_UPDATE_BOOLEAN;
    case TYPE_TINYINT: return IRFunction::HLL_UPDATE_TINYINT;
    case TYPE_SMALLINT: return IRFunction::HLL_UPDATE_SMALLINT;
    case TYPE_INT: return IRFunction::HLL_UPDATE_INT;
    case TYPE_BIGINT: return IRFunction::HLL_UPDATE_BIGINT;
    case TYPE_FLOAT: return IRFunction::HLL_UPDATE_FLOAT;
    case TYPE_DOUBLE: return IRFunction::HLL_UPDATE_DOUBLE;
    case TYPE_STRING: return IRFunction::HLL_UPDATE_STRING;
    case TYPE_DECIMAL: return IRFunction::HLL_UPDATE_DECIMAL;
    default:
      DCHECK(false) << "Unsupported type: " << type;
      return IRFunction::FN_END;
  }
}

// IR Generation for updating a single aggregation slot. Signature is:
// void UpdateSlot(FunctionContext* fn_ctx, AggTuple* agg_tuple, char** row)
//
// The IR for sum(double_col) is:
//   define void @UpdateSlot(%"class.impala_udf::FunctionContext"* %fn_ctx,
//                           { i8, double }* %agg_tuple, i8** %row) {
//   entry:
//     %src_null_ptr = alloca i1
//     %src_value = call double @SlotRef(i8** %row, i8* null, i1* %src_null_ptr)
//     %child_null = load i1* %src_null_ptr
//     br i1 %child_null, label %ret, label %src_not_null
//
//   src_not_null:                                     ; preds = %entry
//     %dst_slot_ptr = getelementptr inbounds { i8, double }* %agg_tuple, i32 0, i32 1
//     call void @SetNotNull({ i8, double }* %agg_tuple)
//     %dst_val = load double* %dst_slot_ptr
//     %0 = fadd double %dst_val, %src_value
//     store double %0, double* %dst_slot_ptr
//     br label %ret
//
//   ret:                                              ; preds = %src_not_null, %entry
//     ret void
//   }
//
// The IR for ndv(double_col) is:
//   define void @UpdateSlot(%"class.impala_udf::FunctionContext"* %fn_ctx,
//                           { i8, %"struct.impala::StringValue" }* %agg_tuple,
//                           i8** %row) {
//   entry:
//     %src_null_ptr = alloca i1
//     %src_value = call double @SlotRef(i8** %row, i8* null, i1* %src_null_ptr)
//     %child_null = load i1* %src_null_ptr
//     br i1 %child_null, label %ret, label %src_not_null
//
//   src_not_null:                                     ; preds = %entry
//     %dst_slot_ptr = getelementptr inbounds
//                     { i8, %"struct.impala::StringValue" }* %agg_tuple, i32 0, i32 1
//     call void @SetNotNull({ i8, %"struct.impala::StringValue" }* %agg_tuple)
//     %dst_val = load %"struct.impala::StringValue"* %dst_slot_ptr
//     %src_anyval = insertvalue { i8, double } zeroinitializer, double %src_value, 1
//     %src_lowered_ptr = alloca { i8, double }
//     store { i8, double } %src_anyval, { i8, double }* %src_lowered_ptr
//     %src_unlowered_ptr = bitcast { i8, double }* %src_lowered_ptr to
//                          %"struct.impala_udf::DoubleVal"*
//     %ptr = extractvalue %"struct.impala::StringValue" %dst_val, 0
//     %dst_stringval = insertvalue { i64, i8* } zeroinitializer, i8* %ptr, 1
//     %len = extractvalue %"struct.impala::StringValue" %dst_val, 1
//     %0 = extractvalue { i64, i8* } %dst_stringval, 0
//     %1 = zext i32 %len to i64
//     %2 = shl i64 %1, 32
//     %3 = and i64 %0, 0
//     %4 = or i64 %3, %2
//     %dst_stringval1 = insertvalue { i64, i8* } %dst_stringval, i64 %4, 0
//     %dst_lowered_ptr = alloca { i64, i8* }
//     store { i64, i8* } %dst_stringval1, { i64, i8* }* %dst_lowered_ptr
//     %dst_unlowered_ptr = bitcast { i64, i8* }* %dst_lowered_ptr to
//                          %"struct.impala_udf::StringVal"*
//     call void @HllUpdate(%"class.impala_udf::FunctionContext"* %fn_ctx,
//                          %"struct.impala_udf::DoubleVal"* %src_unlowered_ptr,
//                          %"struct.impala_udf::StringVal"* %dst_unlowered_ptr)
//     %anyval_result = load { i64, i8* }* %dst_lowered_ptr
//     %5 = extractvalue { i64, i8* } %anyval_result, 1
//     %6 = insertvalue %"struct.impala::StringValue" zeroinitializer, i8* %5, 0
//     %7 = extractvalue { i64, i8* } %anyval_result, 0
//     %8 = ashr i64 %7, 32
//     %9 = trunc i64 %8 to i32
//     %10 = insertvalue %"struct.impala::StringValue" %6, i32 %9, 1
//     store %"struct.impala::StringValue" %10,
//           %"struct.impala::StringValue"* %dst_slot_ptr
//     br label %ret
//
//   ret:                                              ; preds = %src_not_null, %entry
//     ret void
//   }
llvm::Function* AggregationNode::CodegenUpdateSlot(
    LlvmCodeGen* codegen, AggFnEvaluator* evaluator, SlotDescriptor* slot_desc) {
  int field_idx = slot_desc->field_idx();
  DCHECK(slot_desc->is_materialized());

  LLVMContext& context = codegen->context();

  PointerType* fn_ctx_type =
      codegen->GetPtrType(FunctionContextImpl::LLVM_FUNCTIONCONTEXT_NAME);
  StructType* tuple_struct = agg_tuple_desc_->GenerateLlvmStruct(codegen);
  PointerType* tuple_ptr_type = PointerType::get(tuple_struct, 0);
  PointerType* ptr_type = codegen->ptr_type();

  // Create UpdateSlot prototype
  LlvmCodeGen::FnPrototype prototype(codegen, "UpdateSlot", codegen->void_type());
  prototype.AddArgument(LlvmCodeGen::NamedVariable("fn_ctx", fn_ctx_type));
  prototype.AddArgument(LlvmCodeGen::NamedVariable("agg_tuple", tuple_ptr_type));
  prototype.AddArgument(LlvmCodeGen::NamedVariable("row", PointerType::get(ptr_type, 0)));

  LlvmCodeGen::LlvmBuilder builder(context);
  Value* args[3];
  Function* fn = prototype.GeneratePrototype(&builder, &args[0]);

  Value* fn_ctx_arg = args[0];
  Value* agg_tuple_arg = args[1];
  Value* row_arg = args[2];

  LlvmCodeGen::NamedVariable null_var("src_null_ptr", codegen->boolean_type());
  Value* src_is_null_ptr = codegen->CreateEntryBlockAlloca(fn, null_var);

  // Call expr function to get src slot value
  DCHECK_EQ(evaluator->input_exprs().size(), 1);
  Expr* input_expr = evaluator->input_exprs()[0];
  Function* agg_expr_fn = input_expr->codegen_fn();
  int scratch_buffer_size = input_expr->scratch_buffer_size();
  DCHECK_EQ(scratch_buffer_size, 0);
  DCHECK(agg_expr_fn != NULL);
  if (agg_expr_fn == NULL) return NULL;

  BasicBlock* src_not_null_block, *ret_block;
  codegen->CreateIfElseBlocks(fn, "src_not_null", "ret", &src_not_null_block, &ret_block);

  Value* expr_args[] = { row_arg, ConstantPointerNull::get(ptr_type), src_is_null_ptr };
  Value* src_value = input_expr->CodegenGetValue(codegen, builder.GetInsertBlock(),
      expr_args, ret_block, src_not_null_block, "src_value");

  // Src slot is not null, update dst_slot
  builder.SetInsertPoint(src_not_null_block);
  Value* dst_ptr = builder.CreateStructGEP(agg_tuple_arg, field_idx, "dst_slot_ptr");
  Value* result = NULL;

  if (slot_desc->is_nullable()) {
    // Dst is NULL, just update dst slot to src slot and clear null bit
    Function* clear_null_fn = slot_desc->CodegenUpdateNull(codegen, tuple_struct, false);
    builder.CreateCall(clear_null_fn, agg_tuple_arg);
  }

  // Update the slot
  Value* dst_value = builder.CreateLoad(dst_ptr, "dst_val");
  switch (evaluator->agg_op()) {
    case AggFnEvaluator::COUNT:
      result = builder.CreateAdd(dst_value,
          codegen->GetIntConstant(TYPE_BIGINT, 1), "count_inc");
      break;
    case AggFnEvaluator::MIN: {
      Function* min_fn = codegen->CodegenMinMax(slot_desc->type(), true);
      Value* min_args[] = { dst_value, src_value };
      result = builder.CreateCall(min_fn, min_args, "min_value");
      break;
    }
    case AggFnEvaluator::MAX: {
      Function* max_fn = codegen->CodegenMinMax(slot_desc->type(), false);
      Value* max_args[] = { dst_value, src_value };
      result = builder.CreateCall(max_fn, max_args, "max_value");
      break;
    }
    case AggFnEvaluator::SUM:
      if (slot_desc->type().type == TYPE_FLOAT || slot_desc->type().type == TYPE_DOUBLE) {
        result = builder.CreateFAdd(dst_value, src_value);
      } else {
        result = builder.CreateAdd(dst_value, src_value);
      }
      break;
    case AggFnEvaluator::NDV: {
      DCHECK_EQ(slot_desc->type().type, TYPE_STRING);
      IRFunction::Type ir_function_type = is_merge_ ? IRFunction::HLL_MERGE
                                          : GetHllUpdateFunction(input_expr->type());
      Function* hll_fn = codegen->GetFunction(ir_function_type);

      // Convert src_value to *Val src_anyval. src_value is the returned value of a
      // codegen'd expr compute function, so either is a native type or a StringVal*.
      CodegenAnyVal src_anyval = CodegenAnyVal::GetNonNullVal(
          codegen, &builder, input_expr->type(), "src_anyval");
      if (src_value->getType()->isPointerTy()) {
        DCHECK_EQ(src_value->getType(), codegen->GetPtrType(TYPE_STRING))
            << endl << LlvmCodeGen::Print(src_value);
        src_anyval.SetFromRawPtr(src_value);
      } else {
        src_anyval.SetFromRawValue(src_value);
      }

      // Create pointer to src_anyval to pass to HllUpdate() function. We must use the
      // unlowered type.
      Value* src_lowered_ptr = codegen->CreateEntryBlockAlloca(
          fn, LlvmCodeGen::NamedVariable("src_lowered_ptr",
                                         src_anyval.value()->getType()));
      builder.CreateStore(src_anyval.value(), src_lowered_ptr);
      Type* unlowered_ptr_type =
          CodegenAnyVal::GetUnloweredType(codegen, input_expr->type())->getPointerTo();
      Value* src_unlowered_ptr =
          builder.CreateBitCast(src_lowered_ptr, unlowered_ptr_type, "src_unlowered_ptr");

      // Create StringVal* intermediate argument from dst_value
      CodegenAnyVal dst_stringval = CodegenAnyVal::GetNonNullVal(
          codegen, &builder, TYPE_STRING, "dst_stringval");
      dst_stringval.SetFromRawValue(dst_value);
      // Create pointer to dst_stringval to pass to HllUpdate() function. We must use
      // the unlowered type.
      Value* dst_lowered_ptr = codegen->CreateEntryBlockAlloca(
          fn, LlvmCodeGen::NamedVariable("dst_lowered_ptr",
                                         dst_stringval.value()->getType()));
      builder.CreateStore(dst_stringval.value(), dst_lowered_ptr);
      unlowered_ptr_type =
          codegen->GetPtrType(CodegenAnyVal::GetUnloweredType(codegen, TYPE_STRING));
      Value* dst_unlowered_ptr =
          builder.CreateBitCast(dst_lowered_ptr, unlowered_ptr_type, "dst_unlowered_ptr");

      // Call 'hll_fn'
      builder.CreateCall3(hll_fn, fn_ctx_arg, src_unlowered_ptr, dst_unlowered_ptr);

      // Convert StringVal intermediate 'dst_arg' back to StringValue
      Value* anyval_result = builder.CreateLoad(dst_lowered_ptr, "anyval_result");
      result = CodegenAnyVal(codegen, &builder, TYPE_STRING, anyval_result)
               .ToRawValue();
      break;
    }
    default:
      DCHECK(false) << "bad aggregate operator: " << evaluator->agg_op();
  }

  builder.CreateStore(result, dst_ptr);
  builder.CreateBr(ret_block);

  builder.SetInsertPoint(ret_block);
  builder.CreateRetVoid();

  return codegen->FinalizeFunction(fn);
}

// IR codegen for the UpdateAggTuple loop.  This loop is query specific and
// based on the aggregate functions.  The function signature must match the non-
// codegen'd UpdateAggTuple exactly.
// For the query:
// select count(*), count(int_col), sum(double_col) the IR looks like:
//
// define void @UpdateAggTuple(%"class.impala::AggregationNode"* %this_ptr,
//                             %"class.impala::Tuple"* %agg_tuple,
//                             %"class.impala::TupleRow"* %tuple_row) {
// entry:
//   %tuple = bitcast %"class.impala::AggregationNode"* %agg_tuple to
//                                              { i8, i64, i64, double }*
//   %row = bitcast %"class.impala::TupleRow"* %tuple_row to i8**
//   %src_slot = getelementptr inbounds { i8, i64, i64, double }* %tuple, i32 0, i32 2
//   %count_star_val = load i64* %src_slot
//   %count_star_inc = add i64 %count_star_val, 1
//   store i64 %count_star_inc, i64* %src_slot
//   call void @UpdateSlot({ i8, i64, i64, double }* %tuple, i8** %row)
//   call void @UpdateSlot2({ i8, i64, i64, double }* %tuple, i8** %row)
//   ret void
// }
Function* AggregationNode::CodegenUpdateAggTuple(LlvmCodeGen* codegen) {
  SCOPED_TIMER(codegen->codegen_timer());

  for (int i = 0; i < probe_exprs_.size(); ++i) {
    if (probe_exprs_[i]->codegen_fn() == NULL) {
      VLOG_QUERY << "Could not codegen UpdateAggTuple because "
                 << "codegen for the grouping exprs is not yet supported.";
      return NULL;
    }
  }

  int j = probe_exprs_.size();
  for (int i = 0; i < aggregate_evaluators_.size(); ++i, ++j) {
    // skip non-materialized slots; we don't have evaluators instantiated for those
    while (!agg_tuple_desc_->slots()[j]->is_materialized()) {
      DCHECK_LT(j, agg_tuple_desc_->slots().size() - 1);
      ++j;
    }
    SlotDescriptor* slot_desc = agg_tuple_desc_->slots()[j];
    AggFnEvaluator* evaluator = aggregate_evaluators_[i];

    // Timestamp and char are never supported. NDV supports decimal and string but no
    // other functions.
    // TODO: the other aggregate functions might work with decimal as-is
    if (slot_desc->type().type == TYPE_TIMESTAMP || slot_desc->type().type == TYPE_CHAR ||
        (evaluator->agg_op() != AggFnEvaluator::NDV &&
         (slot_desc->type().type == TYPE_DECIMAL ||
          slot_desc->type().type == TYPE_STRING))) {
      VLOG_QUERY << "Could not codegen UpdateAggTuple because "
                 << "string, char, timestamp and decimal are not yet supported.";
      return NULL;
    }

    // If the evaluator can't be generated, bail generating this function
    if (!evaluator->is_count_star() &&
        evaluator->input_exprs()[0]->codegen_fn() == NULL) {
      VLOG_QUERY << "Could not codegen UpdateAggTuple because the "
                 << "underlying exprs cannot be codegened.";
      return NULL;
    }

    // Don't codegen things that aren't builtins (for now)
    if (!evaluator->is_builtin()) return NULL;
  }

  if (agg_tuple_desc_->GenerateLlvmStruct(codegen) == NULL) {
    VLOG_QUERY << "Could not codegen UpdateAggTuple because we could"
               << "not generate a  matching llvm struct for the result agg tuple.";
    return NULL;
  }

  // Get the types to match the UpdateAggTuple signature
  Type* agg_node_type = codegen->GetType(AggregationNode::LLVM_CLASS_NAME);
  Type* agg_tuple_type = codegen->GetType(Tuple::LLVM_CLASS_NAME);
  Type* tuple_row_type = codegen->GetType(TupleRow::LLVM_CLASS_NAME);

  DCHECK(agg_node_type != NULL);
  DCHECK(agg_tuple_type != NULL);
  DCHECK(tuple_row_type != NULL);

  PointerType* agg_node_ptr_type = PointerType::get(agg_node_type, 0);
  PointerType* agg_tuple_ptr_type = PointerType::get(agg_tuple_type, 0);
  PointerType* tuple_row_ptr_type = PointerType::get(tuple_row_type, 0);

  // The signature needs to match the non-codegen'd signature exactly.
  PointerType* ptr_type = codegen->ptr_type();
  StructType* tuple_struct = agg_tuple_desc_->GenerateLlvmStruct(codegen);
  PointerType* tuple_ptr = PointerType::get(tuple_struct, 0);
  LlvmCodeGen::FnPrototype prototype(codegen, "UpdateAggTuple", codegen->void_type());
  prototype.AddArgument(LlvmCodeGen::NamedVariable("this_ptr", agg_node_ptr_type));
  prototype.AddArgument(LlvmCodeGen::NamedVariable("agg_tuple", agg_tuple_ptr_type));
  prototype.AddArgument(LlvmCodeGen::NamedVariable("tuple_row", tuple_row_ptr_type));

  LlvmCodeGen::LlvmBuilder builder(codegen->context());
  Value* args[3];
  Function* fn = prototype.GeneratePrototype(&builder, &args[0]);

  // Cast the parameter types to the internal llvm runtime types.
  args[1] = builder.CreateBitCast(args[1], tuple_ptr, "tuple");
  args[2] = builder.CreateBitCast(args[2], PointerType::get(ptr_type, 0), "row");

  // Loop over each expr and generate the IR for that slot.  If the expr is not
  // count(*), generate a helper IR function to update the slot and call that.
  j = probe_exprs_.size();
  for (int i = 0; i < aggregate_evaluators_.size(); ++i, ++j) {
    // skip non-materialized slots; we don't have evaluators instantiated for those
    while (!agg_tuple_desc_->slots()[j]->is_materialized()) {
      DCHECK_LT(j, agg_tuple_desc_->slots().size() - 1);
      ++j;
    }
    SlotDescriptor* slot_desc = agg_tuple_desc_->slots()[j];
    AggFnEvaluator* evaluator = aggregate_evaluators_[i];
    if (evaluator->is_count_star()) {
      // TODO: we should be able to hoist this up to the loop over the batch and just
      // increment the slot by the number of rows in the batch.
      int field_idx = slot_desc->field_idx();
      Value* const_one = codegen->GetIntConstant(TYPE_BIGINT, 1);
      Value* slot_ptr = builder.CreateStructGEP(args[1], field_idx, "src_slot");
      Value* slot_loaded = builder.CreateLoad(slot_ptr, "count_star_val");
      Value* count_inc = builder.CreateAdd(slot_loaded, const_one, "count_star_inc");
      builder.CreateStore(count_inc, slot_ptr);
    } else {
      Function* update_slot_fn = CodegenUpdateSlot(codegen, evaluator, slot_desc);
      if (update_slot_fn == NULL) return NULL;
      Value* fn_ctx_arg = codegen->CastPtrToLlvmPtr(
          codegen->GetPtrType(FunctionContextImpl::LLVM_FUNCTIONCONTEXT_NAME),
          evaluator->ctx());
      builder.CreateCall3(update_slot_fn, fn_ctx_arg, args[1], args[2]);
    }
  }
  builder.CreateRetVoid();

  // CodegenProcessRowBatch() does the final optimizations.
  return codegen->FinalizeFunction(fn);
}

Function* AggregationNode::CodegenProcessRowBatch(
    LlvmCodeGen* codegen, Function* update_tuple_fn) {
  SCOPED_TIMER(codegen->codegen_timer());
  DCHECK(update_tuple_fn != NULL);

  // Get the cross compiled update row batch function
  IRFunction::Type ir_fn = (!probe_exprs_.empty() ?
      IRFunction::AGG_NODE_PROCESS_ROW_BATCH_WITH_GROUPING :
      IRFunction::AGG_NODE_PROCESS_ROW_BATCH_NO_GROUPING);
  Function* process_batch_fn = codegen->GetFunction(ir_fn);

  if (process_batch_fn == NULL) {
    LOG(ERROR) << "Could not find AggregationNode::ProcessRowBatch in module.";
    return NULL;
  }

  int replaced = 0;
  if (!probe_exprs_.empty()) {
    // Aggregation w/o grouping does not use a hash table.

    // Codegen for hash
    Function* hash_fn = hash_tbl_->CodegenHashCurrentRow(codegen);
    if (hash_fn == NULL) return NULL;

    // Codegen HashTable::Equals
    Function* equals_fn = hash_tbl_->CodegenEquals(codegen);
    if (equals_fn == NULL) return NULL;

    // Codegen for evaluating build rows
    Function* eval_build_row_fn = hash_tbl_->CodegenEvalTupleRow(codegen, true);
    if (eval_build_row_fn == NULL) return NULL;

    // Codegen for evaluating probe rows
    Function* eval_probe_row_fn = hash_tbl_->CodegenEvalTupleRow(codegen, false);
    if (eval_probe_row_fn == NULL) return NULL;

    // Replace call sites
    process_batch_fn = codegen->ReplaceCallSites(process_batch_fn, false,
        eval_build_row_fn, "EvalBuildRow", &replaced);
    DCHECK_EQ(replaced, 1);

    process_batch_fn = codegen->ReplaceCallSites(process_batch_fn, false,
        eval_probe_row_fn, "EvalProbeRow", &replaced);
    DCHECK_EQ(replaced, 1);

    process_batch_fn = codegen->ReplaceCallSites(process_batch_fn, false,
        hash_fn, "HashCurrentRow", &replaced);
    DCHECK_EQ(replaced, 2);

    process_batch_fn = codegen->ReplaceCallSites(process_batch_fn, false,
        equals_fn, "Equals", &replaced);
    DCHECK_EQ(replaced, 1);
  }

  process_batch_fn = codegen->ReplaceCallSites(process_batch_fn, false,
      update_tuple_fn, "UpdateAggTuple", &replaced);
  DCHECK_EQ(replaced, 1) << "One call site should be replaced.";
  DCHECK(process_batch_fn != NULL);

  return codegen->OptimizeFunctionWithExprs(process_batch_fn);
}

}

ext: cc

max_stars_repo_name: mapr/impala

max_issues_repo_name: mapr/impala

max_forks_repo_name: mapr/impala