examples/lucaslehmer.cc

   1 // Check whether a mersenne number is prime,
   2 // using the Lucas-Lehmer test.
   3 // [Donald Ervin Knuth: The Art of Computer Programming, Vol. II:
   4 //  Seminumerical Algorithms, second edition. Section 4.5.4, p. 391.]
   5
   6 // We work with integers.
   7 #include <cln/integer.h>
   8
   9 using namespace std;
  10 using namespace cln;
  11
  12 // Checks whether 2^q-1 is prime, q an odd prime.
  13 bool mersenne_prime_p (int q)
  14 {
  15         cl_I m = ((cl_I)1 << q) - 1;
  16         int i;
  17         cl_I L_i;
  18         for (i = 0, L_i = 4; i < q-2; i++)
  19                 L_i = mod(L_i*L_i - 2, m);
  20         return (L_i==0);
  21 }
  22
  23 // Same thing, but optimized.
  24 bool mersenne_prime_p_opt (int q)
  25 {
  26         cl_I m = ((cl_I)1 << q) - 1;
  27         int i;
  28         cl_I L_i;
  29         for (i = 0, L_i = 4; i < q-2; i++) {
  30                 L_i = square(L_i) - 2;
  31                 L_i = ldb(L_i,cl_byte(q,q)) + ldb(L_i,cl_byte(q,0));
  32                 if (L_i >= m)
  33                         L_i = L_i - m;
  34         }
  35         return (L_i==0);
  36 }
  37
  38 // Now we work with modular integers.
  39 #include <cln/modinteger.h>
  40
  41 // Same thing, but using modular integers.
  42 bool mersenne_prime_p_modint (int q)
  43 {
  44         cl_I m = ((cl_I)1 << q) - 1;
  45         cl_modint_ring R = find_modint_ring(m); // Z/mZ
  46         int i;
  47         cl_MI L_i;
  48         for (i = 0, L_i = R->canonhom(4); i < q-2; i++)
  49                 L_i = R->minus(R->square(L_i),R->canonhom(2));
  50         return R->equal(L_i,R->zero());
  51 }
  52
  53 #include <cln/io.h> // we do I/O
  54 #include <cstdlib>  // declares exit()
  55 #include <cln/timing.h>
  56
  57 int main (int argc, char* argv[])
  58 {
  59         if (!(argc == 2)) {
  60                 cerr << "Usage: lucaslehmer exponent" << endl;
  61                 exit(1);
  62         }
  63         int q = atoi(argv[1]);
  64         if (!(q >= 2 && ((q % 2)==1))) {
  65                 cerr << "Usage: lucaslehmer q  with q odd prime" << endl;
  66                 exit(1);
  67         }
  68         bool isprime;
  69         { CL_TIMING; isprime = mersenne_prime_p(q); }
  70         { CL_TIMING; isprime = mersenne_prime_p_opt(q); }
  71         { CL_TIMING; isprime = mersenne_prime_p_modint(q); }
  72         cout << "2^" << q << "-1 is ";
  73         if (isprime)
  74                 cout << "prime" << endl;
  75         else
  76                 cout << "composite" << endl;
  77 }
  78
  79 // Computing time on a i486, 33 MHz:
  80 //  1279: 2.02 s
  81 //  2281: 8.74 s
  82 // 44497: 14957 s