Project P083
Clock arithmetic
Working modulo m: powers by repeated squaring step by step, inverses by the extended Euclidean algorithm, the Chinese remainder theorem for moduli that need not be coprime, and the Fermat and Miller-Rabin tests - where Carmichael numbers fool the first.
2 modules · 2 recorded sessionstext-menu UI in the terminalupdated 2026-10-15
Every screen below was recorded under CPython. When this page was built, the EML interpreter replayed each session from the same input and printed the same bytes.
About
Working modulo m. Powers by repeated squaring, shown step by step; inverses by the extended Euclidean algorithm; the Chinese remainder theorem; and the Fermat and Miller-Rabin primality tests, where the Carmichael numbers fool the first and not the second.
main.eml- the menu, reading numbers, and the reportsmodular.eml- exact halving and division, modular powers, the extended Euclidean algorithm, inverses, the Chinese remainder theorem, and the two tests
How each part works:
- Numbers go up to 10^18, past 2^53, so nothing goes through a float: halving and division are long division, digit by digit, and a % m is exact at any size.
- Powers, as in the corpus case
modular-exponentiation: go through the exponent's bits from the lowest, square the base each time and multiply it in for each 1 bit, reducing mod m every step. That is one squaring per bit, with every value below m. Done directly, 7^1000 would have 846 digits. - Inverse: the extended Euclidean algorithm finds x and y with a x + m y = gcd(a, m). When the gcd is 1, x is the inverse of a mod m; otherwise there is none.
- Chinese remainders: conditions are combined two at a time. Moduli that share a factor are allowed, and when two conditions disagree on that shared part, there is no answer at all.
- Fermat's test: a prime n has a^(n-1) = 1 (mod n) for every a it does not divide. Carmichael numbers such as 561 pass it for every such a.
- Miller-Rabin: write n - 1 = 2^s x d with d odd. For a prime, a^d is 1, or one of its squarings is n - 1. 561 fails that for a = 2. A few composites pass it for some bases (3215031751 for 2, 3, 5 and 7), so the verdict uses the first 13 primes as bases, which is known to be exact below 3.3 x 10^24 - far past the 10^18 the program accepts.
What is checked: menu choices 1 to 5; numbers from 0 (or 1, or 2) to 10^18; for the Chinese remainder theorem, 2 to 4 remainders and as many moduli from 2 to 1,000,000. An empty answer cancels.
Sessions: sessions/basic.in gives: - 7^1000 mod 13 step by step, and 2^(10^12) mod 1000000007; - the inverses of 17 mod 3120 (2753) and of 3 mod 1000000007; - x = 23 (mod 105) from three conditions, and 53 (mod 210) from four; - 561 (a Carmichael number), the primes 1000000007 and 999999999999999989, and 3215031751, which fools both tests for base 2 and is caught by base 11.
sessions/bad-input.in gives: - menu choices 0 and x; - a base of abc and of 10^18 + 1, an exponent of -1, a modulus of 1; - 6 mod 9, which has no inverse; - one and five remainders; one modulus for two remainders, a modulus over a million; and 1 mod 4 with 2 mod 6, which contradict each other; - 1, and the even number 4.
Built on the verified corpus case modular-exponentiation (repeated squaring, never holding a value above the modulus squared).
Recorded sessions
What the screen shows while someone uses the program. Each typed line appears after its prompt, the way a terminal shows it.
bad-input
interpreter: byte-equal== Clock arithmetic ==
Working modulo m: powers, inverses, the Chinese remainder theorem, and testing primes.
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> 0
Pick a number from 1 to 5.
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> x
Pick a number from 1 to 5.
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> 1
base (0 to 1000000000000000000)> abc
Type a whole number from 0 to 1000000000000000000.
base (0 to 1000000000000000000)> 1000000000000000001
Type a whole number from 0 to 1000000000000000000.
base (0 to 1000000000000000000)> 5
exponent (0 to 1000000000000000000)> -1
Type a whole number from 0 to 1000000000000000000.
exponent (0 to 1000000000000000000)> 3
modulus (2 to 1000000000000000000)> 1
Type a whole number from 2 to 1000000000000000000.
modulus (2 to 1000000000000000000)>
Cancelled.
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> 2
number (1 to 1000000000000000000)> 6
modulus (2 to 1000000000000000000)> 9
6 = 0 x 9 + 6
9 = 1 x 6 + 3
6 = 2 x 3 + 0
6 has no inverse modulo 9: they share the factor 3.
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> 3
remainders (2 to 4, like 2 3 2)> 1
Type 2 to 4 whole numbers.
remainders (2 to 4, like 2 3 2)> 1 2 3 4 5
Type 2 to 4 whole numbers.
remainders (2 to 4, like 2 3 2)> 1 2
moduli, as many, each 2 to 1000000 (like 3 5 7)> 3
Type 2 moduli from 2 to 1000000.
moduli, as many, each 2 to 1000000 (like 3 5 7)> 1 1000001
Type 2 moduli from 2 to 1000000.
moduli, as many, each 2 to 1000000 (like 3 5 7)> 4 6
x = 1 (mod 4), x = 2 (mod 6)
No number does all of that: two of the conditions contradict each other where their moduli overlap.
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> 4
number (2 to 1000000000000000000)> 1
Type a whole number from 2 to 1000000000000000000.
number (2 to 1000000000000000000)> 4
4 is even, so it is not prime.
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> 5
Bye.
What was typed (24 lines)
0
x
1
abc
1000000000000000001
5
-1
3
1
2
6
9
3
1
1 2 3 4 5
1 2
3
1 1000001
4 6
4
1
4
5
basic
interpreter: byte-equal== Clock arithmetic ==
Working modulo m: powers, inverses, the Chinese remainder theorem, and testing primes.
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> 1
base (0 to 1000000000000000000)> 7
exponent (0 to 1000000000000000000)> 1000
modulus (2 to 1000000000000000000)> 13
bit base (squared each step) result so far
0 7 1
0 10 1
0 9 1
1 3 3
0 9 3
1 3 9
1 9 3
1 3 9
1 9 3
1 3 9
7^1000 mod 13 = 9
One squaring for each of the exponent's 10 bits, every value kept below the modulus, instead of multiplying 1000 times.
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> 1
base (0 to 1000000000000000000)> 2
exponent (0 to 1000000000000000000)> 1000000000000
modulus (2 to 1000000000000000000)> 1000000007
2^1000000000000 mod 1000000007 = 959366170
One squaring for each of the exponent's 40 bits, every value kept below the modulus, instead of multiplying 1000000000000 times.
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> 2
number (1 to 1000000000000000000)> 17
modulus (2 to 1000000000000000000)> 3120
17 = 0 x 3120 + 17
3120 = 183 x 17 + 9
17 = 1 x 9 + 8
9 = 1 x 8 + 1
8 = 8 x 1 + 0
Working back up the divisions: 17 x 2753 = 1 (mod 3120), so the inverse is 2753.
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> 2
number (1 to 1000000000000000000)> 3
modulus (2 to 1000000000000000000)> 1000000007
3 = 0 x 1000000007 + 3
1000000007 = 333333335 x 3 + 2
3 = 1 x 2 + 1
2 = 2 x 1 + 0
Working back up the divisions: 3 x 333333336 = 1 (mod 1000000007), so the inverse is 333333336.
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> 3
remainders (2 to 4, like 2 3 2)> 2 3 2
moduli, as many, each 2 to 1000000 (like 3 5 7)> 3 5 7
x = 2 (mod 3), x = 3 (mod 5), x = 2 (mod 7)
x = 23 (mod 105): 23, 128, 233, ...
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> 3
remainders (2 to 4, like 2 3 2)> 1 2 3 4
moduli, as many, each 2 to 1000000 (like 3 5 7)> 2 3 5 7
x = 1 (mod 2), x = 2 (mod 3), x = 3 (mod 5), x = 4 (mod 7)
x = 53 (mod 210): 53, 263, 473, ...
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> 4
number (2 to 1000000000000000000)> 561
Fermat, base 2: 2^(n-1) mod n = 1 - it passes.
Miller-Rabin, base 2: n - 1 = 2^4 x 35; 2^35 and its squarings mod n: 263, 166, 67, 1
A prime would give 1 at the start or n - 1 on the way: so it is composite, and base 2 is the witness.
561 is composite.
Base 2 is a Fermat liar for it.
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> 4
number (2 to 1000000000000000000)> 1000000007
Fermat, base 2: 2^(n-1) mod n = 1 - it passes.
Miller-Rabin, base 2: n - 1 = 2^1 x 500000003; 2^500000003 and its squarings mod n: 1
That is what a prime gives (1 at the start, or n - 1 on the way): it passes.
1000000007 is prime: it passes Miller-Rabin for the 13 bases 2 to 41, which is exact below 3.3 x 10^24.
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> 4
number (2 to 1000000000000000000)> 3215031751
Fermat, base 2: 2^(n-1) mod n = 1 - it passes.
Miller-Rabin, base 2: n - 1 = 2^1 x 1607515875; 2^1607515875 and its squarings mod n: 1
That is what a prime gives (1 at the start, or n - 1 on the way): it passes.
3215031751 is composite.
Base 2 is a Fermat liar for it.
Base 2 is a strong liar too; base 11 is the first of the 13 to catch it.
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> 4
number (2 to 1000000000000000000)> 999999999999999989
Fermat, base 2: 2^(n-1) mod n = 1 - it passes.
Miller-Rabin, base 2: n - 1 = 2^2 x 249999999999999997; 2^249999999999999997 and its squarings mod n: 639232157799794134, 999999999999999988
That is what a prime gives (1 at the start, or n - 1 on the way): it passes.
999999999999999989 is prime: it passes Miller-Rabin for the 13 bases 2 to 41, which is exact below 3.3 x 10^24.
1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit
choice> 5
Bye.
What was typed (29 lines)
1
7
1000
13
1
2
1000000000000
1000000007
2
17
3120
2
3
1000000007
3
2 3 2
3 5 7
3
1 2 3 4
2 3 5 7
4
561
4
1000000007
4
3215031751
4
999999999999999989
5
Modules
The program as written, entry module first. Each module transpiles to its own Python file, which is what eml project run executes.
main.eml(entry)
eml# P083 clock arithmetic: powers modulo m step by step, inverses by the
# extended Euclidean algorithm, the Chinese remainder theorem, and the
# Fermat and Miller-Rabin primality tests - where the Carmichael numbers
# fool the first and not the second.
import modular
def trim(s):
0 => i
len(s) => j
while i < j and s[i] == " ":
i + 1 => i
while j > i and s[j - 1] == " ":
j - 1 => j
return s[i:j]
def whole(s):
# A whole number of up to 18 digits, or -1.
if len(s) == 0 or len(s) > 18:
return -1
0 => n
for c in s:
if not (c >= "0" and c <= "9"):
return -1
n * 10 + int(c) => n
return n
def ask(prompt, low, high):
while True:
trim(input(prompt + " (" + str(low) + " to " + str(high) + ")> ")) => a
if a == "":
return -1
whole(a) => n
if n >= low and n <= high:
return n
("Type a whole number from " + str(low) + " to " + str(high) + ".") ^0
def numbers(s):
[] => out
"" => w
for c in s + " ":
if c == " " or c == ",":
if w != "":
whole(w) => v
if v < 0:
return []
out + [v] => out
"" => w
else:
w + c => w
return out
def lpad(s, n):
while len(s) < n:
" " + s => s
return s
def pow_mod():
ask("base", 0, 1000000000000000000) => a
if a == -1:
"Cancelled." ^0
return
ask("exponent", 0, 1000000000000000000) => e
if e == -1:
"Cancelled." ^0
return
ask("modulus", 2, 1000000000000000000) => m
if m == -1:
"Cancelled." ^0
return
modular.power(a, e, m) => r
if len(r[1]) <= 20:
" bit base (squared each step) result so far" ^0
for st in r[1]:
(" " + lpad(str(st[0]), 3) + " " + lpad(str(st[1]), 24) + " " + lpad(str(st[2]), 19)) ^0
(str(a) + "^" + str(e) + " mod " + str(m) + " = " + str(r[0])) ^0
("One squaring for each of the exponent's " + str(len(r[1])) + " bits, every value kept below the modulus, instead of multiplying " + str(e) + " times.") ^0
def inverse():
ask("number", 1, 1000000000000000000) => a
if a == -1:
"Cancelled." ^0
return
ask("modulus", 2, 1000000000000000000) => m
if m == -1:
"Cancelled." ^0
return
modular.extended(a % m, m) => e
if len(e[3]) <= 20:
for row in e[3]:
(" " + str(row[0]) + " = " + str(row[2]) + " x " + str(row[1]) + " + " + str(row[0] - row[2] * row[1])) ^0
if e[0] != 1:
(str(a) + " has no inverse modulo " + str(m) + ": they share the factor " + str(e[0]) + ".") ^0
return
modular.inverse(a, m) => x
("Working back up the divisions: " + str(a) + " x " + str(x) + " = 1 (mod " + str(m) + "), so the inverse is " + str(x) + ".") ^0
def remainders():
while True:
trim(input("remainders (2 to 4, like 2 3 2)> ")) => s
if s == "":
"Cancelled." ^0
return
numbers(s) => rs
if len(rs) >= 2 and len(rs) <= 4:
while True:
trim(input("moduli, as many, each 2 to 1000000 (like 3 5 7)> ")) => t
if t == "":
"Cancelled." ^0
return
numbers(t) => ms
True => ok
for m in ms:
if m < 2 or m > 1000000:
False => ok
if ok and len(ms) == len(rs):
"" => line
for i in [0:len(ms) - 1]:
if line != "":
line + ", " => line
line + "x = " + str(rs[i]) + " (mod " + str(ms[i]) + ")" => line
line ^0
modular.crt(rs, ms) => r
if len(r) == 0:
"No number does all of that: two of the conditions contradict each other where their moduli overlap." ^0
else:
("x = " + str(r[0]) + " (mod " + str(r[1]) + "): " + str(r[0]) + ", " + str(r[0] + r[1]) + ", " + str(r[0] + 2 * r[1]) + ", ...") ^0
return
("Type " + str(len(rs)) + " moduli from 2 to 1000000.") ^0
else:
"Type 2 to 4 whole numbers." ^0
def prime():
ask("number", 2, 1000000000000000000) => n
if n == -1:
"Cancelled." ^0
return
if n == 2:
"2 is prime." ^0
return
if n % 2 == 0:
(str(n) + " is even, so it is not prime.") ^0
return
modular.power(2, n - 1, n)[0] => f
if f == 1:
("Fermat, base 2: 2^(n-1) mod n = 1 - it passes.") ^0
else:
("Fermat, base 2: 2^(n-1) mod n = " + str(f) + ", not 1 - so it is composite.") ^0
modular.strong(n, 2) => mr
"" => line
for x in mr[3]:
if line != "":
line + ", " => line
line + str(x) => line
("Miller-Rabin, base 2: n - 1 = 2^" + str(mr[1]) + " x " + str(mr[2]) + "; 2^" + str(mr[2]) + " and its squarings mod n: " + line) ^0
if mr[0]:
" That is what a prime gives (1 at the start, or n - 1 on the way): it passes." ^0
else:
" A prime would give 1 at the start or n - 1 on the way: so it is composite, and base 2 is the witness." ^0
modular.is_prime(n) => p
if p:
(str(n) + " is prime: it passes Miller-Rabin for the 13 bases 2 to 41, which is exact below 3.3 x 10^24.") ^0
else:
(str(n) + " is composite.") ^0
if f == 1:
"Base 2 is a Fermat liar for it." ^0
if mr[0]:
0 => w
for b in [3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]:
if w == 0 and n % b != 0 and not modular.strong(n, b)[0]:
b => w
if w > 0:
("Base 2 is a strong liar too; base " + str(w) + " is the first of the 13 to catch it.") ^0
else:
"Base 2 is a strong liar too; a small factor gives it away." ^0
"== Clock arithmetic ==" ^0
"Working modulo m: powers, inverses, the Chinese remainder theorem, and testing primes." ^0
True => running
while running:
"" ^0
"1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit" ^0
trim(input("choice> ")) => choice
if choice == "1":
pow_mod()
elif choice == "2":
inverse()
elif choice == "3":
remainders()
elif choice == "4":
prime()
elif choice == "5":
False => running
else:
"Pick a number from 1 to 5." ^0
"Bye." ^0
Python projection (main.py)
import modular
def trim(s):
i = 0
j = len(s)
while i < j and s[i] == " ":
i = i + 1
while j > i and s[j - 1] == " ":
j = j - 1
return s[i:j]
def whole(s):
if len(s) == 0 or len(s) > 18:
return -1
n = 0
for c in s:
if not (c >= "0" and c <= "9"):
return -1
n = n * 10 + int(c)
return n
def ask(prompt, low, high):
while True:
a = trim(input(prompt + " (" + str(low) + " to " + str(high) + ")> "))
if a == "":
return -1
n = whole(a)
if n >= low and n <= high:
return n
print("Type a whole number from " + str(low) + " to " + str(high) + ".")
def numbers(s):
out = []
w = ""
for c in s + " ":
if c == " " or c == ",":
if w != "":
v = whole(w)
if v < 0:
return []
out = out + [v]
w = ""
else:
w = w + c
return out
def lpad(s, n):
while len(s) < n:
s = " " + s
return s
def pow_mod():
a = ask("base", 0, 1000000000000000000)
if a == -1:
print("Cancelled.")
return
e = ask("exponent", 0, 1000000000000000000)
if e == -1:
print("Cancelled.")
return
m = ask("modulus", 2, 1000000000000000000)
if m == -1:
print("Cancelled.")
return
r = modular.power(a, e, m)
if len(r[1]) <= 20:
print(" bit base (squared each step) result so far")
for st in r[1]:
print(" " + lpad(str(st[0]), 3) + " " + lpad(str(st[1]), 24) + " " + lpad(str(st[2]), 19))
print(str(a) + "^" + str(e) + " mod " + str(m) + " = " + str(r[0]))
print("One squaring for each of the exponent's " + str(len(r[1])) + " bits, every value kept below the modulus, instead of multiplying " + str(e) + " times.")
def inverse():
a = ask("number", 1, 1000000000000000000)
if a == -1:
print("Cancelled.")
return
m = ask("modulus", 2, 1000000000000000000)
if m == -1:
print("Cancelled.")
return
e = modular.extended(a % m, m)
if len(e[3]) <= 20:
for row in e[3]:
print(" " + str(row[0]) + " = " + str(row[2]) + " x " + str(row[1]) + " + " + str(row[0] - row[2] * row[1]))
if e[0] != 1:
print(str(a) + " has no inverse modulo " + str(m) + ": they share the factor " + str(e[0]) + ".")
return
x = modular.inverse(a, m)
print("Working back up the divisions: " + str(a) + " x " + str(x) + " = 1 (mod " + str(m) + "), so the inverse is " + str(x) + ".")
def remainders():
while True:
s = trim(input("remainders (2 to 4, like 2 3 2)> "))
if s == "":
print("Cancelled.")
return
rs = numbers(s)
if len(rs) >= 2 and len(rs) <= 4:
while True:
t = trim(input("moduli, as many, each 2 to 1000000 (like 3 5 7)> "))
if t == "":
print("Cancelled.")
return
ms = numbers(t)
ok = True
for m in ms:
if m < 2 or m > 1000000:
ok = False
if ok and len(ms) == len(rs):
line = ""
for i in range(0, len(ms)):
if line != "":
line = line + ", "
line = line + "x = " + str(rs[i]) + " (mod " + str(ms[i]) + ")"
print(line)
r = modular.crt(rs, ms)
if len(r) == 0:
print("No number does all of that: two of the conditions contradict each other where their moduli overlap.")
else:
print("x = " + str(r[0]) + " (mod " + str(r[1]) + "): " + str(r[0]) + ", " + str(r[0] + r[1]) + ", " + str(r[0] + 2 * r[1]) + ", ...")
return
print("Type " + str(len(rs)) + " moduli from 2 to 1000000.")
else:
print("Type 2 to 4 whole numbers.")
def prime():
n = ask("number", 2, 1000000000000000000)
if n == -1:
print("Cancelled.")
return
if n == 2:
print("2 is prime.")
return
if n % 2 == 0:
print(str(n) + " is even, so it is not prime.")
return
f = modular.power(2, n - 1, n)[0]
if f == 1:
print("Fermat, base 2: 2^(n-1) mod n = 1 - it passes.")
else:
print("Fermat, base 2: 2^(n-1) mod n = " + str(f) + ", not 1 - so it is composite.")
mr = modular.strong(n, 2)
line = ""
for x in mr[3]:
if line != "":
line = line + ", "
line = line + str(x)
print("Miller-Rabin, base 2: n - 1 = 2^" + str(mr[1]) + " x " + str(mr[2]) + "; 2^" + str(mr[2]) + " and its squarings mod n: " + line)
if mr[0]:
print(" That is what a prime gives (1 at the start, or n - 1 on the way): it passes.")
else:
print(" A prime would give 1 at the start or n - 1 on the way: so it is composite, and base 2 is the witness.")
p = modular.is_prime(n)
if p:
print(str(n) + " is prime: it passes Miller-Rabin for the 13 bases 2 to 41, which is exact below 3.3 x 10^24.")
else:
print(str(n) + " is composite.")
if f == 1:
print("Base 2 is a Fermat liar for it.")
if mr[0]:
w = 0
for b in [3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]:
if w == 0 and n % b != 0 and not modular.strong(n, b)[0]:
w = b
if w > 0:
print("Base 2 is a strong liar too; base " + str(w) + " is the first of the 13 to catch it.")
else:
print("Base 2 is a strong liar too; a small factor gives it away.")
print("== Clock arithmetic ==")
print("Working modulo m: powers, inverses, the Chinese remainder theorem, and testing primes.")
running = True
while running:
print("")
print("1) power mod m 2) inverse mod m 3) Chinese remainders 4) is it prime? 5) quit")
choice = trim(input("choice> "))
if choice == "1":
pow_mod()
elif choice == "2":
inverse()
elif choice == "3":
remainders()
elif choice == "4":
prime()
elif choice == "5":
running = False
else:
print("Pick a number from 1 to 5.")
print("Bye.")
modular.eml
eml# P083 clock arithmetic - working modulo m: powers by repeated squaring, the
# inverse by the extended Euclidean algorithm, the Chinese remainder
# theorem, and two primality tests. Numbers go far past 2^53, so halving is
# exact long division and nothing goes through a float; a % b is exact at
# any size.
def half(n):
# n // 2 for n >= 0, digit by digit.
0 => q
0 => r
for ch in str(n):
r * 10 + int(ch) => r
0 => k
while r >= 2:
r - 2 => r
k + 1 => k
q * 10 + k => q
return q
def quotient(a, b):
# a // b for a >= 0, b > 0, digit by digit (at most 9 subtractions each).
0 => q
0 => r
for ch in str(a):
r * 10 + int(ch) => r
0 => k
while r >= b:
r - b => r
k + 1 => k
q * 10 + k => q
return q
def power(a, e, m):
# The corpus case modular-exponentiation: a^e mod m, squaring the base
# and multiplying it in for each 1 bit of e, from the lowest. Returns
# [result, steps], each step [bit, base before squaring, result so far].
1 % m => out
a % m => b
e => left
[] => steps
while left > 0:
left % 2 => bit
if bit == 1:
out * b % m => out
steps + [[bit, b, out]] => steps
half(left) => left
b * b % m => b
return [out, steps]
def extended(a, b):
# [g, x, y] with a x + b y = g = gcd(a, b), and the rows of the
# algorithm: [a, b, q] for each division.
[] => rows
1 => x0
0 => y0
0 => x1
1 => y1
while b != 0:
quotient(a, b) => q
rows + [[a, b, q]] => rows
a - q * b => r
b => a
r => b
x0 - q * x1 => t
x1 => x0
t => x1
y0 - q * y1 => t
y1 => y0
t => y1
return [a, x0, y0, rows]
def inverse(a, m):
# The x in 0..m-1 with a x = 1 (mod m), or -1 when gcd(a, m) is not 1.
extended(a % m, m) => e
if e[0] != 1:
return -1
return e[1] % m
def crt(rs, ms):
# x = r_i (mod m_i) for all i, combined two at a time. The moduli need
# not be coprime; when two conditions contradict there is no answer.
# Returns [x, modulus] or [] if none.
rs[0] % ms[0] => x
ms[0] => m
for i in [1:len(ms) - 1]:
ms[i] => n
rs[i] % n => r
extended(m, n) => e
e[0] => g
if (r - x) % g != 0:
return []
quotient(n, g) => step
# x + m t = r (mod n): t = ((r - x)/g) * (m/g)^-1 (mod n/g)
(r - x) % n => diff
quotient(diff, g) => k
inverse(quotient(m, g) % step, step) => inv
if step == 1:
0 => t
else:
k * inv % step => t
x + m * t => x
m * step => m
x % m => x
return [x, m]
def fermat(n, a):
# True when a^(n-1) = 1 (mod n): every prime passes for every a it does
# not divide; a composite that passes is a Fermat liar for base a.
return power(a, n - 1, n)[0] == 1
def strong(n, a):
# The Miller-Rabin test for base a, n odd > 2: n - 1 = 2^s d with d odd;
# a prime gives a^d = 1, or -1 at some squaring. Returns [passes, s, d,
# the sequence a^d, a^(2d), ... mod n].
n - 1 => d
0 => s
while d % 2 == 0:
half(d) => d
s + 1 => s
power(a, d, n)[0] => x
[x] => seq
(x == 1 or x == n - 1) => ok
for k in [1:s - 1]:
x * x % n => x
seq + [x] => seq
if x == n - 1:
True => ok
return [ok, s, d, seq]
def is_prime(n):
# Miller-Rabin with the first 13 primes as bases, which is known to be
# exact for every n below 3,317,044,064,679,887,385,961,981 (Sorenson
# and Webster, 2015) - far past the 10^18 the program accepts.
if n < 2:
return False
for p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]:
if n == p:
return True
if n % p == 0:
return False
for p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]:
if not strong(n, p)[0]:
return False
return True
Python projection (modular.py)
def half(n):
q = 0
r = 0
for ch in str(n):
r = r * 10 + int(ch)
k = 0
while r >= 2:
r = r - 2
k = k + 1
q = q * 10 + k
return q
def quotient(a, b):
q = 0
r = 0
for ch in str(a):
r = r * 10 + int(ch)
k = 0
while r >= b:
r = r - b
k = k + 1
q = q * 10 + k
return q
def power(a, e, m):
out = 1 % m
b = a % m
left = e
steps = []
while left > 0:
bit = left % 2
if bit == 1:
out = out * b % m
steps = steps + [[bit, b, out]]
left = half(left)
b = b * b % m
return [out, steps]
def extended(a, b):
rows = []
x0 = 1
y0 = 0
x1 = 0
y1 = 1
while b != 0:
q = quotient(a, b)
rows = rows + [[a, b, q]]
r = a - q * b
a = b
b = r
t = x0 - q * x1
x0 = x1
x1 = t
t = y0 - q * y1
y0 = y1
y1 = t
return [a, x0, y0, rows]
def inverse(a, m):
e = extended(a % m, m)
if e[0] != 1:
return -1
return e[1] % m
def crt(rs, ms):
x = rs[0] % ms[0]
m = ms[0]
for i in range(1, len(ms)):
n = ms[i]
r = rs[i] % n
e = extended(m, n)
g = e[0]
if (r - x) % g != 0:
return []
step = quotient(n, g)
diff = (r - x) % n
k = quotient(diff, g)
inv = inverse(quotient(m, g) % step, step)
if step == 1:
t = 0
else:
t = k * inv % step
x = x + m * t
m = m * step
x = x % m
return [x, m]
def fermat(n, a):
return power(a, n - 1, n)[0] == 1
def strong(n, a):
d = n - 1
s = 0
while d % 2 == 0:
d = half(d)
s = s + 1
x = power(a, d, n)[0]
seq = [x]
ok = x == 1 or x == n - 1
for k in range(1, s):
x = x * x % n
seq = seq + [x]
if x == n - 1:
ok = True
return [ok, s, d, seq]
def is_prime(n):
if n < 2:
return False
for p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]:
if n == p:
return True
if n % p == 0:
return False
for p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]:
if not strong(n, p)[0]:
return False
return True