/*
 * implementation of
 * Agrawal, Kayal, and Saxena's
 * deterministic polytime primality test
 *
 * Copyright © 2002  Bart Massey.
 * All Rights Reserved.  See the file COPYING in this directory
 * for licensing information.
 *
 * Bart Massey 2002/8/16
 */


/* The algorithm (from page 4 of the original abstract)
*
* Input: integer n > 1
* 1.  if ( n is of the form a ** b , b > 1 ) output COMPOSITE;
* 2.  r = 2;
* 3.  while(r < n) {
* 4.      if ( gcd(n,r) =/= 1 ) output COMPOSITE;
* 5.      if (r is prime)
* 6.          let q be the largest prime factor of r - 1;
* 7.          if (q >= 4 sqrt(r) log(n)) and ((n ** ((r-1)/q)) =/= 1) (mod r))
* 8.              break;
* 9.      r = r + 1;
* 10. }
* 11. for a = 1 to 2 sqrt(r) log(n)
* 12.     if ( (x - a) ** n =/= x**n - a (mod x**r - 1, n) ) output COMPOSITE;
* 13. output PRIME;
*/



load "numbers.5c";
import Numbers;

load "polynomial.5c";
import Polynomial;

/* Returns the largest factor of n, computed
   the dumbest way. */
int max_factor(int n) {
    int factor = 1;
    for (int i = 2; i * i <= n; i++) {
	if (n % i == 0) {
	    factor = i;
	    n /= factor;
	}
    }
    return factor;
}

/*
 * returns true iff n is of the form
 * a ** b for some integers a > 1, b > 1.
 */
bool is_xn(int n) {
    /* Could just do prime powers, but this is about as cheap
       and easier. */
    int b = 2;
    real a = exp(log(n) / b);
    while (a  > 1.9) {
	if (floor(a) ** b == n || ceil(a) ** b == n)
	    return true;
	a = exp(log(n) / ++b);
    }
    return false;
}


/*
 * The polytime algorithm below is painfully slow
 * in the small case.  It's only a constant, but
 * a big-deal one.
 */
bool is_small_prime(int n) {
    if (n < 2)
	abort("small prime test below 2");
    global int small = 50;
    if (n >= small)
	return false;
    global bool[*] init_sieve() {
	bool[small] sieve = {true ...};
	int k = 2;
	while (k < small) {
	    for (int i = k + k; i < small; i += k)
		sieve[i] = false;
	    k++;
	    while(k < small && !sieve[k])
		k++;
	}
	return sieve;
    }
    global bool[small] sieve = init_sieve();
    return sieve[n];
}


bool is_prime(int n) {
    if (n < 2)
	abort("prime test below 2");
    if (is_small_prime(n))
	return true;
    if (is_xn(n))
	return false;
    int r;
    for (r = 2; r < n; r++) {
	if (gcd(n, r) != 1)
	    return false;
	if (!is_prime(r))
	    continue;
	int q = max_factor(r - 1);
	if (q ** 2 >= 16 * r * (log2(n) ** 2) &&
	    bigpowmod(n, ((r - 1) / q), r) != 1)
	    break;
    }
    int[r + 1] m1 = {-1, 0...};
    m1[r] = 1;
    polynomial mn = m1;
    mn[0] = -n;
    for (int a = 1; a ** 2 <= 4 * r * (log2(n) ** 2); a++) {
	bool eqn(polynomial m) {
	    polynomial lhs = power_mod((polynomial){-a, 1}, n, m);
	    int[n + 1] rhsp = {-a, 0...};
	    rhsp[n] = 1;
	    polynomial rhs = div_mod(rhsp, m).rem;
	    return lhs == rhs;
	}
	if (eqn(m1) && eqn(mn))
	    return false;
    }
    return true;
}

autoload MillerRabin;

void test_is_prime() {
    for (int i = 2; i < 1000; i++) {
	if (is_prime(i) == MillerRabin::composite(i, 24)) {
	    string msg;
	    if (MillerRabin::composite(i, 24))
		msg = sprintf("%d: false prime\n");
	    else
		msg = sprintf("%d: false composite\n");
	    abort(msg);
	}
	if (true || i % 100 == 0)
	    printf("tested to %d\n", i);
    }
}

test_is_prime();