4.1. Introduction
Many exact computations with integers and rationals are done fastest by not
doing them over the integers at all. The multimodular method computes the
answer modulo several primes, where every number fits in a machine word,
combines the images by the Chinese remainder theorem into the answer modulo
the product of the primes, and reads off the integer or rational answer from
that residue. HexModular provides the last two steps: Chinese remaindering,
incrementally or for many moduli at once, and rational reconstruction, which
recovers a fraction p / q from the residue of p · q⁻¹. It also provides the
loop that adds primes until a caller's test accepts the result.
The harmonic number H₅₀ = 1 + 1/2 + ⋯ + 1/50 is a fraction whose numerator
has 74 bits and whose denominator has 72. Modulo a prime p > 50 it can be
computed with arithmetic on numbers below p: keep H_k as a fraction
a / b of residues, and at the end multiply by b^(p-2), which is the
inverse of b by Fermat's little theorem:
open Hex.Modular
namespace HexModularChapter
def harmonicMod (n p : Nat) : Int :=
let (a, b) := (List.range n).foldl
(fun (a, b) k =>
((a * (k + 1) + b) % p, b * (k + 1) % p))
(0, 1)
a * Nat.powMod b (p - 2) p % p
Here are the five largest primes below 2^31. The function recover combines
the images of H_n modulo the first j of them with
Hex.Modular.Crt.push, and then asks Hex.Modular.ratReconWide?
for a fraction whose numerator, in absolute value, and denominator are both at
most ⌊√((M - 1) / 2)⌋, where M is the product of the primes used:
def primes : List Nat :=
[2147483647, 2147483629, 2147483587,
2147483579, 2147483563]
def recover (n j : Nat) : Option Rat := do
let c ← (primes.take j).foldlM
(fun c p => c.push (harmonicMod n p) p) Crt.init
ratReconWide? c.value c.modulus
#eval recover 50 5
This is H₅₀, as an exact computation with Lean's Rat confirms:
def harmonic (n : Nat) : Rat :=
(List.range n).foldl (fun s k => s + 1 / (k + 1 : Rat)) 0
#guard recover 50 5 = some (harmonic 50)
With these primes, five suffice and four do not. The product of five is about
2^155, so the bound is about 2^77, which exceeds both parts of H₅₀. With three
primes the bound is about 2^46. H₅₀ is not within it, and there is a
different fraction within it with the same residue, which reconstruction
duly returns. With four primes there is none:
#eval (recover 50 3, recover 50 4)
So a multimodular computation needs a reason to stop: a bound on the size of the answer, or a check that the answer is right. The section The multimodular loop returns to this example with a bound.
Import HexModular to use the library; the Mathlib examples near the end of
this chapter also import Mathlib.Data.Int.ModEq and Mathlib.Data.ZMod.Basic.
HexModular does not depend on Mathlib, and has no theory companion. Its
theorems state congruences with % on Int, which is how Mathlib defines
Int.ModEq, so they can be used in Mathlib proofs directly; see
Using the theorems with Mathlib. Its only
dependency is HexArith, whose extended Euclidean algorithm
uses GMP in compiled code.