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.