tommath/etc/mersenne.c

145 lines
2.2 KiB
C
Raw Normal View History

2003-02-28 11:06:56 -05:00
/* Finds Mersenne primes using the Lucas-Lehmer test
*
2006-04-06 15:49:59 -04:00
* Tom St Denis, tomstdenis@gmail.com
2003-02-28 11:06:56 -05:00
*/
#include <time.h>
2003-03-12 21:11:11 -05:00
#include <tommath.h>
2003-02-28 11:06:56 -05:00
2003-02-28 11:08:34 -05:00
int
is_mersenne (long s, int *pp)
2003-02-28 11:06:56 -05:00
{
2003-05-29 09:35:26 -04:00
mp_int n, u;
2003-03-12 21:11:11 -05:00
int res, k;
2003-05-29 09:35:26 -04:00
2003-02-28 11:06:56 -05:00
*pp = 0;
2003-02-28 11:08:34 -05:00
if ((res = mp_init (&n)) != MP_OKAY) {
return res;
2003-02-28 11:06:56 -05:00
}
2003-02-28 11:08:34 -05:00
if ((res = mp_init (&u)) != MP_OKAY) {
2004-12-22 21:40:37 -05:00
goto LBL_N;
2003-02-28 11:06:56 -05:00
}
2003-02-28 11:08:34 -05:00
2003-02-28 11:06:56 -05:00
/* n = 2^s - 1 */
2003-05-29 09:35:26 -04:00
if ((res = mp_2expt(&n, s)) != MP_OKAY) {
2004-12-22 21:40:37 -05:00
goto LBL_MU;
2003-02-28 11:06:56 -05:00
}
2003-02-28 11:08:34 -05:00
if ((res = mp_sub_d (&n, 1, &n)) != MP_OKAY) {
2004-12-22 21:40:37 -05:00
goto LBL_MU;
2003-02-28 11:06:56 -05:00
}
2003-02-28 11:08:34 -05:00
2003-02-28 11:06:56 -05:00
/* set u=4 */
2003-02-28 11:08:34 -05:00
mp_set (&u, 4);
2003-02-28 11:06:56 -05:00
/* for k=1 to s-2 do */
for (k = 1; k <= s - 2; k++) {
2003-02-28 11:08:34 -05:00
/* u = u^2 - 2 mod n */
if ((res = mp_sqr (&u, &u)) != MP_OKAY) {
2004-12-22 21:40:37 -05:00
goto LBL_MU;
2003-02-28 11:08:34 -05:00
}
if ((res = mp_sub_d (&u, 2, &u)) != MP_OKAY) {
2004-12-22 21:40:37 -05:00
goto LBL_MU;
2003-02-28 11:08:34 -05:00
}
/* make sure u is positive */
2003-05-29 09:35:26 -04:00
while (u.sign == MP_NEG) {
2003-02-28 11:08:34 -05:00
if ((res = mp_add (&u, &n, &u)) != MP_OKAY) {
2004-12-22 21:40:37 -05:00
goto LBL_MU;
2003-02-28 11:06:56 -05:00
}
2003-02-28 11:08:34 -05:00
}
/* reduce */
2003-05-29 09:35:26 -04:00
if ((res = mp_reduce_2k (&u, &n, 1)) != MP_OKAY) {
2004-12-22 21:40:37 -05:00
goto LBL_MU;
2003-02-28 11:08:34 -05:00
}
2003-02-28 11:06:56 -05:00
}
2003-02-28 11:08:34 -05:00
2003-02-28 11:06:56 -05:00
/* if u == 0 then its prime */
2003-02-28 11:08:34 -05:00
if (mp_iszero (&u) == 1) {
2003-06-06 15:35:48 -04:00
mp_prime_is_prime(&n, 8, pp);
2003-05-29 09:35:26 -04:00
if (*pp != 1) printf("FAILURE\n");
2003-02-28 11:06:56 -05:00
}
2003-02-28 11:08:34 -05:00
2003-02-28 11:06:56 -05:00
res = MP_OKAY;
2004-12-22 21:40:37 -05:00
LBL_MU:mp_clear (&u);
LBL_N:mp_clear (&n);
2003-02-28 11:06:56 -05:00
return res;
}
/* square root of a long < 65536 */
2003-02-28 11:08:34 -05:00
long
i_sqrt (long x)
2003-02-28 11:06:56 -05:00
{
2003-03-12 21:11:11 -05:00
long x1, x2;
2003-02-28 11:08:34 -05:00
x2 = 16;
do {
x1 = x2;
x2 = x1 - ((x1 * x1) - x) / (2 * x1);
} while (x1 != x2);
if (x1 * x1 > x) {
--x1;
}
return x1;
2003-02-28 11:06:56 -05:00
}
/* is the long prime by brute force */
2003-02-28 11:08:34 -05:00
int
isprime (long k)
2003-02-28 11:06:56 -05:00
{
2003-03-12 21:11:11 -05:00
long y, z;
2003-02-28 11:08:34 -05:00
y = i_sqrt (k);
for (z = 2; z <= y; z++) {
if ((k % z) == 0)
return 0;
}
return 1;
}
int
main (void)
2003-02-28 11:06:56 -05:00
{
2003-03-12 21:11:11 -05:00
int pp;
long k;
clock_t tt;
2003-02-28 11:08:34 -05:00
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) {
2003-02-28 11:06:56 -05:00
k += 2;
2003-02-28 11:08:34 -05:00
}
}
return 0;
2003-02-28 11:06:56 -05:00
}
2005-08-01 12:37:28 -04:00
/* $Source$ */
/* $Revision$ */
/* $Date$ */