Faster bitwise modulus for Lucas-Lehmer primality test
Question: The [Lucas-Lehmer primality
test](http://en.wikipedia.org/wiki/Lucas%E2%80%93Lehmer_primality_test) tests
prime numbers to determine whether they are also [Mersenne
primes](http://en.wikipedia.org/wiki/Mersenne_prime). One of the bottlenecks
is the modulus operation in the calculation of `(s**2 − 2) % (2**p - 1)`.

Using bitwise operations can speed things up considerably (see the L-L link),
the best I have so far being:

    
    
    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 is where `p` has 5-9 digits and `s` has 10,000+ digits (or
more; not important what they are). Solutions can be tested by `mod((s**2 -
2), p) == (s**2 - 2) % (2**p -1)`. Keep in mind that p - 2 iterations of this
modulus operation are required in the L-L test, each with exponentially
increasing `s`, hence the need for optimization.

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


Answer: In the case where n is much longer than 2^p, you can avoid some quadratic-time
pain 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 screwed up the formatting before; thanks to Benjamin for
pointing this out. Moral: don't copy-and-paste from an Idle window into SO.
Sorry!]

(Note: the criterion for halving the length of n rather than taking p off it,
and the exact choice of k1, are both a bit wrong, but it doesn't matter so I
haven't bothered fixing them.)

If I take p=12345 and n=9**200000 (yes, I know p then isn't prime, but that
doesn't matter here) then this is about 13 times faster.

**Unfortunately this will not help you** , because in the L-L test n is never
bigger than about (2^p)^2. Sorry.

