Solving rational linear systems

The solvers of this module find exact solutions of $AX = B$ over $\mathbb{Q}$ for rational or integer matrices $A$ and $B$, returning the solution as an fmpq_mat_t in canonical form. All of them reduce the problem to integer matrices and call the solvers of fmpz_mat.h, which work modulo primes and reconstruct the rational solution; see Solving linear systems over the rationals.

Reduction to integer systems#

Multiplying row $i$ of the augmented matrix $(A \mid B)$ by the least common denominator $d_i$ of its entries gives an integer system $A'X = B'$ with the same solutions, since $A' = \operatorname{diag}(d) A$ and $B' = \operatorname{diag}(d) B$ with $\operatorname{diag}(d)$ invertible. fmpq_mat_get_fmpz_mat_rowwise_2 computes $A'$ and $B'$. The integer solver returns $X$ as an integer matrix and a denominator, which fmpq_mat_set_fmpz_mat_div_fmpz turns into a rational matrix in canonical form. The variants whose names contain fmpz_mat take $A$ and $B$ as integer matrices directly.

Clearing denominators row by row keeps the integer entries as small as possible: for a matrix with a common structure of denominators, such as the Hilbert matrix, the integers $d_i$ are the least common multiples of the denominators in each row, not of the whole matrix.

Choice of algorithm#

For a square $n \times n$ matrix $A$ and $k$ right-hand sides, fmpq_mat_solve and fmpq_mat_solve_fmpz_mat use

• for $n \le 15$: fmpz_mat_solve, that is Cramer's rule for $n \le 3$ and fraction-free elimination otherwise (the fraction_free variants);
• for $n \ge 16$ and $2k \ge n$: multimodular solving (fmpz_mat_solve_multi_mod_den, the multi_mod variants);
• otherwise: $p$-adic lifting (fmpz_mat_solve_dixon_den, the dixon variants).

All of them return 1 if $A$ is nonsingular and 0 otherwise, and in the latter case leave $X$ unchanged. The modular solvers are output-sensitive: their running time depends on the size of the solution rather than on a worst-case bound.

fmpq_mat_can_solve accepts any $m \times n$ matrix $A$. It returns 1 with a solution if the system is consistent, and 0 if it is not; the answer is proven, by exact verification of the solution or by an exact certificate of inconsistency. For $m \le 15$ it uses fraction-free elimination (fmpq_mat_can_solve_fraction_free), otherwise the modular method of fmpz_mat_can_solve; fmpq_mat_can_solve_dixon and fmpq_mat_can_solve_multi_mod force the solver used for the nonsingular pivot block.

fmpq_mat_inv solves $A' X = \operatorname{diag}(d)$ with fmpq_mat_solve_fmpz_mat, where $A' = \operatorname{diag}(d) A$, so that $X = A^{-1}$; it uses explicit formulas for dimensions up to 2.

Example#

The $10 \times 10$ Hilbert matrix $H$, with entries $1/(i + j + 1)$, is notoriously ill-conditioned but has an integral inverse. Solving $Hx = (1, \ldots, 1)^T$ exactly, inverting $H$, and testing two singular systems:

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

int main(void)
{
    fmpq_mat_t H, B, X, Hinv, A, Y;
    slong i, n = 10;

    fmpq_mat_init(H, n, n);
    fmpq_mat_init(B, n, 1);
    fmpq_mat_init(X, n, 1);
    fmpq_mat_init(Hinv, n, n);

    fmpq_mat_hilbert_matrix(H);
    for (i = 0; i < n; i++)
        fmpq_one(fmpq_mat_entry(B, i, 0));

    flint_printf("solvable: %d\n", fmpq_mat_solve(X, H, B));
    flint_printf("x_0 = %{fmpq}, x_9 = %{fmpq}\n",
                 fmpq_mat_entry(X, 0, 0), fmpq_mat_entry(X, 9, 0));

    flint_printf("invertible: %d\n", fmpq_mat_inv(Hinv, H));
    flint_printf("inverse is integral: %d\n", fmpq_mat_is_integral(Hinv));
    flint_printf("entry (4, 4) of the inverse: %{fmpq}\n", fmpq_mat_entry(Hinv, 4, 4));

    /* a singular system: rows 0, 1 and their sum */
    fmpq_mat_init(A, 3, 2);
    fmpq_mat_init(Y, 2, 1);
    fmpq_mat_clear(B);
    fmpq_mat_init(B, 3, 1);
    fmpq_set_si(fmpq_mat_entry(A, 0, 0), 1, 2);
    fmpq_set_si(fmpq_mat_entry(A, 0, 1), 1, 3);
    fmpq_set_si(fmpq_mat_entry(A, 1, 0), 1, 4);
    fmpq_set_si(fmpq_mat_entry(A, 1, 1), -1, 5);
    fmpq_add(fmpq_mat_entry(A, 2, 0), fmpq_mat_entry(A, 0, 0), fmpq_mat_entry(A, 1, 0));
    fmpq_add(fmpq_mat_entry(A, 2, 1), fmpq_mat_entry(A, 0, 1), fmpq_mat_entry(A, 1, 1));
    fmpq_set_si(fmpq_mat_entry(B, 0, 0), 1, 1);
    fmpq_set_si(fmpq_mat_entry(B, 1, 0), 2, 1);
    fmpq_set_si(fmpq_mat_entry(B, 2, 0), 3, 1);

    flint_printf("consistent: %d, ", fmpq_mat_can_solve(Y, A, B));
    flint_printf("y = (%{fmpq}, %{fmpq})\n",
                 fmpq_mat_entry(Y, 0, 0), fmpq_mat_entry(Y, 1, 0));

    fmpq_set_si(fmpq_mat_entry(B, 2, 0), 4, 1);
    flint_printf("consistent: %d\n", fmpq_mat_can_solve(Y, A, B));

    fmpq_mat_clear(H);
    fmpq_mat_clear(B);
    fmpq_mat_clear(X);
    fmpq_mat_clear(Hinv);
    fmpq_mat_clear(A);
    fmpq_mat_clear(Y);
    return 0;
}

solvable: 1
x_0 = -10, x_9 = 923780
invertible: 1
inverse is integral: 1
entry (4, 4) of the inverse: 176752976400
consistent: 1, y = (52 / 11, -45 / 11)
consistent: 0