Rational reconstruction

Multimodular and $p$-adic algorithms compute in $\mathbb{Z}/m\mathbb{Z}$, where a rational number $n/d$ with $\gcd(d, m) = 1$ is represented by the residue $a \equiv n d^{-1} \pmod m$. When the answer of a computation over $\mathbb{Q}$ — the solution of a linear system, the coefficients of a polynomial — has been found modulo a large enough $m$, rational reconstruction recovers it from $a$. The method goes back to Wang [WGD1982]; see also [CollinsEncarnacion1995], [GG2003, §5.10] and [Sho2009, §4.6].

Throughout, $m > 2$, $0 \le a < m$, and $N, D \ge 1$ are bounds with $2ND < m$.

Uniqueness#

The condition $2ND < m$ cannot be weakened much: the number of fractions with $|n| \le N$ and $0 < d \le D$ in lowest terms is about $(12/\pi^2) ND$, so their residues cannot all be distinct if $ND$ is much larger than $m$. With the balanced choice $N = D \approx \sqrt{m/2}$, about $6/\pi^2 \approx 61\%$ of all residues are reconstructible.

The extended Euclidean algorithm#

Apply the extended Euclidean algorithm to $r_0 = m$ and $r_1 = a$: with the quotients $q_i = \lfloor r_{i-1} / r_i \rfloor$,

starting from $(s_0, t_0) = (1, 0)$ and $(s_1, t_1) = (0, 1)$, until $r_{\ell + 1} = 0$. Then $r_i = s_i m + t_i a$, so $r_i \equiv t_i a \pmod m$ for every $i$. We use the following standard facts [GG2003, §3.2, §4.6].

1. The $t_i$ alternate in sign and $|t_1| < |t_2| < \cdots < |t_{\ell+1}|$ (with $|t_2| = q_1 \ge 1 = |t_1|$ the only possible equality).
2. $s_i t_{i+1} - s_{i+1} t_i = (-1)^i$, so $\gcd(s_i, t_i) = 1$.
3. $r_i |t_{i+1}| + r_{i+1} |t_i| = m$ for $0 \le i \le \ell$.
4. The convergents of the continued fraction of $a/m$ are the fractions $-s_i/t_i$, $1 \le i \le \ell + 1$.

Fact 3 follows from $r_i t_{i+1} - r_{i+1} t_i = (s_i t_{i+1} - s_{i+1} t_i) m = \pm m$ and the alternating signs of the $t_i$.

The remainders $r_i$ shrink and the cofactors $t_i$ grow as the algorithm proceeds, with $r_{i-1} |t_i| \le m$; the bound $N$ decides where the algorithm stops, trading the size of the numerator against the size of the denominator. Monagan's maximal quotient variant [Mon2004] does not fix $N$ in advance but stops after the largest quotient $q_i$, which succeeds with a modulus only slightly larger than $2|n|d$ with high probability; FLINT implements the bounded version, which is a proof when the bounds are known.

Complexity#

The Euclidean algorithm on $b$-bit numbers costs $O(b^2)$ bit operations. Lehmer's method [Leh1938] simulates several steps with word-size leading parts of $r_{i-1}$ and $r_i$, so that most of the work is done with single-word $2 \times 2$ matrices. The half-GCD algorithm [Sch1971], [ThullYap1990] computes the quotient sequence up to a given size recursively from the leading halves and costs $O(M(b) \log b)$. The stopping rule $r_j \le N$ is a partial GCD, which the half-GCD algorithm handles naturally: it reduces $(m, a)$ to a pair straddling a given size, then continues with Lehmer steps.

Several entries with a common denominator#

When a vector or matrix of rationals is reconstructed, the entries often share most of their denominator. Multiplying each residue by the denominator $\delta$ of the entries found so far, the next entry is usually an integer $\delta x \bmod m$ in the symmetric range $|\cdot| \le N$, recognised without any GCD; only otherwise is the entry reconstructed, with the denominator bound $D/\delta$, and $\delta$ multiplied by the new denominator.

Use in FLINT#

_fmpq_reconstruct_fmpz_2 implements Algorithm 1. After the quick integer tests of step 1, it chooses by the size of $m$:

• for $m$ of one or two limbs, a Euclidean algorithm in one or two words, with a binary GCD for the final coprimality test;
• up to 12 limbs, Lehmer's method on stack arrays of limbs;
• for larger $m$, Lehmer's method on fmpz, switching to the half-GCD algorithm _fmpq_hgcd of the continued fraction code (Continued fractions of rational numbers) when $m$ has at least FMPQ_RECONSTRUCT_HGCD_CUTOFF (500) more limbs than $N$, and to plain Euclidean steps when they differ by fewer than 3 limbs.

_fmpq_reconstruct_fmpz_2_naive is the plain algorithm on fmpz, used to check the fast version in tests. fmpq_reconstruct_fmpz uses the balanced bounds $N = D = \lfloor \sqrt{\lfloor (m-1)/2 \rfloor} \rfloor$.

The reconstruction returns rational results in fmpq_mat_set_fmpz_mat_mod_fmpz, fmpz_mat_det_divisor and the Dixon solver fmpz_mat_solve_dixon_den (through _fmpz_mat_reconstruct_matwise), which all use the common denominator technique above, and in fmpq_poly_interpolate_multi_mod and fmpz_poly_roots_fmpq. How the multimodular and $p$-adic linear algebra algorithms choose the bounds and decide when to stop is described in multimodular and $p$-adic linear algebra.