/*
 * Mic (10/17): slightly modified this to work with BigInt class.
 * Would be easier if BigInt had "int operator%(BigInt,int)",
 * and some "truncating" conversion back to int.
 * BigInt has y*(x/y)+(x%y)=x, even when y is negative, which is good.
*/

/* ----------------------------------------------------------------------
 *
 *		author	:  Mark a. Weiss
 *		adapted	:  Mic Grigni, Sarah E. Chodrow
 *		project	:  CS 253, Spring 1998
 *		created	:  15 October 1997
 *		edited	:  15 February 1998
 *
 *			numlib.C
 *
 *	Templated recursive numerical algorithms on HugeInts.
 *
 *	HugeInt: must have copy constructor, operator=,
 *		conversion from int, *, /, %, ==, and !=
 *
 * ---------------------------------------------------------------------- */

//#include <iostream.h>
#include <stdio.h>

// Compute X^n ( mod P)
//     assumes that P is not zero and power( 0, 0, P) is 0

template <class HugeInt>
HugeInt power( const HugeInt & X, const HugeInt & n , const HugeInt & P)
{
  if (n == 0) return 1;

  HugeInt tmp = power( (X * X) % P, n / 2, P);
    
  if (n % 2 != 0)
    tmp = (tmp * X) % P;

  return tmp;
}


// Cheap RandInt
// Obviously I should have used the Random class
// but some professors may want to cover this at the
// the same time as the rest of Chapter 7, so I avoided
// the reference to it.

extern "C" double drand48();
int rand_int( int low, int high)
{
  return int (drand48() * (high - low + 1) + low);
}

// Prime routine
// If witness does not return 1, n is definitely composite
// Do this by computing a^i ( mod n) and looking for
// non-trivial square roots of 1 along the way
// HugeInt: must have copy constructor, operator=,
//     conversion from int, *, /, -, %, ==, and !=

template <class HugeInt>
HugeInt witness( const HugeInt & a, const HugeInt & i, const HugeInt & n)
{
  if (i == 0) return 1;

  HugeInt x = witness( a, i / 2, n);
  if (x == 0)    // If n is recursively composite, stop
    return 0;

  // n is not prime if we find a non-trivial square root of 1
  HugeInt y  = ( x * x) % n;
  if (y == 1 && x != 1 && x != n - 1)
    return 0;

  if (i % 2 != 0)
    y = ( a * y) % n;

  return y;
}

// Make numTrials call to witness to check if n is prime

template <class HugeInt>
int is_prime( const HugeInt & n)
{
  const int num_trials = 5;

  for (int counter = 0; counter < num_trials; counter++)
    if (witness( rand_int( 2, 65535)%n, n - 1, n) != 1)
      return 0;

  return 1;
}

template <class HugeInt>
HugeInt gcd( const HugeInt & a, const HugeInt & b)
{
  if (b == 0)
    return a;
  else
    return gcd( b, a % b);
}

// Given a and b, assume gcd( a, b) = 1 
// Find x and y such that a x + b y = 1
// HugeInt: must have copy constructor,
//     zero-parameter constructor, operator=,
//     conversion from int, *, /, +, -, %, ==, and >

template <class HugeInt>
void full_gcd( const HugeInt & a, const HugeInt & b, HugeInt & x,
	       HugeInt & y)
{
  HugeInt x1, y1;

  if (b == 0) {
    x = 1;       // If a != 1, there is no Inverse
    y = 0;       // We omit this check
  }
  else {
    full_gcd( b, a % b, x1, y1);
    x = y1;
    y = x1 - ( a / b) * y1;
  }
}

// Solve a x == 1 ( mod n)
// assume that gcd( a, n) = 1

template <class HugeInt>
HugeInt inverse( const HugeInt & a, const HugeInt & n)
{
  HugeInt x, y;

  full_gcd( a, n, x, y);
  return x > 0 ? x : x + n;
}