1f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project/* Finds Mersenne primes using the Lucas-Lehmer test 2f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project * 3f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project * Tom St Denis, tomstdenis@gmail.com 4f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project */ 5f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project#include <time.h> 6f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project#include <tommath.h> 7f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 8f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Projectint 9f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Projectis_mersenne (long s, int *pp) 10f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project{ 11f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project mp_int n, u; 12f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project int res, k; 13f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 14f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project *pp = 0; 15f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 16f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project if ((res = mp_init (&n)) != MP_OKAY) { 17f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project return res; 18f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } 19f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 20f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project if ((res = mp_init (&u)) != MP_OKAY) { 21f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project goto LBL_N; 22f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } 23f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 24f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project /* n = 2^s - 1 */ 25f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project if ((res = mp_2expt(&n, s)) != MP_OKAY) { 26f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project goto LBL_MU; 27f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } 28f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project if ((res = mp_sub_d (&n, 1, &n)) != MP_OKAY) { 29f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project goto LBL_MU; 30f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } 31f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 32f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project /* set u=4 */ 33f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project mp_set (&u, 4); 34f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 35f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project /* for k=1 to s-2 do */ 36f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project for (k = 1; k <= s - 2; k++) { 37f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project /* u = u^2 - 2 mod n */ 38f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project if ((res = mp_sqr (&u, &u)) != MP_OKAY) { 39f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project goto LBL_MU; 40f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } 41f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project if ((res = mp_sub_d (&u, 2, &u)) != MP_OKAY) { 42f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project goto LBL_MU; 43f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } 44f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 45f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project /* make sure u is positive */ 46f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project while (u.sign == MP_NEG) { 47f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project if ((res = mp_add (&u, &n, &u)) != MP_OKAY) { 48f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project goto LBL_MU; 49f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } 50f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } 51f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 52f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project /* reduce */ 53f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project if ((res = mp_reduce_2k (&u, &n, 1)) != MP_OKAY) { 54f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project goto LBL_MU; 55f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } 56f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } 57f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 58f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project /* if u == 0 then its prime */ 59f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project if (mp_iszero (&u) == 1) { 60f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project mp_prime_is_prime(&n, 8, pp); 61f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project if (*pp != 1) printf("FAILURE\n"); 62f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } 63f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 64f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project res = MP_OKAY; 65f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source ProjectLBL_MU:mp_clear (&u); 66f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source ProjectLBL_N:mp_clear (&n); 67f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project return res; 68f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project} 69f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 70f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project/* square root of a long < 65536 */ 71f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Projectlong 72f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Projecti_sqrt (long x) 73f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project{ 74f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project long x1, x2; 75f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 76f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project x2 = 16; 77f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project do { 78f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project x1 = x2; 79f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project x2 = x1 - ((x1 * x1) - x) / (2 * x1); 80f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } while (x1 != x2); 81f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 82f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project if (x1 * x1 > x) { 83f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project --x1; 84f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } 85f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 86f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project return x1; 87f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project} 88f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 89f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project/* is the long prime by brute force */ 90f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Projectint 91f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Projectisprime (long k) 92f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project{ 93f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project long y, z; 94f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 95f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project y = i_sqrt (k); 96f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project for (z = 2; z <= y; z++) { 97f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project if ((k % z) == 0) 98f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project return 0; 99f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } 100f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project return 1; 101f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project} 102f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 103f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 104f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Projectint 105f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Projectmain (void) 106f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project{ 107f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project int pp; 108f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project long k; 109f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project clock_t tt; 110f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 111f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project k = 3; 112f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 113f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project for (;;) { 114f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project /* start time */ 115f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project tt = clock (); 116f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 117f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project /* test if 2^k - 1 is prime */ 118f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project if (is_mersenne (k, &pp) != MP_OKAY) { 119f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project printf ("Whoa error\n"); 120f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project return -1; 121f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } 122f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 123f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project if (pp == 1) { 124f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project /* count time */ 125f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project tt = clock () - tt; 126f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 127f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project /* display if prime */ 128f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project printf ("2^%-5ld - 1 is prime, test took %ld ticks\n", k, tt); 129f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } 130f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 131f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project /* goto next odd exponent */ 132f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project k += 2; 133f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 134f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project /* but make sure its prime */ 135f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project while (isprime (k) == 0) { 136f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project k += 2; 137f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } 138f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project } 139f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project return 0; 140f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project} 141f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project 142f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project/* $Source: /cvs/libtom/libtommath/etc/mersenne.c,v $ */ 143f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project/* $Revision: 1.3 $ */ 144f7fc46c63fdc8f39234fea409b8dbe116d73ebf8The Android Open Source Project/* $Date: 2006/03/31 14:18:47 $ */ 145