Rational reconstruction

A multimodular or $p$-adic algorithm whose answer is a rational number computes it as a residue $a$ modulo a large integer $m$, a product of primes or a prime power. Rational reconstruction recovers the fraction $n/d$ from $a = n d^{-1} \bmod m$, provided that $|n|$ and $d$ are small compared with $m$. fmpq_mod_fmpz computes the residue of a fraction, and the reconstruction functions invert it:

• fmpq_reconstruct_fmpz_2 takes explicit bounds $N$ and $D$ with $2ND < m$ and finds the unique $n/d$ with $|n| \le N$ and $0 < d \le D$, if it exists;
• fmpq_reconstruct_fmpz uses the balanced bounds $N = D = \lfloor\sqrt{\lfloor (m-1)/2 \rfloor}\rfloor$.

The residue must be reduced, $0 \le a < m$, and $m > 2$. The result is canonical: $\gcd(n, d) = 1$ and $d > 0$.

Success and failure#

The bound $2ND < m$ makes the answer unique: two fractions $n_1/d_1 \ne n_2/d_2$ within the bounds have different residues, because $m$ would otherwise divide $n_1 d_2 - n_2 d_1$, a nonzero number of absolute value at most $2ND$. A return value of 0 therefore proves that no fraction within the bounds has residue $a$.

A return value of 1 proves only that the fraction found is the unique one within the bounds. If the true answer of a computation is larger than the bounds, the function either fails or returns a different, small fraction with the same residue. About $6/\pi^2 \approx 61\%$ of all residues have a reconstruction with the balanced bounds, so a successful reconstruction of a residue that is not known to come from a small fraction is weak evidence. Multimodular algorithms therefore either use a priori bounds on the answer (Hadamard's bound for determinants and solutions of linear systems, for example) and choose $m$ larger than $2ND$, or verify the reconstructed result, or add primes until the reconstruction stabilises.

When the shape of the answer is known, unbalanced bounds help. For a fraction whose denominator divides a known $d_0$, take $D = d_0$ and $N = \lfloor (m-1)/(2 d_0) \rfloor$: almost all of $m$ is then available for the numerator.

Algorithm#

The fraction is found among the remainders and cofactors of the extended Euclidean algorithm applied to $m$ and $a$, stopping at the first remainder $r_j \le N$; then $n/d = \pm r_j / t_j$ where $t_j$ is the cofactor of $a$, and the reconstruction succeeds if $|t_j| \le D$ and $\gcd(r_j, t_j) = 1$ (Wang's algorithm [WGD1982]). The implementation works in one or two words or on fixed arrays of up to 12 limbs for small $m$, with Lehmer's method for larger $m$, and with the half-GCD algorithm when the modulus has at least FMPQ_RECONSTRUCT_HGCD_CUTOFF (500) more limbs than $N$, for a total cost of $O(M(\log m) \log\log m)$. Since the search stops at a remainder of the size of $N$, its cost depends on the size of $m$ relative to $N$, not only on the size of $m$. _fmpq_reconstruct_fmpz_2_naive is the plain Euclidean algorithm on fmpz, used for testing.

Example#

The modulus is the prime $m = 10^9 + 7$, for which the balanced bound is $N = D = 22360$. The fraction $-355/113$ is recovered. The fraction $123456/789 = 41152/263$ has a numerator above the balanced bound; the balanced reconstruction fails, and unbalanced bounds $N = 400000$, $D = 1000$ (with $2ND < m$) recover it. The residue of $1234567/89$ is reconstructed with the balanced bounds, but to a different fraction of small height, which shows that a successful reconstruction has to be checked when the answer is not known to be within the bounds.

#include <flint/flint.h>
#include <flint/fmpz.h>
#include <flint/fmpq.h>

static void
show(const fmpq_t x, const fmpz_t m, const fmpz_t N, const fmpz_t D)
{
    fmpz_t a;
    fmpq_t y;
    int ok;

    fmpz_init(a);
    fmpq_init(y);
    fmpq_mod_fmpz(a, x, m);
    if (N == NULL)
        ok = fmpq_reconstruct_fmpz(y, a, m);
    else
        ok = fmpq_reconstruct_fmpz_2(y, a, m, N, D);

    fmpq_print(x);
    flint_printf(" -> ");
    fmpz_print(a);
    if (ok)
    {
        flint_printf(" -> ");
        fmpq_print(y);
        flint_printf("\n");
    }
    else
        flint_printf(" -> no fraction within the bounds\n");

    fmpz_clear(a);
    fmpq_clear(y);
}

int main(void)
{
    fmpz_t m, N, D;
    fmpq_t x;

    fmpz_init_set_ui(m, 1000000007);
    fmpz_init_set_ui(N, 400000);
    fmpz_init_set_ui(D, 1000);
    fmpq_init(x);

    fmpq_set_si(x, -355, 113);
    show(x, m, NULL, NULL);

    fmpq_set_si(x, 123456, 789);
    show(x, m, NULL, NULL);
    show(x, m, N, D);

    fmpq_set_si(x, 1234567, 89);
    show(x, m, NULL, NULL);

    fmpz_clear(m);
    fmpz_clear(N);
    fmpz_clear(D);
    fmpq_clear(x);
    return 0;
}

-355/113 -> 699115046 -> -355/113
41152/263 -> 726235903 -> no fraction within the bounds
41152/263 -> 726235903 -> 41152/263
1234567/89 -> 955070058 -> -14062/18629