/usr/share/nickle/examples/miller-rabin.5c is in nickle 2.81-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 | /*
* Miller-Rabin probabilistic test for composite numbers
* as described in Corman/Leiserson/Rivest
*
* composite(n,d): returns false for all prime n,
* and true for most composite n
* (failure probability 1/2**d)
* primebits(n,d): returns a random n-bit number k where !composite(k,d)
*
* Copyright © 1999 Bart Massey.
* All Rights Reserved. See the file COPYING in this directory
* for licensing information.
*
* Bart Massey 1999/1
*/
autoload PRNG;
namespace MillerRabin {
int[*] primes = {
2, 3, 5, 7, 11,
13, 17, 19, 23, 29
};
int nprimes = dim(primes);
typedef struct { int pow, wit; } witness_result;
/*
* Modified version of bigpowmod() from numbers.5c.
* Computes core of Miller-Rabin test
* as suggested by Cormen/Leiserson/Rivest.
*/
witness_result witnessexp(int b, int e, int m) {
switch (e) {
case 0:
return (witness_result){ .pow = 0, .wit = 1};
case 1:
return (witness_result){ .pow = b % m, .wit = 0};
}
witness_result tmp = witnessexp(b, e // 2, m);
if (tmp.wit != 0)
return tmp;
int t = (tmp.pow * tmp.pow) % m;
if (t == 1 && tmp.pow != 1 && tmp.pow != m - 1) {
tmp.wit = tmp.pow;
tmp.pow = t;
return tmp;
}
if (e % 2 == 0)
tmp.pow = t;
else
tmp.pow = (t * b) % m;
return tmp;
}
/* Rest of Miller-Rabin test */
bool witness(int a, int n) {
witness_result we = witnessexp(a, n - 1, n);
if (we.wit != 0)
return true;
if (we.pow != 1)
return true;
return false;
}
/* Try small primes, then Miller-Rabin */
public bool composite(int n, int d) {
int i, j;
for (i = 0; i < nprimes && primes[i] < n; i++)
if (n % primes[i] == 0)
return true;
for (j = 0; j < d; j++) {
int a = 1 + PRNG::randint(n - 1);
if (witness(a, n))
return true;
}
return false;
}
/* generate an n-bit prime (with probability 1-(2**-d)) number */
public int primebits(int n, int d) {
while (true) {
int q = PRNG::randbits(n - 1) + 2**(n - 1);
bool why = composite(q, d);
if (!why)
return q;
# printf("*\n");
}
}
}
|