/usr/share/nickle/examples/is-prime.5c is in nickle 2.71-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 | /*
* 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();
|