143 lines
1.9 KiB
C++
143 lines
1.9 KiB
C++
// Bignum prime test (returns 1 if prime, 0 if not)
|
|
|
|
// Uses Algorithm P (probabilistic primality test) from p. 395 of
|
|
// "The Art of Computer Programming, Volume 2" by Donald E. Knuth.
|
|
|
|
#include "stdafx.h"
|
|
#include "defs.h"
|
|
|
|
static int mprimef(unsigned int *, unsigned int *, int);
|
|
|
|
int
|
|
mprime(unsigned int *n)
|
|
{
|
|
int i, k;
|
|
unsigned int *q;
|
|
|
|
// 1?
|
|
|
|
if (MLENGTH(n) == 1 && n[0] == 1)
|
|
return 0;
|
|
|
|
// 2?
|
|
|
|
if (MLENGTH(n) == 1 && n[0] == 2)
|
|
return 1;
|
|
|
|
// even?
|
|
|
|
if ((n[0] & 1) == 0)
|
|
return 0;
|
|
|
|
// n = 1 + (2 ^ k) q
|
|
|
|
q = mcopy(n);
|
|
|
|
k = 0;
|
|
do {
|
|
mshiftright(q);
|
|
k++;
|
|
} while ((q[0] & 1) == 0);
|
|
|
|
// try 25 times
|
|
|
|
for (i = 0; i < 25; i++)
|
|
if (mprimef(n, q, k) == 0)
|
|
break;
|
|
|
|
mfree(q);
|
|
|
|
if (i < 25)
|
|
return 0;
|
|
else
|
|
return 1;
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
//
|
|
// This is the actual implementation of Algorithm P.
|
|
//
|
|
// Input: n The number in question.
|
|
//
|
|
// q n = 1 + (2 ^ k) q
|
|
//
|
|
// k
|
|
//
|
|
// Output: 1 when n is probably prime
|
|
//
|
|
// 0 when n is definitely not prime
|
|
//
|
|
//-----------------------------------------------------------------------------
|
|
|
|
static int
|
|
mprimef(unsigned int *n, unsigned int *q, int k)
|
|
{
|
|
int i, j;
|
|
unsigned int *t, *x, *y;
|
|
|
|
// generate x
|
|
|
|
t = mcopy(n);
|
|
|
|
while (1) {
|
|
for (i = 0; i < MLENGTH(t); i++)
|
|
t[i] = rand();
|
|
x = mmod(t, n);
|
|
if (!MZERO(x) && !MEQUAL(x, 1))
|
|
break;
|
|
mfree(x);
|
|
}
|
|
|
|
mfree(t);
|
|
|
|
// exponentiate
|
|
|
|
y = mmodpow(x, q, n);
|
|
|
|
// done?
|
|
|
|
if (MEQUAL(y, 1)) {
|
|
mfree(x);
|
|
mfree(y);
|
|
return 1;
|
|
}
|
|
|
|
j = 0;
|
|
|
|
while (1) {
|
|
|
|
// y = n - 1?
|
|
|
|
t = msub(n, y);
|
|
|
|
if (MEQUAL(t, 1)) {
|
|
mfree(t);
|
|
mfree(x);
|
|
mfree(y);
|
|
return 1;
|
|
}
|
|
|
|
mfree(t);
|
|
|
|
if (++j == k) {
|
|
mfree(x);
|
|
mfree(y);
|
|
return 0;
|
|
}
|
|
|
|
// y = (y ^ 2) mod n
|
|
|
|
t = mmul(y, y);
|
|
mfree(y);
|
|
y = mmod(t, n);
|
|
mfree(t);
|
|
|
|
// y = 1?
|
|
|
|
if (MEQUAL(y, 1)) {
|
|
mfree(x);
|
|
mfree(y);
|
|
return 0;
|
|
}
|
|
}
|
|
} |