Created
November 20, 2017 07:45
-
-
Save w32zhong/db1c35fded0f9f4188a8b4211a0ade57 to your computer and use it in GitHub Desktop.
primality_test
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#include <stdio.h> | |
#include <stdlib.h> | |
#include <sys/time.h> | |
void rand_init() | |
{ | |
unsigned ticks; | |
struct timeval tv; | |
gettimeofday(&tv, NULL); | |
ticks = tv.tv_sec + tv.tv_usec; | |
srand(ticks); | |
} | |
int pow_mod(int base, int power, int modulo) | |
{ | |
int result = 1; | |
while (power > 0) { | |
if ((power & 1) == 1) { | |
result *= base; | |
result %= modulo; | |
} | |
power >>= 1; | |
base *= base; | |
base %= modulo; | |
} | |
return result; | |
} | |
int mul_reminder(int even, int *k) | |
{ | |
int cnt = 0, r = even; | |
while (r % 2 == 0) { | |
r = r / 2; | |
cnt ++; | |
} | |
*k = cnt; | |
return r; | |
} | |
int miller_rabin(int a, int n) | |
{ | |
int j, k; | |
int r = mul_reminder(n - 1, &k); | |
int z = pow_mod(a, r, n); | |
//printf("choose random base: %d (r=%d, k=%d) \n", a, r, k); | |
if (z == 1 || z == n - 1) | |
return 1; /* likely a prime */ | |
/* repeat k - 1 times */ | |
for (j = 0; j < k - 1; j++) { | |
z = pow_mod(z, 2, n); | |
if (z == 1) | |
return 0; /* not prime witout -1 prior to 1 */ | |
else if (z == n - 1) | |
return 1; /* likely a prime */ | |
} | |
return 0; /* not prime because no -1 prior to a^{2^k·r} */ | |
} | |
int primality_test(int n) | |
{ | |
/* conner cases for very small n */ | |
switch (n) { | |
case 0: | |
case 1: | |
return 0; | |
case 2: | |
case 3: | |
return 1; | |
default: | |
/* pass */; | |
} | |
/* make sure n is odd */ | |
if (n % 2 == 0) | |
return 0; | |
/* pick a random base value in [2, n - 2] */ | |
int base = 2 + rand() % (n - 3); | |
return miller_rabin(base, n); | |
} | |
int test_primes[] = {2, 3, 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, 947, 953, 967, 971, 977, 983, 991, 997}; | |
int true_primality(int n) | |
{ | |
int i; | |
for (i = 0; i < sizeof(test_primes) / sizeof(int); i++) | |
if (test_primes[i] == n) | |
return 1; | |
return 0; | |
} | |
int main() | |
{ | |
int n; | |
rand_init(); | |
for (n = 0; n <= 997; n++) { | |
if (true_primality(n)) { | |
if (!primality_test(n)) | |
printf("judge wrong: %d (which is prime)\n", n); | |
} else { | |
if (primality_test(n)) | |
printf("guess wrong: %d (which is not prime)\n", n); | |
} | |
} | |
return 0; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment