#include <stdio.h>
#include <limits.h>
#include <float.h>
#include <stdlib.h>
#include <math.h>


#include <NTL/version.h>


#if (defined(__GNUC__) && (defined(__i386__) || defined(__i486__) || defined(__i586__)))


#define AutoFix (1)

#else

#define AutoFix (0)

#endif


int val_int(int x);
unsigned int val_uint(unsigned int x);

long val_long(long x);
unsigned long val_ulong(unsigned long x);

size_t val_size_t(size_t x);

double val_double(double x);

void touch_int(int* x);
void touch_uint(unsigned int* x);

void touch_long(long* x);
void touch_ulong(unsigned long* x);

void touch_size_t(size_t* x);

void touch_double(double* x);




double power2(long k)
{
   long i;
   double res;

   res = 1;

   for (i = 1; i <= k; i++)
      res = res * 2;

   return res;
}


long DoubleRounding(long dp)
{
   double a = power2(dp-1) + 1;
   double b = (power2(dp)-1)/power2(dp+1);
   register double x = a + b;
   double y = x;

   touch_double(&y);

   if (y != power2(dp-1) + 1)
      return 1;
   else 
      return 0; 
}



long DoublePrecision(void)
{
   double eps, one, res;
   long k;

   one = val_double(1.0);
   eps = val_double(1.0);

   k = 0;

   do {
      double tmp;

      k++;
      eps *= 1.0/2.0;
      tmp = 1.0 + eps;
      touch_double(&tmp);
      res = tmp - one;
   } while (res == eps);

   return k;
}

long DoublePrecision1(void)
{
   double eps, one, res;
   long k;

   one = val_double(1.0);
   eps = val_double(1.0);

   k = 0;

   do {
      register double tmp;

      k++;
      eps *= 1.0/2.0;
      tmp = 1.0 + eps;
      res = tmp - one;
   } while (res == eps);

   return k;
}


union d_or_rep {
   double d;
   unsigned long rep[2];
};

long RepTest(void)
{
   union d_or_rep v;

   if (sizeof(double) != 2*sizeof(long))
      return 0;

   v.rep[0] = v.rep[1] = 0;

   v.d = 565656565656.0;

   if (v.rep[0] == 0x42607678 && v.rep[1] == 0x46f30000)
      return 1;
   else if (v.rep[1] == 0x42607678 && v.rep[0] == 0x46f30000)
      return -1;
   else
      return 0;
}

void print2k(FILE *f, long k, long bpl)
{
   long m, l;
   long first;

   if (k <= 0) {
      fprintf(f, "((double) 1.0)");
      return;
   }

   m = bpl - 2;
   first = 1;

   fprintf(f, "(");

   while (k > 0) {
      if (k > m)
         l = m;
      else
         l = k;

      k = k - l;

      if (first)
         first = 0;
      else 
         fprintf(f, "*");

      fprintf(f, "((double)(1L<<%ld))", l);
   }

   fprintf(f, ")");
}






