<!-- canonical: https://efficientnewlanguage.org/eml-p/projects/P083-clock-arithmetic/ | updated: 2026-10-15 -->

# 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.

EML-P project `projects/clock-arithmetic` in the EML language repo: 2 module(s), entry `main.eml`, terminal UI. There, `eml project run projects/clock-arithmetic` runs it and `eml project verify projects/clock-arithmetic` replays every session under CPython (two hash seeds) and in the interpreter; the site build replays every session in the interpreter again and publishes a session only if its screen matches.

Built on verified corpus cases: modular-exponentiation (https://efficientnewlanguage.org/cases/117-modular-exponentiation/).

## Sessions

### bad-input - interpreter: byte-equal to the golden

Input:

```text
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
```

Screen:

```text
== 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.
```

### basic - interpreter: byte-equal to the golden

Input:

```text
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
```

Screen:

```text
== 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.
```

## Modules

### main.eml

```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 of main.eml:

```python
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 of modular.eml:

```python
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
```

## README

# P083 - Clock arithmetic

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 reports
- `modular.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).