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.