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.

`def prime_sieve(): # postponed sieve, by Will Ness yield 2; yield 3; yield 5; yield 7; # original code David Eppstein, 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: yield c # a prime continue 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_prime): """ Calculates n % 2^prime-1 where mersenne_prime=2**prime-1 """ while n > mersenne_prime: n = (n & mersenne_prime) + (n >> prime) return n if n != mersenne_prime else 0 def pow_mod_mersenne(base, exp, prime, mersenne_prime): """ 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_prime) exp >>= 1 base = mod_mersenne(base * base, prime, mersenne_prime) return number def pminus1(prime, mersenne_prime, 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) primes = prime_sieve() for p in primes: if p > B1: break M *= p**int(log_B1/log(p)) M = pow_mod_mersenne(3, 2*M*prime, prime, mersenne_prime) g = gcd(mersenne_prime, M-1) if 1 < g < mersenne_prime: return True if g == mersenne_prime: return False # Start part 2. cache = {0:M} S = mod_mersenne(M*M, prime, mersenne_prime) for d in range(2, int(log(B2)**2), 2): cache[d] = mod_mersenne(cache[d-2] * S, prime, mersenne_prime) HQ = M for k, q in enumerate(primes): if q > B2: break d = q - p if d not in cache: raise HQ = mod_mersenne(HQ*cache[d], prime, mersenne_prime) M = mod_mersenne(M*(HQ-1), prime, mersenne_prime) p = q if k % 100 == 0: g = gcd(M, mersenne_prime) if 1 < g < mersenne_prime: return True if g == mersenne_prime: return False if 1 < gcd(M, mersenne_prime) < mersenne_prime: return True return False `

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