/* Finds Mersenne primes using the Lucas-Lehmer test
*
* Tom St Denis, tomstdenis@iahu.ca
*/
#include <time.h>
#include <bn.h>
int is_mersenne(long s, int *pp)
{
mp_int n, u, mu;
int res, k;
long ss;
*pp = 0;
if ((res = mp_init(&n)) != MP_OKAY) {
return res;
}
if ((res = mp_init(&u)) != MP_OKAY) {
goto __N;
}
if ((res = mp_init(&mu)) != MP_OKAY) {
goto __U;
}
/* n = 2^s - 1 */
mp_set(&n, 1);
ss = s;
while (ss--) {
if ((res = mp_mul_2(&n, &n)) != MP_OKAY) {
goto __MU;
}
}
if ((res = mp_sub_d(&n, 1, &n)) != MP_OKAY) {
goto __MU;
}
/* setup mu */
if ((res = mp_reduce_setup(&mu, &n)) != MP_OKAY) {
goto __MU;
}
/* set u=4 */
mp_set(&u, 4);
/* for k=1 to s-2 do */
for (k = 1; k <= s - 2; k++) {
/* u = u^2 - 2 mod n */
if ((res = mp_sqr(&u, &u)) != MP_OKAY) {
goto __MU;
}
if ((res = mp_sub_d(&u, 2, &u)) != MP_OKAY) {
goto __MU;
}
/* make sure u is positive */
if (u.sign == MP_NEG) {
if ((res = mp_add(&u, &n, &u)) != MP_OKAY) {
goto __MU;
}
}
/* reduce */
if ((res = mp_reduce(&u, &n, &mu)) != MP_OKAY) {
goto __MU;
}
}
/* if u == 0 then its prime */
if (mp_iszero(&u) == 1) {
*pp = 1;
}
res = MP_OKAY;
__MU: mp_clear(&mu);
__U: mp_clear(&u);
__N: mp_clear(&n);
return res;
}
/* square root of a long < 65536 */
long i_sqrt(long x)
{
long x1, x2;
x2 = 16;
do {
x1 = x2;
x2 = x1 - ((x1 * x1) - x)/(2*x1);
} while (x1 != x2);
if (x1*x1 > x) {
--x1;
}
return x1;
}
/* is the long prime by brute force */
int isprime(long k)
{
long y, z;
y = i_sqrt(k);
for (z = 2; z <= y; z++) {
if ((k % z) == 0) return 0;
}
return 1;
}
int main(void)
{
int pp;
long k;
clock_t tt;
k = 3;
for (;;) {
/* start time */
tt = clock();
/* test if 2^k - 1 is prime */
if (is_mersenne(k, &pp) != MP_OKAY) {
printf("Whoa error\n");
return -1;
}
if (pp == 1) {
/* count time */
tt = clock() - tt;
/* display if prime */
printf("2^%-5ld - 1 is prime, test took %ld ticks\n", k, tt);
}
/* goto next odd exponent */
k += 2;
/* but make sure its prime */
while (isprime(k) == 0) {
k += 2;
}
}
return 0;
}
