6
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ Have you tried different Python interpreters? \$\endgroup\$ – Dair Dec 25 '17 at 6:08
  • 3
    \$\begingroup\$ How is mpz defined? Could you also add your imports to make this code runnable? Also, how do you use these functions? \$\endgroup\$ – Graipher Dec 25 '17 at 9:58
  • \$\begingroup\$ 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. \$\endgroup\$ – Oscar Smith Dec 25 '17 at 14:53
  • \$\begingroup\$ Once you use gmpy2, pypy is slower. \$\endgroup\$ – Oscar Smith Dec 25 '17 at 15:31
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ – Oscar Smith Dec 25 '17 at 20:02
1
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ yeah, gmpy2 makes it 7x faster once limit gets to 10000 \$\endgroup\$ – Oscar Smith Jan 1 '18 at 16:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.