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 <cl_integer.h>
   8
   9 // Checks whether 2^q-1 is prime, q an odd prime.
  10 bool mersenne_prime_p (int q)
  11 {
  12         cl_I m = ((cl_I)1 << q) - 1;
  13         int i;
  14         cl_I L_i;
  15         for (i = 0, L_i = 4; i < q-2; i++)
  16                 L_i = mod(L_i*L_i - 2, m);
  17         return (L_i==0);
  18 }
  19
  20 // Same thing, but optimized.
  21 bool mersenne_prime_p_opt (int q)
  22 {
  23         cl_I m = ((cl_I)1 << q) - 1;
  24         int i;
  25         cl_I L_i;
  26         for (i = 0, L_i = 4; i < q-2; i++) {
  27                 L_i = square(L_i) - 2;
  28                 L_i = ldb(L_i,cl_byte(q,q)) + ldb(L_i,cl_byte(q,0));
  29                 if (L_i >= m)
  30                         L_i = L_i - m;
  31         }
  32         return (L_i==0);
  33 }
  34
  35 // Now we work with modular integers.
  36 #include <cl_modinteger.h>
  37
  38 // Same thing, but using modular integers.
  39 bool mersenne_prime_p_modint (int q)
  40 {
  41         cl_I m = ((cl_I)1 << q) - 1;
  42         cl_modint_ring R = cl_find_modint_ring(m); // Z/mZ
  43         int i;
  44         cl_MI L_i;
  45         for (i = 0, L_i = R->canonhom(4); i < q-2; i++)
  46                 L_i = R->minus(R->square(L_i),R->canonhom(2));
  47         return R->equal(L_i,R->zero());
  48 }
  49
  50 #include <cl_io.h> // we do I/O
  51 #include <stdlib.h> // declares exit()
  52 #include <cl_timing.h>
  53
  54 int main (int argc, char* argv[])
  55 {
  56         if (!(argc == 2)) {
  57                 fprint(cl_stderr, "Usage: lucaslehmer exponent\n");
  58                 exit(1);
  59         }
  60         int q = atoi(argv[1]);
  61         if (!(q >= 2 && ((q % 2)==1))) {
  62                 fprint(cl_stderr, "Usage: lucaslehmer q  with q odd prime\n");
  63                 exit(1);
  64         }
  65         bool isprime;
  66         { CL_TIMING; isprime = mersenne_prime_p(q); }
  67         { CL_TIMING; isprime = mersenne_prime_p_opt(q); }
  68         { CL_TIMING; isprime = mersenne_prime_p_modint(q); }
  69         fprint(cl_stdout, "2^");
  70         fprintdecimal(cl_stdout, q);
  71         fprint(cl_stdout, "-1 is ");
  72         if (isprime)
  73                 fprint(cl_stdout, "prime");
  74         else
  75                 fprint(cl_stdout, "composite");
  76         fprint(cl_stdout, "\n");
  77 }
  78
  79 // Computing time on a i486, 33 MHz:
  80 //  1279: 2.02 s
  81 //  2281: 8.74 s
  82 // 44497: 14957 s