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.]
6 // We work with integers.
7 #include <cln/integer.h>
12 // Checks whether 2^q-1 is prime, q an odd prime.
13 bool mersenne_prime_p (int q)
15 cl_I m = ((cl_I)1 << q) - 1;
18 for (i = 0, L_i = 4; i < q-2; i++)
19 L_i = mod(L_i*L_i - 2, m);
23 // Same thing, but optimized.
24 bool mersenne_prime_p_opt (int q)
26 cl_I m = ((cl_I)1 << q) - 1;
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));
38 // Now we work with modular integers.
39 #include <cln/modinteger.h>
41 // Same thing, but using modular integers.
42 bool mersenne_prime_p_modint (int q)
44 cl_I m = ((cl_I)1 << q) - 1;
45 cl_modint_ring R = find_modint_ring(m); // Z/mZ
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());
53 #include <cln/io.h> // we do I/O
54 #include <cstdlib> // declares exit()
55 #include <cln/timing.h>
57 int main (int argc, char* argv[])
60 cerr << "Usage: lucaslehmer exponent" << endl;
63 int q = atoi(argv[1]);
64 if (!(q >= 2 && ((q % 2)==1))) {
65 cerr << "Usage: lucaslehmer q with q odd prime" << endl;
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 ";
74 cout << "prime" << endl;
76 cout << "composite" << endl;
79 // Computing time on a i486, 33 MHz: