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, 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, 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

static long ran_initialized = 0;
static _ZZ_arc4_key ran_key;

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

   typedef unsigned char u_char;
   txt = NTL_NEW_OP u_char[nb + 68];
   if (!txt) Error("out of memory");

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

   delete [] txt;
}


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) SetSeed(ZZ::zero());
   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) 
      Error("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) 
      Error("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) 
      Error("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))
      Error("RandomBits: length too big");

   long nb = (l+7)/8;

   static unsigned char *buf = 0;
   static long buf_len = 0;

   if (nb > buf_len) {
      if (buf) delete [] buf;
      buf_len = ((nb + 1023)/1024)*1024; // allocate in 1024-byte lots
      typedef unsigned char u_char;
      buf = NTL_NEW_OP u_char[buf_len];
      if (!buf) Error("out of memory");
   }

   ran_bytes(buf, nb);

   static ZZ 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))
      Error("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;

   static ZZ t, r, 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) {
      static ZZ Bnd, 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)
{
   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) Error("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) Error("GenPrime: bad length");

   if (k > (1L << 20)) Error("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) Error("GenPrime: bad length");

   if (k >= NTL_BITS_PER_LONG) Error("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) Error("GenGermainPrime: bad length");

   if (k > (1L << 20)) Error("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)
      Error("GenGermainPrime_long: length too long");

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


NTL_END_IMPL