68 lines
1.4 KiB
Text
68 lines
1.4 KiB
Text
seed = 1 /* seed of the random number generator */
|
|
scale = 0
|
|
|
|
/* Random number from 0 to 32767. */
|
|
define rand() {
|
|
/* Cheap formula (from POSIX) for random numbers of low quality. */
|
|
seed = (seed * 1103515245 + 12345) % 4294967296
|
|
return ((seed / 65536) % 32768)
|
|
}
|
|
|
|
/* Random number in range [from, to]. */
|
|
define rangerand(from, to) {
|
|
auto b, h, i, m, n, r
|
|
|
|
m = to - from + 1
|
|
h = length(m) / 2 + 1 /* want h iterations of rand() % 100 */
|
|
b = 100 ^ h % m /* want n >= b */
|
|
while (1) {
|
|
n = 0 /* pick n in range [b, 100 ^ h) */
|
|
for (i = h; i > 0; i--) {
|
|
r = rand()
|
|
while (r < 68) { r = rand(); } /* loop if the modulo bias */
|
|
n = (n * 100) + (r % 100) /* append 2 digits to n */
|
|
}
|
|
if (n >= b) { break; } /* break unless the modulo bias */
|
|
}
|
|
return (from + (n % m))
|
|
}
|
|
|
|
|
|
|
|
/* n is probably prime? */
|
|
define miller_rabin_test(n, k) {
|
|
auto d, r, a, x, s
|
|
|
|
if (n <= 3) { return (1); }
|
|
if ((n % 2) == 0) { return (0); }
|
|
|
|
/* find s and d so that d * 2^s = n - 1 */
|
|
d = n - 1
|
|
s = 0
|
|
while((d % 2) == 0) {
|
|
d /= 2
|
|
s += 1
|
|
}
|
|
|
|
while (k-- > 0) {
|
|
a = rangerand(2, n - 2)
|
|
x = (a ^ d) % n
|
|
if (x != 1) {
|
|
for (r = 0; r < s; r++) {
|
|
if (x == (n - 1)) { break; }
|
|
x = (x * x) % n
|
|
}
|
|
if (x != (n - 1)) {
|
|
return (0)
|
|
}
|
|
}
|
|
}
|
|
return (1)
|
|
}
|
|
|
|
for (i = 1; i < 1000; i++) {
|
|
if (miller_rabin_test(i, 10) == 1) {
|
|
i
|
|
}
|
|
}
|
|
quit
|