# Modular Arithmetic and Check Digits {#sec-modular}
Here is a question that sounds trivial: if today is Wednesday, what day
will it be 100 days from now? You do not need to count 100 boxes on a
calendar. Compute $100 \bmod 7 = 2$, add 2 days to Wednesday, and the
answer is Friday. That single step --- dividing and keeping only the
remainder --- is **modular arithmetic**.
Its reach goes far beyond day-of-the-week puzzles. The check digit at
the end of every book ISBN, every grocery barcode, and every credit card
number is a modular sum designed to catch typing errors. The
cryptographic algorithm protecting your internet passwords rests on a
theorem about prime remainders that Fermat wrote in a letter in 1640 ---
without a computer. In this chapter we build the machinery from scratch
and watch it work.
## Clocks and Remainders {#sec-modular-clocks}
A clock is the original modular arithmetic machine. After 12 o'clock
comes 1, not 13 --- the count wraps. The calendar works the same way:
days of the week cycle with period 7, so "100 days from Wednesday" is
just a wrapping addition.
### Research Example: How Does a Clock Do Arithmetic? {.unnumbered .unlisted}
Can a picture of a clock show why $100 \bmod 7 = 2$ lands us on Friday? We plot the seven days as dots on a circle and draw the arrow that jumps exactly two steps — turning an abstract remainder into a concrete visual leap.
```{python}
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
days = ['Sun', 'Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat']
n = len(days)
# angles: place Sun at top, go clockwise
angles = [2 * np.pi * i / n - np.pi / 2 for i in range(n)]
r_dot, r_lbl = 1.0, 1.35
fig, ax = plt.subplots(figsize=(5, 5))
ax.set_aspect('equal')
ax.axis('off')
ax.add_patch(plt.Circle((0, 0), r_dot,
fill=False, lw=1.5, color='#AAAAAA'))
for i, (day, ang) in enumerate(zip(days, angles)):
xd = r_dot * np.cos(ang)
yd = r_dot * np.sin(ang)
xl = r_lbl * np.cos(ang)
yl = r_lbl * np.sin(ang)
hi = day in ('Wed', 'Fri')
color = '#2C6FAC' if hi else '#888888'
ax.plot(xd, yd, 'o', color=color, ms=11, zorder=3)
ax.text(xl, yl, day, ha='center', va='center',
fontsize=10, color=color,
fontweight='bold' if hi else 'normal')
# arc arrow from Wed (index 3) to Fri (index 5)
wed, fri = angles[3], angles[5]
ax.annotate(
'', xy=(r_dot * np.cos(fri), r_dot * np.sin(fri)),
xytext=(r_dot * np.cos(wed), r_dot * np.sin(wed)),
arrowprops=dict(arrowstyle='->', color='#2C6FAC',
lw=2, connectionstyle='arc3,rad=-0.35')
)
ax.text(0.28, 0.05, '+2', fontsize=12, color='#2C6FAC',
ha='center', va='center', fontweight='bold')
ax.set_xlim(-1.8, 1.8)
ax.set_ylim(-1.8, 1.8)
ax.set_title('100 mod 7 = 2: Wed + 2 = Fri', fontsize=12)
fig.tight_layout()
plt.show()
```
The arrow from Wed to Fri says everything: instead of counting 100 boxes on a calendar, you compute one remainder and jump that many steps. That is all modular arithmetic ever is — and 25 lines of Python just made it visible.
For any integer $a$ and any positive integer $m$, there exist unique
integers $q$ (the quotient) and $r$ (the remainder) with
$$a = q \cdot m + r, \qquad 0 \leq r < m.$$
We say $a$ and $b$ are **congruent modulo $m$**, written
$a \equiv b \pmod{m}$, when they leave the same remainder after dividing
by $m$ --- equivalently, when $m$ divides $a - b$. The number $m$ is
called the **modulus**. Every integer is congruent to exactly one value
in $\{0, 1, 2, \ldots, m-1\}$, so these $m$ classes cover all the
integers without overlap.
The key algebraic fact is that congruence survives addition and
multiplication. If $a \equiv b$ and $c \equiv d \pmod{m}$, then
$$a + c \equiv b + d \pmod{m}
\qquad \text{and} \qquad
ac \equiv bd \pmod{m}.$$
This means we can reduce intermediate results *before* multiplying,
keeping numbers small throughout a computation. Carl Friedrich Gauss
introduced the congruence notation $\equiv$ in his 1801
*Disquisitiones Arithmeticae* [@gauss1801] — the work that first organized number
theory into a coherent mathematical system.
```{python}
# Division algorithm: a = q*m + r
examples = [
(17, 5), (100, 12), (365, 7), (-7, 3)
]
for a, m in examples:
q, r = divmod(a, m)
print(f"{a:5d} = {q:4d} * {m} + {r}")
```
The last row shows an important convention: Python's `divmod` (and the
`%` operator) always returns a *non-negative* remainder when $m > 0$.
So $-7 = (-3) \cdot 3 + 2$, not $(-2) \cdot 3 + (-1)$.
### Research Example: Do Integers Tile into Neat Bands? {.unnumbered .unlisted}
A visual check: every integer belongs to exactly one residue class, and the classes repeat with period $m$. We color integers 0–34 by their remainder mod 3, 5, and 7 to see whether that claim holds in three different moduli at once.
```{python}
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('default')
fig, axes = plt.subplots(3, 1, figsize=(8, 3),
sharex=True)
n_vals = list(range(35))
for ax, m in zip(axes, [3, 5, 7]):
ax.scatter(
n_vals, [0] * 35,
c=[n % m for n in n_vals],
cmap='tab10', vmin=0, vmax=9, s=60
)
ax.set_yticks([])
ax.set_ylabel(f"mod {m}", fontsize=9)
axes[0].set_title(
"Integers 0-34 colored by residue class"
)
axes[-1].set_xlabel("n")
fig.tight_layout()
plt.show()
```
Each row repeats with period equal to its modulus, and the colors confirm that every integer falls into exactly one residue class — no overlap, no gaps. You just proved the Quotient-Remainder Theorem with a scatter plot.
## The `%` Operator: Python's Remainder {#sec-modular-operator}
The expression `a % m` returns the remainder of `a` divided by `m`.
For repeated modular multiplication, Python provides a three-argument
form: `pow(a, k, m)` computes $a^k \bmod m$ via **fast modular
exponentiation**. At each squaring step it reduces modulo $m$,
so intermediate values never exceed $m^2$.
Without this shortcut, computing $3^{200000} \bmod (10^9+7)$ would
first build a roughly 95,000-digit integer before taking the remainder.
With `pow(3, 200000, 10**9 + 7)` the calculation finishes in
microseconds.
```{python}
import time
a, k, m = 3, 200_000, 10**9 + 7
t0 = time.perf_counter()
fast = pow(a, k, m)
t1 = time.perf_counter()
slow = (a ** k) % m # builds a 95k-digit int first
t2 = time.perf_counter()
print(f"Result: {fast}")
print(f"Results match: {fast == slow}")
print(f"Fast pow (ms): {(t1 - t0) * 1000:.3f}")
print(f"Slow pow (ms): {(t2 - t1) * 1000:.3f}")
```
Both give the same answer. The speedup comes entirely from when the
modular reduction happens: `pow(a, k, m)` reduces at each squaring
step while `a ** k` builds the full power first.
The three-argument `pow` will reappear in @sec-modular-fermat when
we verify Fermat's Little Theorem across many primes.
## Divisibility Tests {#sec-modular-divisibility}
Grade-school rules like "a number is divisible by 9 if its digit sum
is divisible by 9" are theorems in modular arithmetic, not lucky
observations. The key is the base: in base 10, each digit is multiplied
by a power of 10, and those powers behave differently modulo different
primes.
**Divisibility by 2 and 5.** Since $10 \equiv 0 \pmod 2$ and
$10 \equiv 0 \pmod 5$, every higher power of 10 also vanishes. A
number's remainder mod 2 or mod 5 equals the remainder of its last
digit alone.
**Divisibility by 3 and 9.** Since $10 \equiv 1 \pmod 9$ (and
$\pmod 3$), every power $10^k \equiv 1$. For a number with digits
$d_k \cdots d_0$:
$$N = d_k \cdot 10^k + \cdots + d_0
\equiv d_k + \cdots + d_0 \pmod{9}.$$
$N$ is divisible by 9 (or 3) if and only if its **digit sum** is.
**Divisibility by 11.** Since $10 \equiv -1 \pmod{11}$, we get
$10^k \equiv (-1)^k$, and the digits alternate in sign:
$$N \equiv d_0 - d_1 + d_2 - d_3 + \cdots \pmod{11}.$$
$N$ is divisible by 11 if and only if its **alternating digit sum**
(starting from the rightmost digit with a $+$ sign) is divisible by
11.
```{python}
def digit_sum(n):
return sum(int(d) for d in str(abs(n)))
def alt_digit_sum(n):
digits = [int(d) for d in str(abs(n))]
return sum(
d if i % 2 == 0 else -d
for i, d in enumerate(reversed(digits))
)
for n in [123456789, 9801, 12345678901234]:
ds = digit_sum(n)
ads = alt_digit_sum(n)
print(f"n = {n}")
print(f" digit sum {ds:3d}: "
f"div by 9? {n%9==0}, "
f"rule: {ds%9==0}")
print(f" alt sum {ads:4d}: "
f"div by 11? {n%11==0}, "
f"rule: {ads%11==0}")
print()
```
Let us confirm the rules never fail on 10,000 random numbers:
```{python}
# uses: digit_sum(), alt_digit_sum()
import random
random.seed(42)
samples = [random.randint(1, 10**15)
for _ in range(10_000)]
errors = sum(
1 for n in samples
if ((n % 3 == 0) != (digit_sum(n) % 3 == 0)
or (n % 9 == 0) != (digit_sum(n) % 9 == 0)
or (n % 11 == 0) != (alt_digit_sum(n) % 11 == 0))
)
print(f"Errors in 10,000 trials: {errors}")
```
## Check Digits: ISBN, Barcodes, Credit Cards {#sec-modular-check}
```{python}
#| echo: false
from pathlib import Path; import urllib.request
_d = Path('images'); _d.mkdir(exist_ok=True)
_p = _d / 'shannon.jpg'
if not _p.exists():
try: urllib.request.urlretrieve('https://upload.wikimedia.org/wikipedia/commons/9/99/ClaudeShannon_MFO3807.jpg', _p)
except Exception: pass
```
::: {.content-visible when-format="pdf"}
```{=latex}
\begin{center}
\begin{minipage}[c]{0.22\textwidth}
\includegraphics[width=\textwidth]{images/shannon.jpg}
\end{minipage}%
\hspace{0.03\textwidth}%
\begin{minipage}[c]{0.35\textwidth}
\small\textit{Claude Shannon (1916--2001)}\\[2pt]
\tiny MFO, CC BY-SA 2.0
\end{minipage}
\end{center}
```
:::
::: {.content-visible when-format="html"}
<div style="display:flex; align-items:center; margin:1em 0; gap:12px;">
<img src="images/shannon.jpg" style="width:100px; flex-shrink:0;" alt="Claude Shannon">
<div style="font-size:0.82em;"><em>Claude Shannon (1916–2001)</em><br><span style="font-size:0.85em;">MFO, CC BY-SA 2.0</span></div>
</div>
:::
A **check digit** is an extra digit appended to a code so the entire
string satisfies a modular constraint. A single transcription error
shifts the modular sum by a nonzero amount, and if the modulus is
chosen carefully that shift cannot bring the sum back to zero ---
so the error is always caught.
### ISBN-10
A ten-digit ISBN (used until 2006) has digits $d_1 d_2 \cdots d_{10}$
satisfying the weighted sum:
$$\sum_{i=1}^{10} i \cdot d_i \equiv 0 \pmod{11}$$
Given the first nine digits, the check digit $d_{10}$ is the unique
value in $\{0,\ldots,10\}$ making the sum divisible by 11. (When that
value is 10 the letter X is printed.) The prime modulus 11 guarantees
detection of every single-digit error and every transposition of two
adjacent digits.
<!-- TODO companion-03: Algebraic derivation showing that the ISBN-10
weighted sum constraint 10*d10 ≡ -S9 (mod 11) simplifies to
d10 ≡ S9 (mod 11) because 10 ≡ -1 (mod 11). Full derivation belongs in
"The Mathematics Behind the Tour." -->
```{python}
def isbn10_check(first9):
"""
Compute ISBN-10 check digit from first 9 digits.
Returns 0-10 (10 means 'X').
"""
s = sum((i + 1) * d
for i, d in enumerate(first9))
return s % 11
def isbn10_valid(s):
"""Validate a 10-character ISBN string."""
s = s.replace('-', '')
if len(s) != 10:
return False
digits = []
for ch in s:
if ch == 'X':
digits.append(10)
elif ch.isdigit():
digits.append(int(ch))
else:
return False
total = sum(
(i + 1) * d for i, d in enumerate(digits)
)
return total % 11 == 0
# 0-306-40615-?: compute check digit
first9 = [0, 3, 0, 6, 4, 0, 6, 1, 5]
cd = isbn10_check(first9)
label = str(cd) if cd < 10 else 'X'
print(f"Check digit: {label}")
tests = [
("0-306-40615-2", True),
("0-306-40615-3", False),
("0-201-53082-1", True),
]
for isbn, expected in tests:
ok = isbn10_valid(isbn)
tag = "OK" if ok == expected else "FAIL"
print(f"{isbn}: valid={ok} [{tag}]")
```
### Research Example: What Does a Book's Check Digit Actually Weigh? {.unnumbered .unlisted}
The chart below makes the weighting scheme visible for ISBN 0-306-40615-2. Each bar shows one digit multiplied by its position weight; the total, read mod 11, must be zero.
```{python}
# uses: isbn10_check()
import matplotlib.pyplot as plt
isbn_str = '0306406152'
digits = [10 if c == 'X' else int(c) for c in isbn_str]
weights = list(range(1, 11))
products = [w * d for w, d in zip(weights, digits)]
fig, ax = plt.subplots(figsize=(9, 3.5))
bars = ax.bar(range(1, 11), products,
color='#5B9BD5', edgecolor='white', linewidth=0.5)
for pos, (d, w, p) in enumerate(zip(digits, weights, products), start=1):
label = 'X' if d == 10 else str(d)
ax.text(pos, p + 0.3, f'd={label}', ha='center', va='bottom',
fontsize=8, color='#333333')
total = sum(products)
ax.axhline(total / 10, color='#ED7D31', lw=1.2, ls='--')
ax.text(10.4, total / 10,
f'avg\n({total} mod 11 = {total % 11})',
va='center', fontsize=8, color='#ED7D31')
ax.set_xlabel('position (weight)')
ax.set_ylabel('weight x digit')
ax.set_xticks(range(1, 11))
ax.set_title('ISBN-10 weighted sum: 0-306-40615-2', fontsize=11)
fig.tight_layout()
plt.show()
```
The bars reveal that the check constraint is a carefully balanced weighted sum — not an afterthought, but a design that guarantees any single-digit typo shifts the total by a nonzero amount mod 11. You just replicated, in 20 lines of Python, the calculation every library scanner performs in milliseconds.
### ISBN-13 and EAN-13
Books published since 2007 carry ISBN-13, identical to the EAN-13
barcode standard on retail products. The 13 digits satisfy:
$$d_1 + 3d_2 + d_3 + 3d_4 + \cdots + 3d_{12} + d_{13}
\equiv 0 \pmod{10}$$
Weights alternate between 1 and 3. The modulus 10 is convenient (the
check digit is a single decimal digit) but weaker than 11: some
adjacent transpositions are undetected, as a research topic below
explores.
```{python}
def ean13_check(first12):
"""Compute EAN-13 / ISBN-13 check digit."""
weights = [1, 3] * 6
s = sum(w * d for w, d in zip(weights, first12))
return (10 - s % 10) % 10
def ean13_valid(s):
"""Validate a 13-digit EAN / ISBN string."""
s = s.replace('-', '')
if len(s) != 13:
return False
digits = [int(c) for c in s]
weights = [1, 3] * 6 + [1]
return (sum(w * d for w, d in
zip(weights, digits)) % 10 == 0)
for s in [
"978-0-306-40615-7", # valid
"978-0-306-40615-8", # wrong check digit
]:
print(f"{s}: {ean13_valid(s)}")
```
### The Luhn Algorithm
Credit card numbers use the **Luhn algorithm** (modulus 10) to
validate a card number before a transaction reaches the bank:
1. From the second-to-last digit moving left, double every other digit.
2. If doubling yields $\geq 10$, subtract 9.
3. Sum all digits. A valid card gives a sum divisible by 10.
```{python}
def luhn_valid(s):
"""Validate a credit card number string."""
digits = [int(c) for c in s if c.isdigit()]
total = 0
for i, d in enumerate(reversed(digits)):
if i % 2 == 1:
d *= 2
if d >= 10:
d -= 9
total += d
return total % 10 == 0
cards = [
("4532015112830366", True),
("4532015112830367", False),
("79927398713", True),
("79927398714", False),
]
for card, expected in cards:
ok = luhn_valid(card)
tag = "OK" if ok == expected else "FAIL"
print(f"{card}: valid={ok} [{tag}]")
```
### Research Example: Can One Diagram Trace a Credit Card's Entire Check? {.unnumbered .unlisted}
The diagram below traces the algorithm step by step on a single card number, making the doubling and subtraction concrete. Can 16 digits fit in one picture that shows exactly which digits get doubled, which get reduced, and how the total lands on a multiple of 10?
```{python}
#| label: fig-luhn-steps
#| fig-cap: "Luhn algorithm applied to card 4532 0151 1283 0366. Shaded (blue) columns are doubled; the bottom row shows each digit's contribution to the total. Total = 70, and $70 \\bmod 10 = 0$, so the card is VALID."
#| fig-align: center
#| out-width: "100%"
import matplotlib.pyplot as plt
import matplotlib.patches as mpatch
BLUE = '#2266AA'
def luhn_table(s):
digits = [int(c) for c in s if c.isdigit()]
rows = []
for i, d in enumerate(reversed(digits)):
if i % 2 == 1:
doubled = d * 2
contrib = doubled - 9 if doubled >= 10 else doubled
rows.append((d, True, doubled, contrib))
else:
rows.append((d, False, d, d))
rows.reverse()
return rows
card = '4532015112830366'
steps = luhn_table(card)
total = sum(r[3] for r in steps)
fig, ax = plt.subplots(figsize=(13, 3.2))
ax.axis('off')
for i, (orig, dbl, dv, contrib) in enumerate(steps):
bg = '#CCDDEE' if dbl else '#F5F5F5'
ax.add_patch(mpatch.Rectangle(
(i * 0.9, 0.6), 0.82, 1.6,
color=bg, ec='#AAAAAA', lw=0.5, zorder=1))
ax.text(i * 0.9 + 0.41, 1.9, str(orig),
ha='center', va='center', fontsize=13, fontweight='bold')
if dbl:
note = f'x2={dv}' + ('-9' if dv >= 10 else '')
ax.text(i * 0.9 + 0.41, 1.35, note,
ha='center', va='center', fontsize=7, color=BLUE)
ax.text(i * 0.9 + 0.41, 0.85, str(contrib),
ha='center', va='center', fontsize=13, fontweight='bold',
color=BLUE if dbl else '#333333')
ax.text(-0.6, 1.9, 'digit:', ha='right', va='center', fontsize=9, color='#555555')
ax.text(-0.6, 0.85, 'contrib:', ha='right', va='center', fontsize=9, color='#555555')
ax.text(len(steps) * 0.9 / 2, 0.2,
f'Total = {total} {total} mod 10 = {total % 10} -> VALID',
ha='center', va='center', fontsize=11, color='#003377', fontweight='bold')
ax.set_xlim(-1.2, len(steps) * 0.9 + 0.3)
ax.set_ylim(0.0, 2.4)
ax.set_title(
f'Luhn algorithm: {card[:4]} {card[4:8]} {card[8:12]} {card[12:]}',
fontsize=11, pad=6)
fig.tight_layout()
plt.show()
```
Every blue column reveals a doubled digit; the orange note shows where 9 is subtracted to keep the value single-digit. The total 70 is divisible by 10, confirming the card is valid — the same check your phone completes in a fraction of a second before a tap payment goes through.
```{python}
#| echo: false
from pathlib import Path
import urllib.request
_d = Path('images'); _d.mkdir(exist_ok=True)
_p = _d / 'thumb_kMpEj7roKxM.jpg'
try:
if not _p.exists():
urllib.request.urlretrieve('https://img.youtube.com/vi/kMpEj7roKxM/0.jpg', _p)
except Exception:
pass
```
::: {.content-visible when-format="pdf"}
```{=latex}
\begin{center}
\begin{minipage}[c]{0.28\textwidth}
\centering
\href{https://youtu.be/kMpEj7roKxM}{\includegraphics[width=\textwidth]{images/thumb_kMpEj7roKxM.jpg}}
\end{minipage}%
\hspace{0.02\textwidth}%
\begin{minipage}[c]{0.28\textwidth}
\small\textbf{Barry Brown}\\[3pt]
\small Checksums: The Luhn Algorithm for Verifying Credit Card Numbers\\[3pt]
\small\href{https://youtu.be/kMpEj7roKxM}{\texttt{youtu.be/kMpEj7roKxM}}
\end{minipage}%
\hspace{0.02\textwidth}%
\begin{minipage}[c]{0.36\textwidth}
\small Walks through the Luhn algorithm step by step, showing how credit card check digits catch common transcription errors.
\end{minipage}
\end{center}
```
:::
::: {.content-visible when-format="html"}
<div style="display:flex; align-items:flex-start; margin:1em 0; gap:12px; width:100%;">
<div style="flex:0 0 200px;"><a href="https://youtu.be/kMpEj7roKxM" target="_blank"><img src="https://img.youtube.com/vi/kMpEj7roKxM/0.jpg" style="width:100%;display:block;" alt="Checksums: The Luhn Algorithm for Verifying Credit Card Numbers"></a></div>
<div style="flex:1; font-size:0.85em;"><strong>Barry Brown</strong><br>Checksums: The Luhn Algorithm for Verifying Credit Card Numbers<br><a href="https://youtu.be/kMpEj7roKxM" target="_blank" style="font-family:Consolas,monospace;">youtu.be/kMpEj7roKxM</a></div>
<div style="flex:1; font-size:0.85em;">Walks through the Luhn algorithm step by step, showing how credit card check digits catch common transcription errors.</div>
</div>
:::
All three systems share the same design: encode a modular constraint
in an extra digit so any single transcription error is detected
automatically.
## Fermat's Little Theorem {#sec-modular-fermat}
```{python}
#| echo: false
from pathlib import Path; import urllib.request
_d = Path('images'); _d.mkdir(exist_ok=True)
_p = _d / 'fermat.jpg'
if not _p.exists():
try: urllib.request.urlretrieve('https://upload.wikimedia.org/wikipedia/commons/f/f3/Pierre_de_Fermat.jpg', _p)
except Exception: pass
```
::: {.content-visible when-format="pdf"}
```{=latex}
\begin{center}
\begin{minipage}[c]{0.22\textwidth}
\includegraphics[width=\textwidth]{images/fermat.jpg}
\end{minipage}%
\hspace{0.03\textwidth}%
\begin{minipage}[c]{0.35\textwidth}
\small\textit{Pierre de Fermat (1601--1665)}\\[2pt]
\tiny Public domain
\end{minipage}
\end{center}
```
:::
::: {.content-visible when-format="html"}
<div style="display:flex; align-items:center; margin:1em 0; gap:12px;">
<img src="images/fermat.jpg" style="width:100px; flex-shrink:0;" alt="Pierre de Fermat">
<div style="font-size:0.82em;"><em>Pierre de Fermat (1601–1665)</em><br><span style="font-size:0.85em;">Public domain</span></div>
</div>
:::
Let $p$ be prime and let $a$ be any integer with $\gcd(a, p) = 1$.
Then:
$$a^{p-1} \equiv 1 \pmod{p}.$$
This is **Fermat's Little Theorem**, stated without proof in 1640 (Euler supplied the first proof in 1736). A more general
form states $a^p \equiv a \pmod p$ for *every* integer $a$ (including
multiples of $p$, where both sides equal 0).
<!-- TODO companion-03: Full proof of Fermat's Little Theorem via the
residue-shuffling argument: the set {1,2,...,p-1} is closed under
multiplication mod p, multiplying by a permutes it, comparing the two
products forces a^(p-1) ≡ 1. Belongs in "The Mathematics Behind the Tour." -->
The three-argument `pow` from @sec-modular-operator makes the
verification instant. `isprime` and `nextprime` were introduced in
@sec-primes-what.
```{python}
from sympy import nextprime
# Start at p=7 so a=1..5 are all coprime to p
p = 7
print("p a=1 a=2 a=3 a=4 a=5")
for _ in range(10):
vals = [pow(a, p - 1, p) for a in range(1, 6)]
print(f"{p:5d} "
+ " ".join(str(v) for v in vals))
p = nextprime(p)
```
Every entry is 1, as the theorem guarantees for primes $p > 5$ with
bases $a = 1, \ldots, 5$ (all coprime to $p$).
**The pitfall: Fermat pseudoprimes.** The converse fails. A composite
$n$ satisfying $a^{n-1} \equiv 1 \pmod n$ is called a **Fermat
pseudoprime** in base $a$. Let us find them.
```{python}
from sympy import isprime
pseudoprimes = [
n for n in range(3, 10_000, 2)
if not isprime(n)
and pow(2, n - 1, n) == 1
]
print("Pseudoprimes (base 2) below 10,000:")
print(pseudoprimes)
```
The **Carmichael numbers** --- 561, 1105, 1729, and others in this
list --- are composites that fool the Fermat test for *every* base
$a$ coprime to $n$. The number 1729 appears here and is also the
smallest "taxicab number": $1^3 + 12^3 = 10^3 + 9^3 = 1729$.
```{python}
#| echo: false
from pathlib import Path
import urllib.request
_d = Path('images'); _d.mkdir(exist_ok=True)
_p = _d / 'thumb_S9JGmA5_unY.jpg'
try:
if not _p.exists():
urllib.request.urlretrieve('https://img.youtube.com/vi/S9JGmA5_unY/0.jpg', _p)
except Exception:
pass
```
::: {.content-visible when-format="pdf"}
```{=latex}
\begin{center}
\begin{minipage}[c]{0.28\textwidth}
\centering
\href{https://youtu.be/S9JGmA5_unY}{\includegraphics[width=\textwidth]{images/thumb_S9JGmA5_unY.jpg}}
\end{minipage}%
\hspace{0.02\textwidth}%
\begin{minipage}[c]{0.28\textwidth}
\small\textbf{3Blue1Brown}\\[3pt]
\small How secure is 256 bit security?\\[3pt]
\small\href{https://youtu.be/S9JGmA5_unY}{\texttt{youtu.be/S9JGmA5\_unY}}
\end{minipage}%
\hspace{0.02\textwidth}%
\begin{minipage}[c]{0.36\textwidth}
\small Visualizes why a 256-bit key space is so astronomically large that exhaustive search is physically impossible even with all the energy in the solar system.
\end{minipage}
\end{center}
```
:::
::: {.content-visible when-format="html"}
<div style="display:flex; align-items:flex-start; margin:1em 0; gap:12px; width:100%;">
<div style="flex:0 0 200px;"><a href="https://youtu.be/S9JGmA5_unY" target="_blank"><img src="https://img.youtube.com/vi/S9JGmA5_unY/0.jpg" style="width:100%;display:block;" alt="How secure is 256 bit security?"></a></div>
<div style="flex:1; font-size:0.85em;"><strong>3Blue1Brown</strong><br>How secure is 256 bit security?<br><a href="https://youtu.be/S9JGmA5_unY" target="_blank" style="font-family:Consolas,monospace;">youtu.be/S9JGmA5_unY</a></div>
<div style="flex:1; font-size:0.85em;">Visualizes why a 256-bit key space is so astronomically large that exhaustive search is physically impossible even with all the energy in the solar system.</div>
</div>
:::
Despite the caveat, Fermat's theorem is the seed from which modern
cryptography grows. RSA encryption [@rivest1978] encodes a message $m$ as
$c = m^e \bmod n$ and recovers it using Fermat's theorem to find an
inverse exponent. Stronger primality tests --- Miller-Rabin and the
AKS algorithm [@agrawal2002] --- close the Carmichael loophole, but both trace their
logic back to this 1640 theorem.
## Further Research Topics {#sec-modular-research}
The following topics are listed roughly in order of increasing
difficulty. Begin by experimenting with small cases, look for patterns,
form a conjecture, then try to explain what you find.
**1. Perpetual calendar.**
If January 1, 2000 was a Saturday, then any future (or past) date can
be reached by counting total days elapsed and taking that count mod 7.
Write a function `day_of_week(year, month, day)` that returns the
weekday using only integer arithmetic and `%` --- no importing
`datetime`. Test it against `datetime.date(year, month, day).weekday()`
for at least 50 dates spanning several centuries. Hint: you need to
handle leap years (divisible by 4 but not 100, or divisible by 400).
*(Problem proposed by Claude Code.)*
**2. Casting out nines.**
Before calculators, "casting out nines" let students spot arithmetic
errors: compute digit sums of all operands and the claimed answer, then
check that the digit sums satisfy the same arithmetic relation mod 9.
For example, verify $247 \times 389 = 96083$ by checking
$(4)(2) \equiv 8 \pmod 9$ and that the digit sum of 96083 is $\equiv 8
\pmod 9$. Find an example where casting out nines produces a false
positive (the product is wrong but the digit-sum check passes anyway).
How often does this happen? Estimate the probability that a random
single-digit error in a product goes undetected by casting out nines.
*(Problem proposed by Claude Code.)*
**3. Catching every ISBN-10 error.**
Prove algebraically that the ISBN-10 scheme (weighted sum mod 11)
detects every single-digit substitution error. Then prove it detects
every adjacent transposition. Verify both claims computationally: for
"0-306-40615-2", generate every possible single-digit modification and
every adjacent swap, and confirm that each one fails `isbn10_valid`.
*(Problem proposed by Claude Code.)*
**4. Where EAN-13 fails.**
EAN-13 weights are 1 and 3 alternating, with modulus 10. Find which
pairs of adjacent digits $(d_i, d_{i+1})$ can be transposed without
changing the check sum. Explain why the choice of modulus 10 (instead
of a prime like 11) is the root cause. How many undetectable adjacent
swaps exist in a random 13-digit EAN code? How would you fix the
scheme if you could change the weights or the modulus?
*(Problem proposed by Claude Code.)*
**5. Solving linear congruences.**
A linear congruence $ax \equiv b \pmod m$ has a solution if and only
if $\gcd(a, m)$ divides $b$; when solutions exist, there are exactly
$\gcd(a, m)$ distinct solutions in $\{0, \ldots, m-1\}$. Use
`math.gcd` and a brute-force search to find all solutions to
$6x \equiv 4 \pmod{14}$ and $7x \equiv 3 \pmod{15}$. Then implement
the extended Euclidean algorithm to solve the general case without
brute force.
*(Problem proposed by Claude Code.)*
**6. Repeating decimals and multiplicative order.**
The decimal expansion of $1/7 = 0.\overline{142857}$ repeats with
period 6. This period is the smallest positive $k$ such that
$10^k \equiv 1 \pmod 7$ --- called the *multiplicative order* of 10
mod 7. Write Python to compute the multiplicative order of 10 mod $p$
for every prime $p < 200$ (skip 2 and 5). Primes whose order equals
$p - 1$ are called *full reptend primes*: they produce the longest
possible repeating blocks. List all full reptend primes below 200 and
verify each period by converting `1/p` with Python's `fractions`
module. Conjecture: 10 is a *primitive root* mod $p$ if and only if
$p$ is a full reptend prime. Finally, show computationally that the
period of $1/p^2$ is always either the same as, or exactly $p$ times,
the period of $1/p$.
*(Problem proposed by Claude Code.)*
**7. Digit sums in other bases.**
The base-10 digit-sum rule for divisibility by 9 rests on
$10 \equiv 1 \pmod 9$. In base 16 (hexadecimal), which modulus has
an analogous sum-of-hex-digits rule? More generally, for which pairs
$(b, m)$ does a base-$b$ digit-sum rule exist? Write Python code to
compute hexadecimal digit sums and verify the rule. State a general
criterion in terms of $b \bmod m$.
*(Problem proposed by Claude Code.)*
**8. Euler's totient function.**
Euler's totient $\varphi(n)$ counts how many integers in
$\{1, 2, \ldots, n\}$ share no common factor with $n$. Compute
$\varphi(n)$ for $n = 1$ to 50 and look for patterns. Discover:
$\varphi(p) = p - 1$ for prime $p$, and $\varphi(mn) = \varphi(m)\varphi(n)$
whenever $\gcd(m, n) = 1$ (multiplicativity). Use these facts to
derive $\varphi(p^k) = p^k - p^{k-1}$ and verify it in Python.
Then verify **Euler's theorem**: if $\gcd(a, n) = 1$ then
$a^{\varphi(n)} \equiv 1 \pmod n$. (Fermat's Little Theorem is the
special case $n = p$.) Plot $\varphi(n)/n$ against $n$ for
$n \leq 500$ and explain why the graph never drops to zero but comes
arbitrarily close. Finally, verify the identity
$\sum_{d \mid n} \varphi(d) = n$ for $n = 1$ to 50 and give a
combinatorial explanation.
*(Problem proposed by Claude Code.)*
**9. Fermat pseudoprimes and multiple bases.**
How many composites below $10^6$ pass the Fermat test for base
$a = 2$? For $a = 3$? For both simultaneously? Plot the fraction of
composites that are fooled as a function of the upper limit, for each
base separately and for both together. Can you find the smallest
composite that passes for $a = 2, 3, 5, 7$ simultaneously?
*(Problem proposed by Claude Code.)*
**10. Carmichael numbers: fools that fool everyone.**
The Fermat test rejects most composites --- but some composites pass
for *every* base $a$ with $\gcd(a, n) = 1$. These are **Carmichael
numbers**, named after Robert Carmichael, who studied them in 1910.
The smallest is $561 = 3 \times 11 \times 17$.
**Korselt's criterion** [@korselt1899]: $n > 1$ is a Carmichael number
if and only if $n$ is squarefree and $(p - 1) \mid (n - 1)$ for every
prime $p$ dividing $n$. Verify Korselt's criterion for 561 by hand,
then write Python to find all Carmichael numbers below $10^6$.
How does your count compare to the count of Fermat pseudoprimes
(base 2) from topic~9? Finally, look up the Miller-Rabin test and
explain in one paragraph why it cannot be fooled by Carmichael numbers,
even though it generalizes the Fermat test.
*(Problem proposed by Claude Code.)*
**11. Wilson's theorem.**
Wilson's theorem: $(p-1)! \equiv -1 \pmod p$ if and only if $p$ is
prime. Verify this for all primes up to 100 using `math.factorial`.
Test composites: how does $(n-1)! \bmod n$ behave for composite $n$?
Wilson's theorem gives a perfect mathematical characterization of
primality but is computationally useless for large $n$ --- explain
why, and estimate the size of $(n-1)!$ for $n = 100$.
*(Problem proposed by Claude Code.)*
**12. Quadratic residues and Euler's criterion.**
An integer $a$ is a **quadratic residue** mod prime $p$ if
$x^2 \equiv a \pmod p$ has a solution. Compute all quadratic residues
mod 13 by squaring $1, 2, \ldots, 12$ mod 13. How many distinct
values appear? Conjecture the count for general prime $p$ and prove
it by considering collisions in the squaring map. Then verify
**Euler's criterion**: $a^{(p-1)/2} \equiv 1 \pmod p$ if $a$ is a
quadratic residue, and $\equiv -1$ otherwise.
*(Problem proposed by Claude Code.)*
**13. Primitive roots and the discrete logarithm.**
An integer $g$ is a **primitive root** modulo prime $p$ if the powers
$g^1, g^2, \ldots, g^{p-1} \pmod p$ produce every nonzero residue
exactly once (Gauss proved that every prime has a primitive root).
Find all primitive roots mod 7, 11, and 13 by brute force. Then write
Python to check, for each prime $p < 500$, whether 2 is a primitive
root mod $p$. Compute the empirical fraction and compare it to the
constant 0.3739558\ldots --- Artin's conjecture (1927) predicts this is
the limiting density, but the conjecture remains unproven to this day
[@artin1927]. Next, implement the **baby-step giant-step** algorithm
(Shanks, 1971) to solve $g^x \equiv h \pmod p$ in $O(\sqrt{p})$ steps:
precompute $\{g^j \bmod p : 0 \leq j < \lceil\sqrt{p}\rceil\}$, then
for each $i$ test whether $h \cdot (g^{-\lceil\sqrt{p}\rceil})^i$ is
in the table [@shanks1971]. Verify on several examples that your solver
matches a brute-force search.
*(Problem proposed by Claude Code.)*
**14. The Chinese Remainder Theorem.**
If $m_1$ and $m_2$ are coprime, the system $x \equiv a_1 \pmod{m_1}$
and $x \equiv a_2 \pmod{m_2}$ has a unique solution mod $m_1 m_2$.
Implement a solver using the extended Euclidean algorithm and verify
it on several examples. Extend to three congruences. Explain why
$10 \equiv 3 \pmod 7$ and $10 \equiv 1 \pmod 3$ together imply
$10 \equiv 10 \pmod{21}$, and use the theorem to break a large
modular arithmetic problem into smaller pieces.
*(Problem proposed by Claude Code.)*
**15. The music of modular arithmetic.**
The 12 notes in Western music form a clock: $\text{C}=0$,
$\text{C}\sharp=1$, \ldots, $\text{B}=11$. A *perfect fifth* is a
jump of 7 semitones; starting on C and jumping by 7 mod 12 cycles
through all 12 notes before returning to C, because $\gcd(7,12)=1$.
Write Python to find all generators of $\mathbb{Z}_{12}$ (intervals
that visit every note exactly once before repeating). Count them and
verify the count equals $\varphi(12)$. Repeat for $\mathbb{Z}_n$ with
$n = 8, 10, 15, 24$ and confirm the count is always $\varphi(n)$.
The *tritone* (6 semitones) is not a generator: it visits only
$12/\gcd(6,12)=2$ distinct notes before cycling back. Use this
formula to predict the cycle length of any interval $k$ in an $n$-note
scale without computing the cycle.
*(Problem proposed by Claude Code.)*
**16. Pisano periods: Fibonacci numbers on a clock.**
The Fibonacci sequence mod $m$ must eventually repeat, since there are
only $m^2$ possible consecutive pairs $(F_n \bmod m, F_{n+1} \bmod m)$.
The repeat period is the **Pisano period** $\pi(m)$: the smallest
$k>0$ with $F_k \equiv 0$ and $F_{k+1} \equiv 1 \pmod m$. Compute
$\pi(m)$ for $m=2,3,\ldots,50$. Plot the results and observe that
$\pi(2)=3$, $\pi(5)=20$, $\pi(10)=60$. Verify that $\pi(mn)$
divides $\mathrm{lcm}(\pi(m),\pi(n))$ whenever $\gcd(m,n)=1$. Make
a scatter plot of $\pi(p)$ vs. prime $p < 300$: what fraction of
primes satisfy $\pi(p) \mid p-1$? What fraction satisfy
$\pi(p) \mid 2(p+1)$?
*(Problem proposed by Claude Code.)*
**17. Burnside's lemma: counting necklaces.**
How many distinct binary necklaces of length $n$ exist, where two
necklaces are the same if one is a rotation of the other? Burnside's
lemma says the answer is $\frac{1}{n}\sum_{k=0}^{n-1}2^{\gcd(k,n)}$
(because a rotation by $k$ fixes exactly $2^{\gcd(k,n)}$ strings).
Implement this formula and compute the necklace count for $n=1$ to 20.
Verify by brute force for $n \leq 12$: generate all binary strings,
map each to its lexicographically smallest rotation, and count distinct
images. Extend to $q$-color necklaces (replace $2^{\gcd(k,n)}$ with
$q^{\gcd(k,n)}$) and produce a table for $q=2,3,4$ and $n=1$ to 15.
Notice that the modular sum $\sum \varphi(d) \cdot q^{n/d}$ (summed
over divisors $d$ of $n$), divided by $n$, gives the same count:
verify this identity and explain why it is equivalent to Burnside.
*(Problem proposed by Claude Code.)*
**18. Miller-Rabin primality test.**
The Fermat test fails for Carmichael numbers (topic 10). Miller-Rabin
closes the gap. Write $n-1 = 2^s \cdot d$ with $d$ odd. A base $a$
is a *witness* to the compositeness of $n$ if $a^d \not\equiv 1$ and
$a^{2^r d} \not\equiv -1 \pmod n$ for all $0 \leq r < s$.
Implement `miller_rabin(n, a)` that returns `True` if $a$ is such a
witness. Test every Carmichael number below $10^6$ from topic 10: how
many fail for $a=2$? For $a\in\{2,3\}$? For $a\in\{2,3,5,7\}$? Plot
the count of survivors as a function of the base set size. The first
composite that fools $\{2,3,5,7,11,13,17,19,23,29,31,37\}$ is
enormously large --- use Python to verify that every number below
$3.3\times10^{24}$ is correctly classified by this set of bases.
*(Problem proposed by Claude Code.)*
**19. Modular square roots and the Tonelli-Shanks algorithm.**
Euler's criterion tells you *whether* $a$ has a square root mod prime
$p$ (it does iff $a^{(p-1)/2} \equiv 1$), but not *what* the root is.
For $p \equiv 3 \pmod 4$, the root is simply $a^{(p+1)/4} \bmod p$ ---
prove why using Euler's criterion. For other primes the
**Tonelli-Shanks algorithm** (1891) efficiently finds roots. Implement
both: use the simple formula when $p \equiv 3 \pmod 4$, and code the
full Tonelli-Shanks loop otherwise. Test on $a=5$, $p=41$ (answer: 28
since $28^2 = 784 = 19\cdot 41+5$), and on a prime $p > 10^{12}$.
Compare your output to SymPy's `sqrt_mod(a, p)`. Finally, count how
many of the primes $p < 1000$ require the full algorithm (i.e., have
$p \not\equiv 3 \pmod 4$) and plot the distribution of the two cases.
*(Problem proposed by Claude Code.)*
**20. Power towers mod $p$ and eventual convergence.**
Define the sequence $T_1 = a$, $T_2 = a^a \bmod p$,
$T_3 = a^{a^a} \bmod p$, \ldots (each exponent computed exactly before
reducing the final result mod $p$). For $a=2$, $p=13$ the sequence
starts $2, 4, 3, 1, 2, \ldots$ and cycles with period 4. Write Python
to compute $T_k(a,p) \bmod p$ for $k=1,\ldots,30$ and find the cycle
length for each $(a,p)$ pair with prime $p < 50$ and $1 \leq a < p$.
Produce a heat map: $a$ on one axis, $p$ on the other, color = first
$k$ where the sequence stabilizes. It is a theorem that the infinite
power tower $a^{a^{a^{\cdots}}}$ converges mod $p$ for any $a \geq 2$,
because the exponent sequence eventually has period dividing
$\lambda(p) = p - 1$. Verify this convergence for five choices of
$(a, p)$ and explain informally why the tower must settle down.
*(Problem proposed by Claude Code.)*
