4.5. The multimodular loop🔗

Hex.Modular.crtLoop runs the multimodular method. It takes a function computing the image of the answer modulo a given modulus, a test deciding whether the accumulated residues determine the answer, a supply of moduli, and fuel:

def

Hex.Modular.crtLoop.{u_1} {k : ℕ} {α : Type u_1}
  (image : ℕ → Option (Vector ℤ k))
  (accept : Hex.Modular.CrtVec k → Option α) (supply : Array ℕ)
  (fuel : ℕ) : Option α
Hex.Modular.crtLoop.{u_1} {k : ℕ}
  {α : Type u_1}
  (image : ℕ → Option (Vector ℤ k))
  (accept :
    Hex.Modular.CrtVec k → Option α)
  (supply : Array ℕ) (fuel : ℕ) : Option α

Fold moduli into a vector CRT state until the caller-supplied exact check accepts a reconstruction. A rejected image never enters the state, and the fuel bounds the number of supply entries inspected.

The image function returns none for a modulus that does not work. For a determinant computed by elimination, that is a prime modulo which no invertible pivot can be found. The loop skips such a modulus, and also one that push rejects because it is not coprime to those before it. Every modulus inspected, used or skipped, consumes one unit of fuel. The test is called after each successful push, and returns none to ask for another modulus or some of the answer, which the loop returns.

For H_n a bound gives a test. The least common multiple L of 1, …, n is a common denominator of the terms, so the denominator of H_n divides L, and its numerator is at most n L since H_n ≤ n. So once the modulus exceeds 2 · n L · L, ratRecon? with these bounds must find H_n. The image function rejects primes p ≤ n, at which some term 1 / k has no residue:

def lcmUpTo (n : Nat) : Nat :=
  (List.range n).foldl (fun l k => Nat.lcm l (k + 1)) 1

def harmonicLoop (n : Nat) : Option Rat :=
  let Q : Int := lcmUpTo n
  let P := n * Q
  crtLoop
    (fun p => if p ≤ n then none
      else some #v[harmonicMod n p])
    (fun s => if 2 * P * Q < s.modulus then
      ratRecon? s.value[0] s.modulus P Q else none)
    primes.toArray primes.length

some (13943237577224054960759 / 3099044504245996706400)#eval harmonicLoop 50

some (13943237577224054960759 / 3099044504245996706400)

The test first passes after the fifth prime. The theorem about the loop identifies the initial segment of the supply that it used and the state on which the test returned the result; a Hex.Modular.CrtTrace records how that state was reached from the images.

theorem

Hex.Modular.crtLoop_trace.{u_1} {k : ℕ} {α : Type u_1}
  {image : ℕ → Option (Vector ℤ k)}
  {accept : Hex.Modular.CrtVec k → Option α} {supply : Array ℕ}
  {fuel : ℕ} {x : α}
  (h : Hex.Modular.crtLoop image accept supply fuel = some x) :
  ∃ consumed state,
    0 < consumed ∧
      consumed ≤ fuel ∧
        consumed ≤ supply.size ∧
          Hex.Modular.CrtTrace image supply consumed state ∧
            accept state = some x
Hex.Modular.crtLoop_trace.{u_1} {k : ℕ}
  {α : Type u_1}
  {image : ℕ → Option (Vector ℤ k)}
  {accept :
    Hex.Modular.CrtVec k → Option α}
  {supply : Array ℕ} {fuel : ℕ} {x : α}
  (h :
    Hex.Modular.crtLoop image accept
        supply fuel =
      some x) :
  ∃ consumed state,
    0 < consumed ∧
      consumed ≤ fuel ∧
        consumed ≤ supply.size ∧
          Hex.Modular.CrtTrace image
              supply consumed state ∧
            accept state = some x

A successful loop result is accepted on a state reached by replaying the exact consumed supply prefix. The consumed prefix is nonempty because the loop tests accept only after a successful push, and it cannot exceed fuel.

theorem

Hex.Modular.CrtTrace.congr {k : ℕ} {image : ℕ → Option (Vector ℤ k)}
  {supply : Array ℕ} {consumed : ℕ} {state : Hex.Modular.CrtVec k}
  (trace : Hex.Modular.CrtTrace image supply consumed state)
  (x : Vector ℤ k)
  (himage :
    ∀ (j : ℕ) (hj : j < supply.size),
      j < consumed →
        ∀ (r : Vector ℤ k),
          image supply[j] = some r →
            ∀ (i : Fin k), r[i] % ↑supply[j] = x[i] % ↑supply[j])
  (i : Fin k) : state.value[i] % ↑state.modulus = x[i] % ↑state.modulus
Hex.Modular.CrtTrace.congr {k : ℕ}
  {image : ℕ → Option (Vector ℤ k)}
  {supply : Array ℕ} {consumed : ℕ}
  {state : Hex.Modular.CrtVec k}
  (trace :
    Hex.Modular.CrtTrace image supply
      consumed state)
  (x : Vector ℤ k)
  (himage :
    ∀ (j : ℕ) (hj : j < supply.size),
      j < consumed →
        ∀ (r : Vector ℤ k),
          image supply[j] = some r →
            ∀ (i : Fin k),
              r[i] % ↑supply[j] =
                x[i] % ↑supply[j])
  (i : Fin k) :
  state.value[i] % ↑state.modulus =
    x[i] % ↑state.modulus

A traced state agrees modulo its accumulated modulus with every vector whose coordinates agree with all successful images in the consumed prefix.

The correctness of a particular use combines these with the theorems about the test. For harmonicLoop, write H_n = N / D and suppose every image is the residue of H_n. Let M be the product of the primes used, which are larger than n and so coprime to D, and let u be an inverse of D modulo M. Then CrtTrace.congr, applied to the vector #v[N * u], shows that the final value v is congruent to N u modulo M, so D v ≡ N. This is the congruence hypothesis of ratRecon?_complete, which then shows that the result is H_n.