Move some more generally useful functions and objects into a new file

millerrabin.py

@@ -1,46 +1,13 @@
 import random
+import psimp as p
 
-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)
+    r, d = p.trd(n-1)
     for iter in range(k):
         a = random.randint(2, n - 2)
         x = pow(a, d, n)
@@ -60,7 +27,7 @@ def MillerRabinRandom(n, k):
 def MillerRabin(n, base):
     """Carry out a Miller-Rabin test on n with a specified base.
     """
-    r, d = trd(n-1)
+    r, d = p.trd(n-1)
     a = base
     x = pow(a, d, n)
     if x == 1 or x == n - 1:
@@ -78,8 +45,8 @@ for base in range(2, 16):
     for i in range(2, 10000):
         if i % 2 == 0 or base % i == 0:
             continue
-        if MillerRabin(i, base) == (hund_div(i) == 0):
-            print(i, hund_div(i), MillerRabin(i, base))
+        if MillerRabin(i, base) == (p.hund_div(i) == 0):
+            print(i, p.hund_div(i), MillerRabin(i, base))
     print()
 
 MillerRabin(28, 3)

psimp.py

@@ -0,0 +1,34 @@
+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