Cipolla's algorithm in Python
Computes a square root modulo an odd prime by working in a quadratic extension field.
- File
cipolla_algorithm.py, 67 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/cipolla-algorithm/. The authors are listed in prov/W39G65.jsonl.
"""Cipolla's algorithm: a square root modulo an odd prime p.
For a quadratic residue n modulo p, choose a at random until d = a^2 - n is a
quadratic nonresidue. In the field F_p(w) with w^2 = d, the element
(a + w)^((p+1)/2) lies in F_p and is a square root of n. Elements x + y w are
stored as pairs (x, y) of residues in [0, p).
The modulus must be prime. The function does not test primality, but it never
returns a wrong root: for a composite modulus it raises ValueError when it
notices, or returns None, or returns a value whose square is n.
"""
import random
from typing import Optional
MAX_TRIALS = 200 # for a prime p, all trials fail with probability below 10^(-35)
def power_mod(base: int, exponent: int, p: int) -> int:
"""base^exponent mod p by binary exponentiation with a loop."""
result = 1 % p
base %= p
while exponent > 0:
if exponent & 1:
result = result * base % p
base = base * base % p
exponent >>= 1
return result
def cipolla_sqrt(n: int, p: int, rng: Optional[random.Random] = None) -> Optional[int]:
"""Return x in [0, p) with x^2 = n (mod p), or None if n is not a square."""
if p < 2:
raise ValueError("the modulus must be a prime")
if p == 2:
return n % 2 # every residue modulo 2 is its own square root
if rng is None:
rng = random.Random()
n %= p # @line L1
if n == 0:
return 0
euler = power_mod(n, (p - 1) // 2, p) # @line L2
if euler == p - 1:
return None # Euler's criterion: n is a nonresidue
if euler != 1:
raise ValueError("the modulus is not prime")
for _ in range(MAX_TRIALS):
a = rng.randrange(p) # @line L3
d = (a * a - n) % p
if power_mod(d, (p - 1) // 2, p) == p - 1: # @line L4
break
else:
raise ValueError("no nonresidue found; the modulus is probably not prime")
x, y = 1, 0 # @line L5
b, c = a, 1
e = (p + 1) // 2 # @line L6
while e > 0: # @line L7
if e & 1: # @line L8
x, y = (x * b + y * c % p * d) % p, (x * c + y * b) % p
b, c = (b * b + c * c % p * d) % p, 2 * b * c % p # @line L9
e >>= 1 # @line L10
if y != 0 or x * x % p != n:
raise ValueError("the modulus is not prime")
return x # @line L11
test_cipolla_algorithm.py 151 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/cipolla-algorithm/. The authors are listed in prov/W39G65.jsonl.
import random
import pytest
from cipolla_algorithm import cipolla_sqrt, power_mod
def is_probable_prime(n, rng):
"""Miller-Rabin with Python's pow, used only to build reference primes."""
if n < 2:
return False
for q in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37):
if n % q == 0:
return n == q
d, s = n - 1, 0
while d % 2 == 0:
d //= 2
s += 1
for _ in range(30):
a = rng.randrange(2, n - 1)
x = pow(a, d, n)
if x in (1, n - 1):
continue
for _ in range(s - 1):
x = x * x % n
if x == n - 1:
break
else:
return False
return True
def is_square_reference(n, p):
"""Euler's criterion with Python's pow, used only as the reference."""
n %= p
return n == 0 or pow(n, (p - 1) // 2, p) == 1
SMALL_PRIMES = [p for p in range(3, 400) if all(p % q for q in range(2, int(p**0.5) + 1))]
def large_primes():
rng = random.Random(10)
primes = [10**9 + 7, 998244353, 2**61 - 1, 2**127 - 1, 2**255 - 19, 2**521 - 1, 2**607 - 1]
for e in (40, 200, 500, 900):
while True: # a prime p = 1 + 2^e * odd, so p - 1 has a large power of two
p = (rng.randrange(1, 2**64) | 1) * 2**e + 1
if is_probable_prime(p, rng):
primes.append(p)
break
while True: # a random prime of several hundred decimal digits
p = rng.randrange(10**300, 10**301) | 1
if is_probable_prime(p, rng):
primes.append(p)
break
return primes
LARGE_PRIMES = large_primes()
def test_wikipedia_example():
class Fixed:
def randrange(self, p):
return 2
assert cipolla_sqrt(10, 13, Fixed()) == 6
assert cipolla_sqrt(10, 13, random.Random(0)) in (6, 7)
def test_power_mod():
rng = random.Random(1)
for _ in range(300):
m = rng.randrange(1, 10**rng.randrange(1, 200))
b, e = rng.randrange(-10**50, 10**50), rng.randrange(0, 10**30)
assert power_mod(b, e, m) == pow(b, e, m)
def test_every_residue_small_primes():
rng = random.Random(2)
for p in SMALL_PRIMES:
squares = {x * x % p for x in range(p)}
for n in range(p):
x = cipolla_sqrt(n, p, rng)
if n in squares:
assert x is not None and 0 <= x < p and x * x % p == n
else:
assert x is None
def test_edge_values():
rng = random.Random(3)
for p in SMALL_PRIMES + LARGE_PRIMES:
assert cipolla_sqrt(0, p, rng) == 0
assert cipolla_sqrt(p, p, rng) == 0
assert cipolla_sqrt(7 * p, p, rng) == 0
assert cipolla_sqrt(1, p, rng) in (1, p - 1)
x = cipolla_sqrt(p - 1, p, rng)
if p % 4 == 1:
assert x * x % p == p - 1
else:
assert x is None
def test_random_large_primes():
rng = random.Random(4)
count = 0
for p in LARGE_PRIMES:
for _ in range(100 if p < 2**130 else 12):
r = rng.randrange(p)
n = r * r % p
x = cipolla_sqrt(n, p, rng)
assert x * x % p == n and x in (r, (p - r) % p)
m = rng.randrange(p)
x = cipolla_sqrt(m, p, rng)
assert (x is not None) == is_square_reference(m, p)
if x is not None:
assert x * x % p == m % p
count += 2
for _ in range(500):
p = rng.choice(SMALL_PRIMES[5:])
n = rng.randrange(-10**300, 10**300)
x = cipolla_sqrt(n, p, rng)
assert (x is not None) == is_square_reference(n, p)
count += 1
assert count >= 1000
def test_modulus_two():
assert cipolla_sqrt(0, 2) == 0
assert cipolla_sqrt(1, 2) == 1
assert cipolla_sqrt(10**300 + 1, 2) == 1
def test_bad_and_composite_moduli():
for p in (-5, 0, 1):
with pytest.raises(ValueError):
cipolla_sqrt(1, p)
rng = random.Random(5)
composites = [4, 8, 16, 2**64, 9, 15, 21, 341, 561, 1105, 3**5, 10**300 + 1, 2**1000]
for m in composites:
for _ in range(30):
n = rng.randrange(m)
try:
x = cipolla_sqrt(n, m, rng)
except ValueError:
continue
assert x is None or x * x % m == n % m
What the next level needs
- the tests did not run in this build