void print_mul_body(FILE *f, long n1, long k, long fn,
                    long half_flag, long short_flag)
{

  long n, i, chop, r;
  unsigned long mask, mask2;

  if (half_flag)
    n = n1/2;
  else
    n = n1;

  chop = n % k; /* first block */
  if (chop == 0)
    chop = k;
  r = n - k;
  mask = (1UL << k) - 1UL;

  fprintf(f, "\n\n#define NTL_");
  if (half_flag) fprintf(f, "HALF_");
  if (short_flag) fprintf(f, "SHORT_");

  fprintf(f, "BB_MUL_CODE%ld \\\n", fn);

  
  if (fn > 0) /* Mul1/AddMul1 */
    {
      fprintf(f, "   long i;\\\n");
      fprintf(f, "   _ntl_ulong carry = 0, b;\\\n");
    }
  fprintf(f, "   _ntl_ulong hi, lo, t;\\\n");
  fprintf(f, "   _ntl_ulong A[%ld];\\\n", 1L << k);
  fprintf(f, "   A[0] = 0;\\\n");

  fprintf(f, "   A[1] = a;\\\n");

  for (i = 2; i < (1L << k); i++)
    {
      if (i % 2 == 0)
        fprintf(f, "   A[%ld] = A[%ld] << 1;\\\n", i, i / 2);
      else
        fprintf(f, "   A[%ld] = A[%ld] ^ A[1];\\\n", i, i - 1);
    }

  if (fn > 0)
    {
      fprintf(f, "   for (i = 0; i < sb; i++) {\\\n");
      fprintf(f, "      b = bp[i];\\\n");
      fprintf(f, "   ");
    }

  fprintf(f, "   lo = A[b & %lu]; ", mask);
  fprintf(f, "t = A[(b >> %ld) & %lu]; ", k, mask);
  fprintf(f, "hi = t >> %ld; lo ^= t << %ld;\\\n", n1-k, k);

  for (i = 2*k; i < n - chop; i += k) {
    if (fn > 0) fprintf(f, "   ");
    fprintf(f, "   t = A[(b >> %ld) & %lu]; ", i, mask);
    fprintf(f, "hi ^= t >> %ld; lo ^= t << %ld;\\\n", n1-i, i);
  }

  if (fn > 0) fprintf(f, "   ");
  fprintf(f, "   t = A[b >> %ld]; ", n-chop);
  fprintf(f, "hi ^= t >> %ld; lo ^= t << %ld;\\\n", n1-i, i);

  mask = 0;
  for (i = 0; i < n; i += k)
    mask |= 1UL << i;
  mask = ~mask;
  if (half_flag) mask &= (1UL << n) - 1UL;
  mask2 = mask;
  
  if (!short_flag) {
    for (i = 1; i < k; i++)
      {
	if (fn > 0) fprintf(f, "   ");

        if (i == 1)
	   fprintf(f, "   if (a >> %ld) ", n1-i);
        else
           fprintf(f, "   if ((a >> %ld) & 1) ", n1-i);

	/* bit n1-i from a was not considered in blocks of
           k bits from b for index j >= i */
	fprintf(f, "hi ^= ((b & 0x%lxUL) >> %ld);\\\n", mask2, i);
	mask2 = (mask2 << 1) & mask;
      }
   }

  if (fn > 0) fprintf(f, "   ");

  if (fn == 0)
    {
      fprintf(f, "   c[0] = lo; ");
      fprintf(f, "   c[1] = hi;\\\n");
    }
  else if (fn == 1 || fn == 3)
    {
      fprintf(f, "   cp[i] = carry ^ lo; ");
      fprintf(f, "   carry = hi;\\\n");
    }
  else if (fn == 2)
    {
      fprintf(f, "   cp[i] ^= (carry ^ lo); ");
      fprintf(f, "   carry = hi;\\\n");
    }
  if (fn > 0)
    {

      fprintf(f, "   }\\\n");

      if (fn == 1 || fn == 3) 
        fprintf(f, "   cp[sb] = carry;\\\n");
      else
        fprintf(f, "   cp[sb] ^= carry;\\\n");

    }
  fprintf(f, "\n\n\n");


}





/*
 * This generates anternative code that runs significantly faster 
 * on some machines, like a PowerPC (and probably other RISC machines).
 * It makes it easier for the compiler to schedule instrucyions better,
 * and it avoids branches.  It seems like this does not help
 * on x86 machines (and can even make things worse).
 */


void print_alt_mul_body(FILE *f, long n1, long k, long fn, 
                        long half_flag, long short_flag)
{

  long n, i, chop, r;
  unsigned long mask, mask2;

  if (half_flag)
    n = n1/2;
  else
    n = n1;

  chop = n % k; /* first block */
  if (chop == 0)
    chop = k;
  r = n - k;
  mask = (1UL << k) - 1UL;

  fprintf(f, "\n\n#define NTL_ALT_");
  if (half_flag) fprintf(f, "HALF_");
  if (short_flag) fprintf(f, "SHORT_");

  fprintf(f, "BB_MUL_CODE%ld \\\n", fn);
  
  if (fn > 0) /* Mul1/AddMul1 */
    {
      fprintf(f, "   long i;\\\n");
      fprintf(f, "   _ntl_ulong carry = 0;\\\n");
    }
  fprintf(f, "   _ntl_ulong A[%ld];\\\n", 1L << k);
  fprintf(f, "   A[0] = 0;\\\n");

  fprintf(f, "   A[1] = a;\\\n");

  for (i = 2; i < (1L << k); i++)
    {
      if (i % 2 == 0)
        fprintf(f, "   A[%ld] = A[%ld] << 1;\\\n", i, i / 2);
      else
        fprintf(f, "   A[%ld] = A[%ld] ^ A[1];\\\n", i, i - 1);
    }

  if (fn > 0)
    {
      fprintf(f, "   for (i = 0; i < sb; i++) {\\\n");
      fprintf(f, "      const _ntl_ulong b = bp[i];\\\n");
    }

  for (i = k; i < n - chop; i += k) {
    if (fn > 0) fprintf(f, "   ");
    fprintf(f, "   const _ntl_ulong t%ld = A[(b >> %ld) & %lu]; \\\n", i, i, mask);

  }
  if (fn > 0) fprintf(f, "   ");
  fprintf(f, "   const _ntl_ulong t%ld = A[b >> %ld]; \\\n", n-chop, n-chop);

  if (fn > 0) fprintf(f, "   ");
  fprintf(f, "   const _ntl_ulong lo = A[b & %lu] \\\n", mask);

  for (i = k; i < n; i += k) {
    if (fn > 0) fprintf(f, "   ");
    fprintf(f, "      ^ (t%ld << %ld)", i, i);
    if (i == n - chop)
      fprintf(f, ";\\\n");
    else
      fprintf(f, "\\\n");
  }

  for (i = k; i < n; i += k) {
    if (fn > 0) fprintf(f, "   ");
    if (i == k)
      fprintf(f, "   const _ntl_ulong hi = ");
    else
      fprintf(f, "      ^ ");
    fprintf(f, "(t%ld >> %ld)", i, n1-i);
    if (i == n - chop && short_flag)
      fprintf(f, ";\\\n");
    else
      fprintf(f, "\\\n");
      
      
  }
 
  mask = 0;
  for (i = 0; i < n; i += k)
    mask |= 1UL << i;
  mask = ~mask;
  if (half_flag) mask &= (1UL << n) - 1UL;
  mask2 = mask;
  
  if (!short_flag) {
    for (i = 1; i < k; i++)
      {

        /* bit n1-i from a was not considered in blocks of
           k bits from b for index j >= i */

        if (fn > 0) fprintf(f, "   ");

	
        if (i == 1)
          fprintf(f, 
             "      ^ (((b & 0x%lxUL) >> %ld) & (-(a >> %ld)))", 
        	mask2, i, n1-1);
        else {
          fprintf(f, 
             "      ^ (((b & 0x%lxUL) >> %ld) & (-((a >> %ld) & 1UL)))",  
                  mask2, i, n1-i);
        }
        if (i == k-1) 
           fprintf(f, ";\\\n");
        else
            fprintf(f, "\\\n");


	
	mask2 = (mask2 << 1) & mask;
      }
   }

  if (fn > 0) fprintf(f, "   ");

  if (fn == 0)
    {
      fprintf(f, "   c[0] = lo; ");
      fprintf(f, "   c[1] = hi;\\\n");
    }
  else if (fn == 1)
    {
      fprintf(f, "   cp[i] = carry ^ lo; ");
      fprintf(f, "   carry = hi;\\\n");
    }
  else if (fn == 2)
    {
      fprintf(f, "   cp[i] ^= (carry ^ lo); ");
      fprintf(f, "   carry = hi;\\\n");
    }
  if (fn > 0)
    {

      fprintf(f, "   }\\\n");

      if (fn == 1 || fn == 3) 
        fprintf(f, "   cp[sb] = carry;\\\n");
      else
        fprintf(f, "   cp[sb] ^= carry;\\\n");

    }
  fprintf(f, "\n\n\n");


}




