Reconstruct a mod m as a rational with numerator absolute value at most
P and positive denominator at most Q. The truncated Euclidean candidate is
normalized and all output conditions are checked before it is returned.
4.4. Rational reconstruction
A fraction p / q with q coprime to m has a residue modulo m, the
residue of p · q⁻¹. For example 2 / 3 is 68 modulo 101, since
3 · 68 = 204 = 2 · 101 + 2. Rational reconstruction goes back: given a
and m and bounds P and Q, it looks for p / q with
q a ≡ p (mod m), |p| ≤ P, 0 < q ≤ Q.
The algorithm is Wang's. Reduce a modulo m, and run the extended Euclidean
algorithm on m and a.
Its rows are remainders rⱼ = sⱼ m + tⱼ a, so each satisfies
tⱼ a ≡ rⱼ (mod m), and as the remainders decrease the coefficients |tⱼ|
increase. Stop at the first row with rⱼ ≤ P and take the candidate
rⱼ / tⱼ, which Rat.divInt puts in lowest terms with a positive
denominator. Hex.Modular.euclidUntil returns that row. For 68 modulo
101 the remainders are 101, 68, 33, 2, and the row for 2 has t = 3:
#eval euclidUntil 101 68 7
def
def
Symmetric rational reconstruction with
P = Q = ⌊√((m-1)/2)⌋, which guarantees 2 P Q < m.
#eval (ratRecon? 68 101 7 7, ratReconWide? 68 101,
ratReconWide? 33 101)
Since the result is checked before it is returned, it is always a solution:
Hex.Modular.ratRecon?_congr gives the congruence,
Hex.Modular.ratRecon?_bounds the bounds, and
Hex.Modular.ratRecon?_den_coprime shows that the denominator is
coprime to the modulus, so that the fraction really has a residue.
4.4.1. Uniqueness and completeness
Without a condition relating the bounds to the modulus, a residue can be
the image of several fractions within them. Modulo 101, both 2 / 3 and
-23 / 16 reduce to 68, and both have numerator and denominator at most
25:
#guard ratReconCheck 68 101 25 25 (2 / 3)
#guard ratReconCheck 68 101 25 25 (-23 / 16)
The condition is 2 P Q < m. If p₁ / q₁ and p₂ / q₂ both satisfy the
congruence, then m divides p₁ q₂ - p₂ q₁, whose absolute value is at most
2 P Q. So if 2 P Q < m it is zero and the fractions are equal.
Hex.Modular.ratReconWide? chooses P = Q = ⌊√((m - 1) / 2)⌋, the
largest equal bounds satisfying the condition.
theorem
Hex.Modular.ratRecon_unique {a P Q : ℤ} {m : ℕ} {y₁ y₂ : Rat} (hm : 2 * P * Q < ↑m) (h₁ : (Int.ofNat y₁.den * a - y₁.num) % ↑m = 0) (h₂ : (Int.ofNat y₂.den * a - y₂.num) % ↑m = 0) (b₁ : ↑y₁.num.natAbs ≤ P ∧ ↑y₁.den ≤ Q) (b₂ : ↑y₂.num.natAbs ≤ P ∧ ↑y₂.den ≤ Q) : y₁ = y₂Hex.Modular.ratRecon_unique {a P Q : ℤ} {m : ℕ} {y₁ y₂ : Rat} (hm : 2 * P * Q < ↑m) (h₁ : (Int.ofNat y₁.den * a - y₁.num) % ↑m = 0) (h₂ : (Int.ofNat y₂.den * a - y₂.num) % ↑m = 0) (b₁ : ↑y₁.num.natAbs ≤ P ∧ ↑y₁.den ≤ Q) (b₂ : ↑y₂.num.natAbs ≤ P ∧ ↑y₂.den ≤ Q) : y₁ = y₂
Rational reconstruction is unique whenever twice the product of the two bounds is strictly smaller than the modulus.
Under the same condition the Euclidean algorithm finds the fraction whenever
it exists. The key fact is that |tⱼ| rⱼ₋₁ ≤ m for every row, which places
the candidate within the bounds whenever any solution is; see von zur Gathen
and Gerhard, Modern Computer Algebra, and Wang, Guy and Davenport,
"P-adic reconstruction of rational numbers" (1982).
theorem
Hex.Modular.ratRecon?_complete {a P Q : ℤ} {m : ℕ} {y : Rat} (hm : 2 * P * Q < ↑m) (hy : (Int.ofNat y.den * a - y.num) % ↑m = 0) (hb : ↑y.num.natAbs ≤ P ∧ ↑y.den ≤ Q) : Hex.Modular.ratRecon? a m P Q = some yHex.Modular.ratRecon?_complete {a P Q : ℤ} {m : ℕ} {y : Rat} (hm : 2 * P * Q < ↑m) (hy : (Int.ofNat y.den * a - y.num) % ↑m = 0) (hb : ↑y.num.natAbs ≤ P ∧ ↑y.den ≤ Q) : Hex.Modular.ratRecon? a m P Q = some y
Under the uniqueness bound, bounded reconstruction finds every rational that satisfies the congruence and bounds.
So none from ratRecon? with 2 P Q < m means that no fraction within the
bounds has the residue a. In a multimodular computation it usually means
that the modulus is not yet large enough.
4.4.2. Common denominators
By Cramer's rule, the solution of a linear system A x = b, with A a
nonsingular square integer matrix and b an integer vector, is a vector of
fractions whose common denominator divides the determinant of A.
Hex.Modular.ratReconVec? reconstructs such a vector, returning the
numerators and the common denominator d. It reconstructs the first entry
on its own. For each later entry aᵢ it first tries symMod (d aᵢ) m as the
numerator, which is one multiplication, and runs a Euclidean reconstruction
only when that numerator is too large. Then it replaces d by the least
common multiple of d and the new denominator, and scales the numerators
found so far to match. At the end it divides the numerators and d by their
common greatest common divisor. Modulo 101, 51 is 1 / 2
and 76 is 1 / 4:
def halfQuarter : Vector Int 2 := #v[51, 76]
#eval ratReconVec? halfQuarter 101 2 4
The bound Q now bounds the common denominator, and P the numerators over
that common denominator, not those of each fraction in lowest terms.
Completeness needs the numerators and the denominator together to have
greatest common divisor 1, although a single numerator may share a factor
with d:
theorem
Hex.Modular.ratReconVec?_complete {k : ℕ} {a y : Vector ℤ k} {m : ℕ} {P Q d : ℤ} (hm : 2 * P * Q < ↑m) (hP : 0 ≤ P) (hd : 0 < d) (hdQ : d ≤ Q) (hy : ∀ (i : Fin k), (d * a[i] - y[i]) % ↑m = 0 ∧ ↑y[i].natAbs ≤ P) (hred : ∀ (g : ℤ), (∀ (i : Fin k), g ∣ y[i]) → g ∣ d → g ∣ 1) : Hex.Modular.ratReconVec? a m P Q = some (y, d)Hex.Modular.ratReconVec?_complete {k : ℕ} {a y : Vector ℤ k} {m : ℕ} {P Q d : ℤ} (hm : 2 * P * Q < ↑m) (hP : 0 ≤ P) (hd : 0 < d) (hdQ : d ≤ Q) (hy : ∀ (i : Fin k), (d * a[i] - y[i]) % ↑m = 0 ∧ ↑y[i].natAbs ≤ P) (hred : ∀ (g : ℤ), (∀ (i : Fin k), g ∣ y[i]) → g ∣ d → g ∣ 1) : Hex.Modular.ratReconVec? a m P Q = some (y, d)
Under the uniqueness bound, common-denominator reconstruction finds every
pair that is reduced as a whole, positive and bounded in its denominator,
and congruent and bounded in every coordinate. Reducedness is not
decoration: at m = 18, a = (10), the pair (2, 2) satisfies every other
hypothesis, but its reduced rational 1 does not satisfy the congruence and
the reconstruction returns none.
4.4.3. Without bounds
Some computations have no useful bound on the answer, or only one far larger
than the answer usually is. Monagan's maximal quotient rational
reconstruction needs none. If p / q has residue a and |p| q is much
smaller than m, then the quotient in the step of the Euclidean algorithm
after the row for p / q is about m / (|p| q), usually much larger than the
other quotients. Hex.Modular.ratReconMaxQuot? takes its candidate
from the row just before the step with the largest quotient. Here m = 2^40 - 87 and a is the residue of
123456789 / 2. The bound √(m / 2) is 741455, so ratReconWide? fails,
while the maximal quotient finds the fraction:
#eval (ratReconWide? 549817542239 1099511627689,
ratReconMaxQuot? 549817542239 1099511627689)
This is a heuristic, and its theorem promises only the congruence. A computation using it needs its own check of the answer.
theorem
Hex.Modular.ratReconMaxQuot?_congr {a : ℤ} {m : ℕ} {x : Rat} (h : Hex.Modular.ratReconMaxQuot? a m = some x) : (Int.ofNat x.den * a - x.num) % ↑m = 0Hex.Modular.ratReconMaxQuot?_congr {a : ℤ} {m : ℕ} {x : Rat} (h : Hex.Modular.ratReconMaxQuot? a m = some x) : (Int.ofNat x.den * a - x.num) % ↑m = 0
A maximal-quotient candidate always satisfies the modular congruence; no claim is made that it is the intended rational.