🖉 A stupid "real-life" application of quadratic reciprocity

The application

During this year’s MOP, we used the following procedure to divide some of our students into two classes:

Let $p = 7075374838595186541578161$ be prime. Take the letters in your name as it appears on the roster, convert them with A1Z26 and take the sum of cubes to get a number $s$ . For example, EVANCHEN corresponds to $s = 5^3 + 22^3 + \dots + 14^3 = 16926$ . Then you’re in Red 1 (room A155) if $s$ is a quadratic residue modulo $p$ , and Red 2 (room A133) otherwise.

The students were understandably a bit confused why the prime was chosen. It turned out to be a prank: if you ran the calculation on the 30-ish students in this class, it was actually just splitting them by last name. (That is, last names A-P in one class and R-Z in one class).

A few people have asked me how I managed to find a prime with this property, so I thought I’d post the source code here (without the names).

How I did it

After computing $s$ for each student, you first start by taking a prime factorization of each $s$ value. Let $Q$ denote the set of primes dividing some $s$ . Then for each $q \in Q$ you want to commit to a value of the Legendre symbol $(\frac qp) = \pm 1$ . This is an $\mathbb{F}_2$ system of equations in $|Q|$ variables, which one can solve easily, e.g. with z3. (Using sums of cubes was necessary to get enough different primes to show up. If you used say sum of squares, for this year’s students, it turned out there was no solution to the $\mathbb{F}_2$ system, meaning no such prime $p$ could be found.)

Then after choosing a solution $Q \to \{\pm 1\}$ , quadratic reciprocity turns it into a condition on $p \pmod q$ . To make things simpler I restricted to $p \equiv 1 \pmod 4$ . The simplest thing to do would then be to pick a single residue class for each $q$ that has the correct Legendre symbol, use the Chinese remainder theorem to put them all together, then search for a prime in that residue class.

However, if you do this the prime you get will have size like $\prod_{q \in Q} q$ , which I thought was too large — it was something like 80 or 100 digits this way.

So instead, I arbitrarily picked a threshold on $Q$ where I only used Chinese remainder theorem on the $q < 100$ , and then just searched among that residue class for any prime that also passed a coin flip for each $q > 100$ . That brought the size of the prime down to the $25$ -digit one you see here, with only a minute or so of thinking in Python.

Code

Here’s the source code:

import string
from collections import defaultdict

from sympy import isprime, jacobi_symbol, legendre_symbol
from sympy.ntheory.modular import crt
from z3 import Bool, PbEq, Solver, sat

LARGE_PRIME_THRESHOLD = 100

with open("rnames.txt") as f:
    lines = [line.strip() for line in f if line.strip()]

sorted_lines = sorted(lines, key=lambda name: name.split()[1])

n = 0
nsums = []
for line in sorted_lines:
    s = sum(
        (string.ascii_uppercase.index(c) + 1) ** 3
        for c in line.upper()
        if c in string.ascii_uppercase
    )
    nsums.append(s)
    print(line, s)

A = nsums[:14]
B = nsums[14:]
print(A, B)


# Factorization helper
def prime_factors(n):
    i = 2
    factors = defaultdict(int)
    while i * i <= n:
        while n % i == 0:
            n //= i
            factors[i] += 1
        i += 1
    if n > 1:
        factors[n] += 1
    return dict(factors)


def solve_completely_multiplicative_z3(A, B):
    all_nums = A + B
    primes = set()
    factorizations = {}

    # Get all primes involved and factor each number
    for n in all_nums:
        f = prime_factors(n)
        factorizations[n] = f
        primes.update(f)

    primes = sorted(primes)
    f_vars = {p: Bool(f"f_{p}") for p in primes}  # True = +1, False = -1

    solver = Solver()

    def parity_constraint(factors, desired_sign):
        # We only care about odd exponents, since (-1)^even = 1
        relevant = [f_vars[p] for p, e in factors.items() if e % 2 == 1]
        if not relevant:
            # f(n) = 1 for empty product, so desired_sign must be +1
            return desired_sign == 1
        # Count how many of the relevant primes are assigned -1 (i.e., f_p = False)
        # Count number of False values → parity must match desired_sign
        # desired_sign = 1 ⇒ even number of Falses; -1 ⇒ odd number
        k = len(relevant)
        parity = 0 if desired_sign == 1 else 1
        solver.add(PbEq([(v, 1) for v in relevant], k - parity))

    for a in A:
        parity_constraint(factorizations[a], +1)
    for b in B:
        parity_constraint(factorizations[b], -1)

    if solver.check() == sat:
        model = solver.model()
        assignment = {p: 1 if model.eval(f_vars[p]) else -1 for p in primes}
        return assignment
    else:
        return None


result = solve_completely_multiplicative_z3(A, B)
assert result is not None
print(result)
print(len(result))

residues = []
moduli = []

if 2 in result:
    moduli.append(8)
    if result.pop(2) == 1:
        residues.append(1)
    else:
        residues.append(5)
for q in result:
    if q > LARGE_PRIME_THRESHOLD:
        continue
    moduli.append(q)
    if result[q] == 1:
        residues.append(1)
    else:
        for t in range(2, q - 1):
            if jacobi_symbol(t, q) == -1:
                residues.append(t)
                break
print(residues, moduli)
large_primes = [q for q in result.keys() if q > LARGE_PRIME_THRESHOLD]
r, M = crt(moduli, residues)
print(r, M)
print(f"{len(large_primes)} large prime conditions to deal with…")

for p in range(r, r + 10**8 * M, M):
    if not isprime(p):
        continue
    if all(legendre_symbol(q, p) == result[q] for q in large_primes):
        print(p)
        break
else:
    raise Exception("too bad so sad")


with open("prime.txt", "w") as f:
    print(p, file=f)

posted at 13:37 · math number theory · Math

Power Overwhelming

The blog of Evan Chen

Jul 11, 2025

🖉 A stupid "real-life" application of quadratic reciprocity

The application

How I did it

Code