Just for the fun, I submit the following for consideration, which (btw) beats Steven's test:
Code:
#include <stdbool.h>
#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>
typedef unsigned long prime_t;
typedef bool (*PrimalityTest)(prime_t);
bool PrimalityTest0 ( prime_t i );
static PrimalityTest isPrime = PrimalityTest0;
static const int PRIMES[] = { 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 047, 053, 067, 971, 977, 983, 991, 997 };
prime_t ModExp ( prime_t radix, prime_t exponent, prime_t modulus )
{
prime_t result = 1;
while (exponent > 0)
{
if ((exponent & 1) == 1)
{
result = (result * radix) % modulus;
}
exponent >>= 1;
radix = (radix * radix) % modulus;
}
return result;
}
void FP2 ( prime_t* s, prime_t* d, prime_t n )
{
const prime_t ONE = 1;
prime_t idx2;
n -= 1;
idx2 = 0;
while (!(n & (ONE << idx2)))
{
idx2++;
}
*s = idx2;
*d = n >> idx2;
}
bool RandomisedMillerRabin ( prime_t n )
{
int i;
prime_t a, x, s, d, r;
FP2(&s, &d, n);
for (i = 0; i < 6; i++)
{
a = rand() % (n - 2);
a += 2;
x = ModExp(a, d, n);
if (x == 1 || x == n-1)
continue;
for (r = 1; r < s-1; r++)
{
x = (x*x)%n;
if (x == 1) return false;
if (x == n-1) continue;
}
return false;
}
return true;
}
bool MillerRabin1373653 ( prime_t n )
{
prime_t a, x, s, d, r;
FP2(&s, &d, n);
a = 2;
x = ModExp(a, d, n);
if (x == 1 || x == n-1)
goto test3;
for (r = 1; r < s-1; r++)
{
x = (x*x)%n;
if (x == 1) return false;
if (x == n-1) goto test3;
}
test3:
a = 3;
x = ModExp(a, d, n);
if (x == 1 || x == n-1)
return false;
for (r = 1; r < s-1; r++)
{
x = (x*x)%n;
if (x == 1) return false;
if (x == n-1) break;
}
return true;
}
bool MillerRabin9080191 ( prime_t n )
{
prime_t a, x, s, d, r;
FP2(&s, &d, n);
a = 31;
x = ModExp(a, d, n);
if (x == 1 || x == n-1)
goto test73;
for (r = 1; r < s-1; r++)
{
x = (x*x)%n;
if (x == 1) return false;
if (x == n-1) goto test73;
}
test73:
a = 73;
x = ModExp(a, d, n);
if (x == 1 || x == n-1)
return false;
for (r = 1; r < s-1; r++)
{
x = (x*x)%n;
if (x == 1) return false;
if (x == n-1) break;
}
return true;
}
bool MillerRabin ( prime_t n )
{
if (n < 1373653)
return MillerRabin1373653(n);
else if (n < 9080191)
return MillerRabin9080191(n);
else
return RandomisedMillerRabin(n);
}
bool TableTest ( prime_t n )
{
const int N = sizeof(PRIMES) / sizeof(*PRIMES);
int l = 0, r = N;
int avg, prime;
int i;
for (i = 0; i < 8; i++)
{
avg = (l + r) / 2;
prime = PRIMES[avg];
if (prime == n)
return true;
else if (prime > n)
r = avg;
else
l = avg;
if (r - l <= 1)
return false;
}
return false;
}
bool ActualTest ( prime_t i )
{
if (i < 1000)
return TableTest(i);
else
return MillerRabin(i);
}
bool PrimalityTest0M6 ( prime_t i );
bool PrimalityTest5M6 ( prime_t i )
{
isPrime = PrimalityTest0M6;
return ActualTest(i);
}
bool PrimalityTest4M6 ( prime_t i )
{
isPrime = PrimalityTest5M6;
return false;
}
bool PrimalityTest3M6 ( prime_t i )
{
isPrime = PrimalityTest4M6;
return false;
}
bool PrimalityTest2M6 ( prime_t i )
{
isPrime = PrimalityTest3M6;
return false;
}
bool PrimalityTest1M6 ( prime_t i )
{
isPrime = PrimalityTest2M6;
return ActualTest(i);
}
bool PrimalityTest0M6 ( prime_t i )
{
isPrime = PrimalityTest1M6;
return false;
}
bool PrimalityTest3 ( prime_t i )
{
isPrime = PrimalityTest4M6;
return true;
}
bool PrimalityTest2 ( prime_t i )
{
isPrime = PrimalityTest3;
return true;
}
bool PrimalityTest1 ( prime_t i )
{
isPrime = PrimalityTest2;
return false;
}
bool PrimalityTest0 ( prime_t i )
{
isPrime = PrimalityTest1;
return false;
}
prime_t genprime ( prime_t max )
{
prime_t count = 0, i;
isPrime = PrimalityTest0;
for (i = 0; count < max; i++)
{
if (isPrime(i))
{
count++;
}
}
return count;
}
int main()
{
prime_t begin = 12500, end = 100000;
struct timeval start, finish;
for (prime_t now = begin; now <= end; now += begin) {
printf("Discovering %10lu primes... ", now);
gettimeofday(&start, NULL);
genprime(now);
gettimeofday(&finish, NULL);
double seconds = (double)(finish.tv_sec - start.tv_sec) +
((double)(finish.tv_usec) - (double)(start.tv_usec)) / 1000000.0;
printf("%0.3lfs\n", seconds);
}
return 0;
}
For reference, these are both compiled on a MacBook Pro 15-inch (2008) using clang -arch x86_64 -O3 on Snow Leopard:
Code:
Imladris:tmp alynn$ clang -O3 -arch x86_64 -o StevenPrime StevenPrime.c
Imladris:tmp alynn$ ./StevenPrime
Discovering 12500 primes... 0.019s
Discovering 25000 primes... 0.049s
Discovering 37500 primes... 0.091s
Discovering 50000 primes... 0.139s
Discovering 62500 primes... 0.196s
Discovering 75000 primes... 0.258s
Discovering 87500 primes... 0.325s
Discovering 100000 primes... 0.397s
Imladris:tmp alynn$ clang -O3 -arch x86_64 -o MRPrime MRPrime.c
Imladris:tmp alynn$ ./MRPrime
Discovering 12500 primes... 0.009s
Discovering 25000 primes... 0.017s
Discovering 37500 primes... 0.027s
Discovering 50000 primes... 0.037s
Discovering 62500 primes... 0.046s
Discovering 75000 primes... 0.057s
Discovering 87500 primes... 0.067s
Discovering 100000 primes... 0.076s
