Miller Rabin (may still be buggy)

bot.py

@@ -1,3 +1,94 @@
 import random
+from math import floor, log
 
-print(random.randint(0,12))
+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)