1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128

#!/usr/bin/python3
 
# SI 335: Computer Algorithms
# Unit 3
 
import random
 
def leastPrimeFactor(n):
    '''Given a positive integer n, returns the smallest prime
       number p which divides n'''
    i = 2
    while i * i <= n:
        if n % i == 0:
            return i
        i = i + 1
    # If we get here, n must be a prime number itself.
    return n
 
def gcd(a, b):
    if b == 0:
        return a
    else:
        return gcd(b, a % b)
 
def gcditer(a, b):
    while b != 0:
        (a, b) = (b, a % b)
    return a
 
def xgcd(a, b):
    if b == 0:
        return (a, 1, 0)
    else:
        q, r = divmod(a, b)
        (g, s0, t0) = xgcd(b, r)
        return (g, t0, s0 - t0*q)
 
def probably_prime(n):
    # This function returns a random integer from 2 up to n-2.
    a = random.randrange(2, n-1)
    d = n-1
    k = 0
    while d % 2 == 0:
        d = d // 2
        k = k + 1
    # IMPORTANT: This next line should be done more efficiently!!
    x = a**d % n # ** is the exponentiation operator
    if x**2 % n == 1:
        return True
    for r in range(1, k):
        x = x**2 % n
        if x == 1:
            return False
        if x == n-1:
            return True
    return False
 
 
# The rest is just to do some "sanity checks" for myself.
 
primes = {863, 1783, 2801, 5431, 14321, 45053, 7, 19387, 2351}
compos = {196, 1306, 3423, 6160, 10983, 5860, 61273, 56058, 87201}
 
def lpftest(alg):
    for p in primes:
        if alg(p) != p:
            print("Prime FAIL for", p)
            return False
    for n in compos:
        k = alg(n)
        if k <= 1 or k >= n or n%k != 0:
            print("Composite FAIL for", n)
            return False
    return True
 
def prtest(alg):
    for p in primes:
        if not alg(p):
            print("MR FAIL for prime", p)
            return False
    for n in compos:
        if alg(n):
            print("MR FAIL for composite", n)
            return False
    return True
 
def gcdtest(alg):
    for i in range(1000):
        a, b, c = (random.randrange(10**5) for i in range(3))
        x, y = a*c, b*c
        g = alg(x, y)
        if g < 1 or g%c != 0 or x%g != 0 or y%g != 0:
            print("gcd FAIL for {} on x={} and y={}".format(alg.__name__,x,y))
            return False
        if alg(x//g, y//g) != 1:
            print("gcd FAIL2 for {} on x={} and y={}".format(alg.__name__,x,y))
            return False
    return True
 
def xgcdtest(alg):
    for i in range(1000):
        a, b, c = (random.randrange(10**5) for i in range(3))
        x, y = a*c, b*c
        g, s, t = alg(x, y)
        if g != gcd(x, y):
            print("{} FAIL wrong gcd for x={}, y={}".format(alg.__name__,x,y))
            return False
        if s*x + t*y != g or abs(s) > y or abs(t) > x:
            print("{} FAIL wrong s,t for x={}, y={}".format(alg.__name__,x,y))
            return False
    return True
 
if __name__ == '__main__':
    tests = [
        (lpftest, leastPrimeFactor),
        (gcdtest, gcd),
        (gcdtest, gcditer),
        (xgcdtest, xgcd),
        (prtest, probably_prime),
    ]
    allGood = True
    for tester, alg in tests:
        allGood = allGood and tester(alg)
    if allGood:
        print("All {} checks passed!".format(len(tests)))
        exit(0)
    else:
        exit(1)