Lucas-Lehmer Faster Bitmap for Primitiveness Test

Question

The Lucas-Lehmer conformance test checks primes to determine if they are also Mersenne prema. . One of the bottlenecks is the operation of the module when calculating (s**2 − 2) % (2**p - 1) .

Using bitwise operations can significantly speed up the process (see LL link), the best I have so far:

 def mod(n,p): """ Returns the value of (s**2 - 2) % (2**p -1)""" Mp = (1<<p) - 1 while n.bit_length() > p: # For Python < 2.7 use len(bin(n)) - 2 > p n = (n & Mp) + (n >> p) if n == Mp: return 0 else: return n 

A simple test case: where p has 5-9 digits and s has 10,000 + digits (or more, not what they are). Solutions can be tested on mod((s**2 - 2), p) == (s**2 - 2) % (2**p -1) . Keep in mind that in the LL test, p - 2 iterations of this module operation are required, each of which exponentially increases s , therefore optimization is required.

Is there a way to speed this up using pure Python (including Python 3)? Is there a better way?

Answer 1

The best improvement I could find was to remove Mp = (1<<p) - 1 from the module function as a whole and pre-calculate it in the LL function before starting the iterations of the LL test. Using while n > Mp: instead of while n.bit_length() > p: also saved some time.

Answer 2

In the case where n is much longer than 2 ^ p, you can avoid the pain in quadratic time by doing something like this:

 def mod1(n,p): while n.bit_length() > 3*p: k = n.bit_length() // p k1 = k>>1 k1p = k1*p M = (1<<k1p)-1 n = (n & M) + (n >> k1p) Mp = (1<<p)-1 while n.bit_length() > p: n = (n&Mp) + (n>>p) if n==Mp: return 0 return n 

[EDITED because I messed up the formatting before; thanks to Benjamin for that. Moral: Do not copy or paste from the standby window in SO. Excuse me!]

(Note: the criterion for shortening the length of n, rather than taking p from it and choosing k1 exactly, is a little wrong, but that doesn't matter, so I didn't worry about fixing them.)

If I take p = 12345 and n = 9 ** 200000 (yes, I know that p is then not simple, but it does not matter here), then this is about 13 times faster.

Unfortunately, this will not help you , because in the LL test, n is never more than about (2 ^ p) ^ 2. Unfortunately.