import random
from math import floor, log

pto100 = [
    2,3,5,7,11,13,17,19,23,29,31,37,41,
     43,47,53,59,61,67,71,73,79,83,89,97
]

def in_100(x):
    """Check if x is in the primes less than 100 i.e., in the first 25 primes.
    Useful for e.g. the first step in Baillie–PSW.
    """
    return x in pto100

def hund_div(x):
    """Check if x is divisible by any of the 25 primes less than 100. Returns 0
    if there is no such divisor, otherwise returns the first divisor.
    """
    for v in pto100:
        if v >= x:
            return 0
        if x % v == 0:
            return v
    return 0

def trd(x):
    """Express x as 2^r * d. Returns two values, r and d.
    """
    count = 0
    xc = x
    while xc % 2 == 0:
        count = count + 1
        xc = xc // 2
    return count, xc

def MillerRabinRandom(n, k):
    """Carry out a Miller-Rabin test on n with k iterations, choosing random
    values for a.
    """
    r, d = trd(n-1)
    for iter in range(k):
        a = random.randint(2, n - 2)
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        compo = True
        for iterb in range(r - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                compo = False
                break
        if compo:
            return True
    return False

def MillerRabin(n, base):
    """Carry out a Miller-Rabin test on n with a specified base.
    """
    #print(n, base)
    r, d = trd(n-1)
    #print(r, d)
    a = base
    x = pow(a, d, n)
    #print(x)
    if x == 1 or x == n - 1:
        #print("First")
        return False
    for iterb in range(r - 1):
        x = pow(x, 2, n)
        if x == n - 1:
            #print("Second,", iterb, x)
            return False
    #print("True")
    return True

#for i in range(2, 100):
#    print(i, hund_div(i))
#
#print()

#Print pseudoprimes greater than 100 for each base, from base 1 to 15
for base in range(2, 16):
    print("{}:".format(base))
    for i in range(2, 10000):
        if base > i - 2 or i % 2 == 0:
            continue
        if MillerRabin(i, base) == (hund_div(i) == 0):
            print(i, hund_div(i), MillerRabin(i, base))
    print()

MillerRabin(28, 3)
print(pow(3, 27, 28))
print(pow(3, 27))
print(pow(3, 27) % 28)