An example of Miller-Rabin primality test breaking

arnault.sage

#!/usr/bin/env sage

# François Arnault. 1995. Constructing Carmichael Numbers which are Strong Pseudoprimes to Several Bases
# https://doi.org/10.1006/jsco.1995.1042


from sys import stderr
from random import choice, getrandbits


def miller_rabin(bases, n):
    if n == 2 or n == 3:
        return True

    if n % 2 == 0:
        return False

    r, s = 0, n - 1
    while s % 2 == 0:
        r += 1
        s //= 2
    
    ZmodN = Zmod(n)

    for b in map(ZmodN, bases):
        x = b ^ s
        if x == 1 or x == -1:
            continue
        for _ in range(r - 1):
            x = x ^ 2
            if x == -1:
                break
        else:
            return False
    
    return True


def search_pseudoprime(bases, coeffs, rlen, iters, verbose=False):
    modules = [4 * b for b in bases]

    residues = dict()

    for b, m in zip(bases, modules):
        residues[b] = set()
        for p in primes(3, 1024 * max(bases)):
            if kronecker(b, p) == -1:
                residues[b].add(p % m)

    sets = dict()

    for b, m in zip(bases, modules):
        s = []
        for c in coeffs:
            s.append({(inverse_mod(c, m) * (r + c - 1)) % m for r in residues[b]})
        sets[b] = list(set.intersection(*s))

    # only support this
    assert len(coeffs) == 3

    coeffs_inv = [
        1, 
        coeffs[1] - inverse_mod(coeffs[2], coeffs[1]), 
        coeffs[2] - inverse_mod(coeffs[1], coeffs[2])
    ]

    mod = lcm(modules + coeffs)

    while True:
        choices = [choice(sets[b]) for b in bases]

        rem = crt(
            choices + coeffs_inv,
            bases + coeffs
        )

        if verbose:
            print(f'[*] Searching pseudoprime...', file=stderr)

        for i in range(iters):
            if verbose and i % 10000 == 0:
                print(f'{i}...')

            p1 = getrandbits(rlen) * mod + rem
            p2 = (p1 - 1) * coeffs[1] + 1
            p3 = (p1 - 1) * coeffs[2] + 1

            pprime = p1 * p2 * p3
            
            if miller_rabin(bases, pprime):
                break
        else:
            if verbose:
                print(f'[-] Failed to find pseudoprime, trying with another choices...', file=stderr)
            
            continue

        if verbose:
            print(f'[+] Found pseudoprime!', file=stderr)
            print(f'[+] P = {pprime}', file=stderr)

        return pprime, [p1, p2, p3]


if __name__ == '__main__':
    rlen = 256
    iters = 30000
    verbose = True

    bases = list(primes(300))
    coeffs = [1, 313, 353]

    pprime, divisors = search_pseudoprime(bases, coeffs, rlen, iters, verbose)
    
    assert not is_prime(pprime) and \
           miller_rabin(bases, pprime)

'''
$ sage arnault.sage
[*] Searching pseudoprime...
0...
10000...
20000...
[-] Failed to find pseudoprime, trying with another choices...
[*] Searching pseudoprime...
0...
10000...
[+] Found pseudoprime!
[+] P = 99597527340020670697596886062721977401836948352586238797499761849061796816245727295797460642211895009946326533856101876592304488359235447755504083536903673408562244316363452203072868521183142694959128745107323188995740668134018742165409361423628304730379121574707411453909999845745038957688998441109092021094925758212635651445626620045726265831347783805945477368631216031783484978212374792517000073275125176790602508815912876763504846656547041590967709195413101791490627310943998497788944526663960420235802025853374061708569334400472016398343229556656720912631463470998180176325607452843441554359644313713952036867
'''