void print_alt1_mul_body(FILE *f, long n1, long k, long fn,
                    long half_flag, long short_flag)
{

  long n, i, chop, r;
  unsigned long mask, mask2;

  if (half_flag)
    n = n1/2;
  else
    n = n1;

  chop = n % k; /* first block */
  if (chop == 0)
    chop = k;
  r = n - k;
  mask = (1UL << k) - 1UL;

  fprintf(f, "\n\n#define NTL_ALT1_");
  if (half_flag) fprintf(f, "HALF_");
  if (short_flag) fprintf(f, "SHORT_");

  fprintf(f, "BB_MUL_CODE%ld \\\n", fn);

  
  if (fn > 0) /* Mul1/AddMul1 */
    {
      fprintf(f, "   long i;\\\n");
      fprintf(f, "   _ntl_ulong carry = 0, b;\\\n");
    }
  fprintf(f, "   _ntl_ulong hi, lo, t;\\\n");
  fprintf(f, "   _ntl_ulong A[%ld];\\\n", 1L << k);
  fprintf(f, "   A[0] = 0;\\\n");

  fprintf(f, "   A[1] = a;\\\n");

  for (i = 2; i < (1L << k); i++)
    {
      if (i % 2 == 0)
        fprintf(f, "   A[%ld] = A[%ld] << 1;\\\n", i, i / 2);
      else
        fprintf(f, "   A[%ld] = A[%ld] ^ A[1];\\\n", i, i - 1);
    }

  if (fn > 0)
    {
      fprintf(f, "   for (i = 0; i < sb; i++) {\\\n");
      fprintf(f, "      b = bp[i];\\\n");
      fprintf(f, "   ");
    }

  fprintf(f, "   lo = A[b & %lu]; ", mask);
  fprintf(f, "t = A[(b >> %ld) & %lu]; ", k, mask);
  fprintf(f, "hi = t >> %ld; lo ^= t << %ld;\\\n", n1-k, k);

  for (i = 2*k; i < n - chop; i += k) {
    if (fn > 0) fprintf(f, "   ");
    fprintf(f, "   t = A[(b >> %ld) & %lu]; ", i, mask);
    fprintf(f, "hi ^= t >> %ld; lo ^= t << %ld;\\\n", n1-i, i);
  }

  if (fn > 0) fprintf(f, "   ");
  fprintf(f, "   t = A[b >> %ld]; ", n-chop);
  fprintf(f, "hi ^= t >> %ld; lo ^= t << %ld;\\\n", n1-i, i);

  mask = 0;
  for (i = 0; i < n; i += k)
    mask |= 1UL << i;
  mask = ~mask;
  if (half_flag) mask &= (1UL << n) - 1UL;
  mask2 = mask;