Miller-Rabin Primality Test-Witnesses and Liars - Implementing in Python

Question

I have been studying the Miller-Rabin Primality Test, and am interested in implementing a code in Python to count witnesses and liars. In the basic code to determine if a number is probably prime or composite, I would like to incorporate 1) and 2) below to better understand witnesses compared to liars for n values:

$1)$ $a$ to be tested for all values of $a$ from $1<a<n-1$, not random as it is currently coded.

$2)$ Then for every $a$ in $1)$ above, a count of how many of those $a's$ are witnesses and how many are non-witnesses(liars)

My ultimate goal is to use this, I'm sure with more modifications to the code, to compare to the Theorem: If n is odd, composite, and n>9, then at least 75% of the elements in $(\mathbb Z/n\mathbb Z)^x$ are Miller-Rabin witnesses.

The Python code I am using is as follows:

from random import randrange

def probably_prime(n, k):
"""Return True if n passes k rounds of the Miller-Rabin primality
test (and is probably prime). Return False if n is proved to be
composite.

"""
if n < 2: return False

r, m = 0, n - 1
while m % 2 == 0:
    r += 1
    m //= 2
for _ in range(k):
    a = randrange(2, n - 1)
    x = pow(a, m, n)
    if x == 1 or x == n - 1:
        continue
    for _ in range(r - 1):
        x = pow(x, 2, n)
        if x == n - 1:
            break
    else:
        return False
return True

Answer 1

For 1) You can change the code to:

def probably_prime(n, k, a):
"""Return True if n passes k rounds of the Miller-Rabin primality
test (and is probably prime). Return False if n is proved to be
composite.

"""
if n < 2: return False

r, m = 0, n - 1
while m % 2 == 0:
    r += 1
    m //= 2
for _ in range(k):
    #a = randrange(2, n - 1)
    x = pow(a, m, n)
    if x == 1 or x == n - 1:
        continue
    for _ in range(r - 1):
        x = pow(x, 2, n)
        if x == n - 1:
            break
    else:
        return False
return True

and then fix some number n, and some k, which is the number of repetitions. Then add at the end:

Result=[]
For a in range(1,n):
    Result.append(probably_prime(n,k,a))

Result is a list with all outcomes for $a$ in the required range. True means composite (ie witness), and false means liar. To find which $a$'s are liars and which are witnesses you can then add:

witnesses=sum(Result)
liars=n-witnesses

Answer 2

A simple code that returns a list of liars or 'None' if n is prime is:

# Miller-Rabin-Test of n, base a
def miller_rabin(n, a):
    d = n1 = n-1
    s = 0
    while d & 1 == 0:
        d >>= 1
        s += 1
    b = pow(a, d, n)
    if b == 1 or b == n1: return True
    while s > 1:
        s -= 1
        b = (b*b) % n
        if b == n1: return True 
    return False

# perform Miller-Rabin-Test of n with bases 1 < a < n-1 and return a list of strong liars (false witnesses) or 'None' if n is prime
def get_liars(n):
    result = []
    passed = 0
    for a in range(2, n-1):
        if miller_rabin(n, a):
            passed += 1
            if passed > n//4:
                return None    # n is prime
            result += [a]
    return result

# main program for testing
if __name__ == "__main__":
    while True:
        n = int(input("n = "))
        if n == 0: break
        if n < 3: continue
        liars = get_liars(n)
        if liars is None:
            print (n, "is prime")
        else:
            print (n, f"is composite, {len(liars)} liars:", liars)

The output for $n = 65$ would be

65 is composite, 4 liars: [8, 18, 47, 57]

The function get_liars returns None if the test has been passed with more than $n/4$ bases because $n$ is certainly a prime then.