# Maximum performance for Pollard's P-1 function

I am trying to find mersenne primes using python, and following the advice from https://www.mersenne.org/various/math.php, I have an intermediate step where I am using Pollard's P-1 factoring algorithm. The faster I can make this step, the higher bounds I can use, which will make my program much faster. Everything other than the prime sieve was written by me, and the sieve is not a bottleneck of the program at all.

from itertools import count, compress, takewhile, cycle
from math import gcd, fsum, log
try:
from gmpy2 import mpz
except ImportError:
mpz = int

def prime_sieve(start=1, end=float('inf')): # postponed sieve, by Will Ness
for c in (2,3,5,7):                     # original code David Eppstein,
if start <= c <= end:               # start and end added by Oscar Smith
yield c
elif end < c:
return
sieve = {}                              # Alex Martelli, ActiveState Recipe 2002
ps = prime_sieve()                      # a separate base Primes Supply:
p = next(ps) and next(ps)               # (3) a Prime to add to dict
q = p*p                                 # (9) its sQuare
for c in count(9,2):                    # the Candidate
if c in sieve:                      # c’s a multiple of some base prime
s = sieve.pop(c)                # i.e. a composite ; or
elif c < q:
if start <= c <= end:
yield c                     # a prime
continue
elif end < c:
return
else:   # (c==q):            # or the next base prime’s square:
s=count(q+2*p,2*p)       #    (9+6, by 6 : 15,21,27,33,...)
p=next(ps)               #    (5)
q=p*p                    #    (25)
for m in s:                  # the next multiple
if m not in sieve:       # no duplicates
break
sieve[m] = s                 # original test entry: ideone.com/WFv4f

def mod_mersenne(n, prime, mersenne):
""" Calculates n % 2^prime-1 where mersenne_prime=2**prime-1 """
while n > mersenne:
n = (n & mersenne) + (n >> prime)
return n if n != mersenne else 0

def pow_mod_mersenne(base, exp, prime, mersenne):
""" Calculates base^exp % 2^prime-1 where mersenne_prime=2**prime-1 """
number = 1
while exp:
if exp & 1:
number = mod_mersenne(number * base, prime, mersenne)
exp >>= 1
base = mod_mersenne(base * base, prime, mersenne)
return number

def p_minus1(prime, mersenne, B1, B2):
""" Does 2**prime-1 have a factor 2*k*prime+1?
such that the largest prime factor of k is less than limit"""
log_B1 = log(B1)
M = mpz(1)
for p in prime_sieve(end=B1):
M *= p**int(log_B1/log(p))
M = pow_mod_mersenne(3, 2*M*prime, prime, mersenne)
g = gcd(mersenne, M-1)
if 1 < g < mersenne:
return True
if g == mersenne:
return False
return p_minus1_stage_2(prime, mersenne, M, p, B2)

def p_minus1_stage_2(prime, mersenne, M, B1, B2):
delta_cache = {0:M}
delta_max = 0
S = mod_mersenne(M*M, prime, mersenne)
p_old = B1
HQ = M
for k, p in enumerate(prime_sieve(start=B1, end=B2)):
delta = p - p_old
if delta not in delta_cache:
for d in range(delta_max, delta, 2):
delta_cache[d+2] = mod_mersenne(delta_cache[d] * S, prime, mersenne)
delta_max = delta
HQ = mod_mersenne(HQ*delta_cache[delta], prime, mersenne)
M = mod_mersenne(M*(HQ-1), prime, mersenne)
p_old = p
if k % 100 == 0:
g = gcd(M, mersenne)
if 1 < g < mersenne:
return True
if g == mersenne:
return False
return 1 < gcd(M, mersenne) < mersenne

if __name__ == '__main__':
primes = prime_sieve()
next(primes)
for prime in primes:
B1 = int(10*prime.bit_length())
p_minus1(prime, 2**mpz(prime)-1, B1, 10*B1)
if prime>2000:
break


Performance improvements would be great. Style is nice too, although definitely not the primary goal here.

• Have you tried different Python interpreters? – Dair Dec 25 '17 at 6:08
• How is mpz defined? Could you also add your imports to make this code runnable? Also, how do you use these functions? – Graipher Dec 25 '17 at 9:58
• Oops,I meant to take out mpz. It's from gmpy2, and is just a faster version of int. I'll add a run script. – Oscar Smith Dec 25 '17 at 14:53
• Once you use gmpy2, pypy is slower. – Oscar Smith Dec 25 '17 at 15:31
• This should fix it. There was a hidden assumption that you wouldn't try to pass B1 less than 3, as otherwise the difference between consecutive primes can be odd. – Oscar Smith Dec 25 '17 at 20:02

This is not much help, but I checked that most of the runtime is spent in integer multiplication: about 90 % if I increase the main loop limit to 10000. (This was with builtin int not gmpy2.) So, unless you can improve the algorithm to do fewer/smaller multiplications, or to perform modular multiplication in a faster way, there is not much to improve.
• yeah, gmpy2 makes it 7x faster once limit gets to 10000 – Oscar Smith Jan 1 '18 at 16:52