/cln-1.3.2/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
 9using namespace std;
10using namespace cln;
11
12// Checks whether 2^q-1 is prime, q an odd prime.
13bool 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.
24bool 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.
42bool 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
57int 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