#include "sys.h"
#include "debug.h"
#include <cassert>
#include <libecc/polynomial.h>
#include "primes.h"

using std::cout;
using std::flush;
using std::endl;

unsigned int const m = LIBECC_M;

unsigned int gcd(unsigned int x, unsigned int y)
{
  if (y <= x)
  {
    if (y == 0)
    {
      assert(x);
      return x;
    }
    x = x % y;
  }
  for(;;)
  {
    if (x == 0)
      return y;
    y %= x;
    if (y == 0)
      return x;
    x %= y;
  }
}

// Choose the size of the sieve such that prime_sieve will include all primes <= sqrt(2^(m+1)).
// sieve_size = sqrt(2^(m+1))/2 (+ 1 for rounding up).
int const sieve_size = (int)((1 << ((m + 1) / 2 - 1)) * ((m & 1) ? 1 : 1.414213562373096) + 1);
int prime_sieve[sieve_size];			// Corresponding prime == 2 * index + 3

int const largest_prime = 2 * (sieve_size - 1) + 3;		// Actually this is larger than exp(log_largest_prime).
// Use a macro because using a 'double const' shouldn't compile according to the C++ standard(?!?).
// See http://gcc.gnu.org/bugzilla/show_bug.cgi?id=29905
#define log_largest_prime ((m + 1) * 0.34657359027997265471)	// (m + 1) * ln(2)/2
// See http://primes.utm.edu/howmany.shtml#1, + 1 for rounding up.
int const max_number_of_primes = (int)((largest_prime / log_largest_prime) * (1 + 1.2762 / log_largest_prime) + 1);
unsigned int primes[max_number_of_primes];
int number_of_primes;

void generate_primes(void)
{
  cout << "Generating about " << max_number_of_primes << " primes... " << flush;

  primes[0] = 2;
  number_of_primes = 1;
  int last_index = -1;
  while(number_of_primes < max_number_of_primes)
  {
    do
    {
      ++last_index;
      if (last_index == sieve_size)
      {
	cout << "done (" << number_of_primes << " primes generated)" << endl;
        return;
      }
    }
    while(prime_sieve[last_index]);
    int prime = last_index * 2 + 3;
    primes[number_of_primes++] = prime;
    for(int i = last_index + prime; i < sieve_size; i += prime)
      prime_sieve[i] = 1;
  }

  cout << "done" << endl;
}

void factorize(unsigned int n, std::vector<unsigned int>& product)
{
  assert(n <= primes[number_of_primes - 1] * primes[number_of_primes - 1]);
  for (int i = 0; i < number_of_primes; ++i)
  {
    if (n % primes[i] == 0)
    {
      int count = 0;
      do
      {
        n /= primes[i];
	++count;
      }
      while (n % primes[i] == 0);
      product.push_back(primes[i]);
      product.push_back(count);
    }
  }
  if (n > 1)
  {
    product.push_back(n);
    product.push_back(1);
  }
}

// Slow but correct.
mpz_class foc_reference(unsigned int n)
{
  if (n == 1)
    return 2;
  mpz_class result(1);
  result <<= n;
  for (unsigned int k = 1; k < n; ++k)
    if (n % k == 0)
      result -= foc_reference(k);
  return result;
}

unsigned int foc_table[] = {
    0, 2, 2, 6, 12, 30, 54, 126, 240, 504, 990, 2046, 4020, 8190, 16254, 32730,
    65280, 131070, 261576, 524286, 1047540, 2097018, 4192254, 8388606, 16772880,
    33554400, 67100670, 134217216, 268419060, 536870910, 1073708010, 2147483646 };

// Really, we only need a cache for every divisor of m, but well - m is rather small usually.
mpz_class foc_cache[10000 /* m */];

// Frobenius order count.
// foc(n), n > 0,  returns the number of elements x in a field with
// characteristic 2 such that x^(2^n) = x and x^(2^k) != x for any 0 < k < n.
// In other words, n is the number of times one has to apply Frobenius
// in order to get x again.
mpz_class const& foc(unsigned int n)
{
  mpz_class& result(foc_cache[n]);
  if (result == 0)
  {
    if (n < 32)
      result = foc_table[n];
    else
    {
      result = 1;
      result <<= n;				// Calculate 2^n - 2, if n would be prime then this
      result -= 2;				// would be the number of elements having order n.

      for (unsigned int k = 2; k < 32; ++k)
	if (n % k == 0)
	  result -= foc_table[k];

      for (unsigned int k = 32; k < n; ++k)
	if (n % k == 0)
	  result -= foc(k);
    }
  }
  return result;
}