# gr_mat.h – dense matrices over generic rings

A `gr_mat_t` represents a matrix implemented as a dense array of entries in a generic ring R.

- In this module, the context object `ctx` always represents the coefficient ring R unless otherwise stated. Creating a context object representing a matrix space only becomes necessary when one wants to manipulate matrices using generic ring methods like `gr_add` instead of the designated matrix methods like `gr_mat_add`.

- Matrix functions generally assume that input as well as output operands have compatible shapes. Some functions return `GR_DOMAIN` for matrices with the wrong shape, but this is not always consistent.

- Some operations (like rank, LU factorization) generally only make sense when the base ring is an integral domain. Typically the algorithms designed for integral domains also work over non-integral domains as long as all inversions of nonzero elements succeed. If an inversion fails, the algorithm will return the `GR_DOMAIN` or `GR_UNABLE` flag. This might not yet be entirely consistent.

## Type compatibility

The `gr_mat` type has the same data layout as most FLINT, Arb and Calcium matrix types. Methods in this module can therefore be mixed freely with methods in the corresponding FLINT, Arb and Calcium modules when the underlying coefficient type is the same.

It is not directly compatible with the `nmod_mat` type, which stores modulus data as part of the matrix object.

## Types, macros and constants

type gr_mat_struct

type gr_mat_t

  Contains a pointer to an array of coefficients (`entries`), the number of rows (`r`), the number of columns (`c`), and an array to pointers marking the start of each row (`rows`).

  A `gr_mat_t` is defined as an array of length one of type `gr_mat_struct`, permitting a `gr_mat_t` to be passed by reference.

GR_MAT_ENTRY(mat, i, j, sz)

  Macro to access the entry at row i and column j of the matrix mat whose entries have size sz bytes.

gr_ptr gr_mat_entry_ptr(gr_mat_t mat, slong i, slong j, gr_ctx_t ctx)

  Function returning a pointer to the entry at row i and column j of the matrix mat. The indices must be in bounds.

gr_mat_nrows(mat, ctx)

  Macro accessing the number of rows of mat.

gr_mat_ncols(mat, ctx)

  Macro accessing the number of columns of mat.

## Memory management

void gr_mat_init(gr_mat_t mat, slong rows, slong cols, gr_ctx_t ctx)

  Initializes mat to a matrix with the given number of rows and columns.

int gr_mat_init_set(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)

  Initializes res to a copy of the matrix mat.

void gr_mat_clear(gr_mat_t mat, gr_ctx_t ctx)

  Clears the matrix.

void gr_mat_swap(gr_mat_t mat1, gr_mat_t mat2, gr_ctx_t ctx)

  Swaps mat1 and mat12 efficiently.

int gr_mat_swap_entrywise(gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)

  Performs a deep swap of mat1 and mat2, swapping the individual entries rather than the top-level structures.

## Window matrices

void gr_mat_window_init(gr_mat_t window, const gr_mat_t mat, slong r1, slong c1, slong r2, slong c2, gr_ctx_t ctx)

  Initializes window to a window matrix into the submatrix of mat starting at the corner at row r1 and column c1 (inclusive) and ending at row r2 and column c2 (exclusive). The indices must be within bounds.

void gr_mat_window_clear(gr_mat_t window, gr_ctx_t ctx)

  Frees the window matrix.

GR_MAT_TMP_INIT_SHALLOW_TRANSPOSE(AT, A, ctx)

  Initialize AT to a shallow transpose of A for temporary use. This macro uses stack allocation if A is sufficiently small. The matrix AT can be used similarly to a window matrix for reading. It can also be used for writing, provided that one finishes the operation by performing `GR_MAT_SHALLOW_TRANSPOSE(A, AT, ctx)` to write back any changes to the shallow data.

GR_MAT_SHALLOW_TRANSPOSE(AT, A, ctx)

  Sets AT to the transpose of A, copying entries shallowly. Assumes that AT and A are not aliased (if aliased, the normal `gr_mat_transpose()` already uses shallow operations).

GR_MAT_TMP_CLEAR_SHALLOW_TRANSPOSE(AT, ctx)

  Free the shallow transpose allocated by `GR_MAT_TMP_INIT_SHALLOW_TRANSPOSE`.

## Input and output

int gr_mat_write(gr_stream_t out, const gr_mat_t mat, gr_ctx_t ctx)

  Write mat to the stream out.

int gr_mat_print(const gr_mat_t mat, gr_ctx_t ctx)

  Prints mat to standard output.

## Comparisons

truth_t gr_mat_equal(const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)

  Returns whether mat1 and mat2 are equal.

## Assignment and special values

truth_t gr_mat_is_zero(const gr_mat_t mat, gr_ctx_t ctx)
truth_t gr_mat_is_one(const gr_mat_t mat, gr_ctx_t ctx)
truth_t gr_mat_is_neg_one(const gr_mat_t mat, gr_ctx_t ctx)

  Returns whether mat respectively is the zero matrix or the scalar matrix with 1 or -1 on the main diagonal.

truth_t gr_mat_is_scalar(const gr_mat_t mat, gr_ctx_t ctx)

  Returns whether mat is a scalar matrix, being a diagonal matrix with identical elements on the main diagonal.

int gr_mat_zero(gr_mat_t res, gr_ctx_t ctx)

  Sets res to the zero matrix.

int gr_mat_one(gr_mat_t res, gr_ctx_t ctx)

  Sets res to the scalar matrix with 1 on the main diagonal and zero elsewhere.

int gr_mat_set(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_set_fmpz_mat(gr_mat_t res, const fmpz_mat_t mat, gr_ctx_t ctx)
int gr_mat_set_fmpq_mat(gr_mat_t res, const fmpq_mat_t mat, gr_ctx_t ctx)
int gr_mat_set_gr_mat_other(gr_mat_t res, const gr_mat_t mat, gr_ctx_t mat_ctx, gr_ctx_t ctx)

  Sets res to the value of mat.

int gr_mat_set_scalar(gr_mat_t res, gr_srcptr c, gr_ctx_t ctx)
int gr_mat_set_ui(gr_mat_t res, ulong c, gr_ctx_t ctx)
int gr_mat_set_si(gr_mat_t res, slong c, gr_ctx_t ctx)
int gr_mat_set_fmpz(gr_mat_t res, const fmpz_t c, gr_ctx_t ctx)
int gr_mat_set_fmpq(gr_mat_t res, const fmpq_t c, gr_ctx_t ctx)

  Set res to the scalar matrix with c on the main diagonal and zero elsewhere.

## Basic row, column and entry operations

int gr_mat_concat_horizontal(gr_mat_t res, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)

int gr_mat_concat_vertical(gr_mat_t res, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)

int gr_mat_transpose(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)

  Sets `res` to the transpose of `mat`. Dimensions must be compatible. Aliasing is allowed for square matrices.

int gr_mat_swap_rows(gr_mat_t mat, slong *perm, slong r, slong s, gr_ctx_t ctx)

  Swaps rows `r` and `s` of `mat`. If `perm` is non-`NULL`, the permutation of the rows will also be applied to `perm`.

int gr_mat_swap_cols(gr_mat_t mat, slong *perm, slong r, slong s, gr_ctx_t ctx)

  Swaps columns `r` and `s` of `mat`. If `perm` is non-`NULL`, the permutation of the columns will also be applied to `perm`.

int gr_mat_invert_rows(gr_mat_t mat, slong *perm, gr_ctx_t ctx)

  Swaps rows `i` and `r - i` of `mat` for `0 <= i < r/2`, where `r` is the number of rows of `mat`. If `perm` is non-`NULL`, the permutation of the rows will also be applied to `perm`.

int gr_mat_invert_cols(gr_mat_t mat, slong *perm, gr_ctx_t ctx)

  Swaps columns `i` and `c - i` of `mat` for `0 <= i < c/2`, where `c` is the number of columns of `mat`. If `perm` is non-`NULL`, the permutation of the columns will also be applied to `perm`.

int gr_mat_move_row(gr_mat_t A, slong i, slong new_i, gr_ctx_t ctx)

  Moves row `i` to the new position `new_i`, displacing all intervening rows. For example, with `i = 4` and `new_i = 7` this replaces the block of rows `A4, A5, A6, A7` by `A5, A6, A7, A4`. With `i = 7` and `new_i = 4` replaces `A4, A5, A6, A7` by `A7, A4, A5, A6`. Returns `GR_DOMAIN` if either index is not in bounds, otherwise is guaranteed to succeed.

truth_t gr_mat_is_empty(const gr_mat_t mat, gr_ctx_t ctx)

  Returns whether mat is an empty matrix, having either zero rows or zero column. This predicate is always decidable (even if the underlying ring is not computable), returning `T_TRUE` or `T_FALSE`.

truth_t gr_mat_is_square(const gr_mat_t mat, gr_ctx_t ctx)

  Returns whether mat is a square matrix, having the same number of rows as columns (not the same thing as being a perfect square!). This predicate is always decidable (even if the underlying ring is not computable), returning `T_TRUE` or `T_FALSE`.

## Entrywise operations

int gr_mat_entrywise_unary_op(gr_mat_t res, gr_method_unary_op f, const gr_mat_t mat, gr_ctx_t ctx)

  Sets res to the application of the function f to the entries of matrix mat. Returns `GR_DOMAIN` if the matrix dimensions do not match.

int gr_mat_entrywise_binary_op(gr_mat_t res, gr_method_binary_op f, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)

  Sets res to the application of the function f to the entries of mat1 as first argument and the entries of mat2 as second argument. Returns `GR_DOMAIN` if the matrix dimensions do not match.

int gr_mat_entrywise_binary_op_scalar(gr_mat_t res, gr_method_binary_op f, const gr_mat_t mat, gr_srcptr c, gr_ctx_t ctx)

  Sets res to the application of the function f to the entries of mat as first argument and the scalar c as second argument. Returns `GR_DOMAIN` if the matrix dimensions do not match.

truth_t gr_mat_entrywise_unary_predicate_all(gr_method_unary_predicate f, const gr_mat_t mat, gr_ctx_t ctx)
truth_t gr_mat_entrywise_unary_predicate_any(gr_method_unary_predicate f, const gr_mat_t mat, gr_ctx_t ctx)

  Returns whether the predicate f is true for all entries, respectively for any entry, in the matrix mat.

truth_t gr_mat_entrywise_binary_predicate_all(gr_method_binary_predicate f, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)

  Returns whether the binary predicate f is true for all entries in mat1 paired with the corresponding entries in mat2. Returns `T_FALSE` if the matrix dimensions are not compatible.

## Norms

int gr_mat_norm_max(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)

  Max norm: \(\max_{i,j} |a_{i,j}|\).

int gr_mat_norm_1(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)

  1-norm (largest absolute column sum): \(\max_{1\le j \le n} \sum_{i=1}^m |a_{i,j}|\).

int gr_mat_norm_inf(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)

  Infinity-norm (largest absolute row sum): \(\max_{1\le i \le m} \sum_{j=1}^n |a_{i,j}|\).

int gr_mat_norm_frobenius(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)

  Frobenius norm: \(\sqrt{\sum_{i,j} |a_{i,j}|^2}\).

## Addition and scalar arithmetic

int gr_mat_neg(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)

int gr_mat_add(gr_mat_t res, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)

int gr_mat_sub(gr_mat_t res, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)

int gr_mat_add_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t ctx)
int gr_mat_scalar_add(gr_mat_t res, gr_srcptr x, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_add_ui(gr_mat_t res, const gr_mat_t mat, ulong x, gr_ctx_t ctx)
int gr_mat_add_si(gr_mat_t res, const gr_mat_t mat, slong x, gr_ctx_t ctx)
int gr_mat_add_fmpz(gr_mat_t res, const gr_mat_t mat, const fmpz_t x, gr_ctx_t ctx)
int gr_mat_add_fmpq(gr_mat_t res, const gr_mat_t mat, const fmpq_t x, gr_ctx_t ctx)
int gr_mat_add_scalar_other(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t x_ctx, gr_ctx_t ctx)
int gr_mat_scalar_other_add(gr_mat_t res, gr_srcptr x, gr_ctx_t x_ctx, const gr_mat_t mat, gr_ctx_t ctx)

  Perform the matrix-scalar or scalar-matrix operation \(A + Ix\) or \(Ix + A\).

int gr_mat_sub_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t ctx)
int gr_mat_scalar_sub(gr_mat_t res, gr_srcptr x, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_sub_ui(gr_mat_t res, const gr_mat_t mat, ulong x, gr_ctx_t ctx)
int gr_mat_sub_si(gr_mat_t res, const gr_mat_t mat, slong x, gr_ctx_t ctx)
int gr_mat_sub_fmpz(gr_mat_t res, const gr_mat_t mat, const fmpz_t x, gr_ctx_t ctx)
int gr_mat_sub_fmpq(gr_mat_t res, const gr_mat_t mat, const fmpq_t x, gr_ctx_t ctx)
int gr_mat_sub_scalar_other(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t x_ctx, gr_ctx_t ctx)
int gr_mat_scalar_other_sub(gr_mat_t res, gr_srcptr x, gr_ctx_t x_ctx, const gr_mat_t mat, gr_ctx_t ctx)

  Perform the matrix-scalar or scalar-matrix operation \(A - Ix\) or \(Ix - A\).

int gr_mat_mul_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t ctx)
int gr_mat_scalar_mul(gr_mat_t res, gr_srcptr x, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_mul_ui(gr_mat_t res, const gr_mat_t mat, ulong x, gr_ctx_t ctx)
int gr_mat_mul_si(gr_mat_t res, const gr_mat_t mat, slong x, gr_ctx_t ctx)
int gr_mat_mul_fmpz(gr_mat_t res, const gr_mat_t mat, const fmpz_t x, gr_ctx_t ctx)
int gr_mat_mul_fmpq(gr_mat_t res, const gr_mat_t mat, const fmpq_t x, gr_ctx_t ctx)
int gr_mat_mul_scalar_other(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t x_ctx, gr_ctx_t ctx)
int gr_mat_scalar_other_mul(gr_mat_t res, gr_srcptr x, gr_ctx_t x_ctx, const gr_mat_t mat, gr_ctx_t ctx)

  Perform the matrix-scalar or scalar-matrix operation \(A x\) or \(x A\).

int gr_mat_div_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t ctx)
int gr_mat_div_scalar_other(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t x_ctx, gr_ctx_t ctx)
int gr_mat_div_ui(gr_mat_t res, const gr_mat_t mat, ulong x, gr_ctx_t ctx)
int gr_mat_div_si(gr_mat_t res, const gr_mat_t mat, slong x, gr_ctx_t ctx)
int gr_mat_div_fmpz(gr_mat_t res, const gr_mat_t mat, const fmpz_t x, gr_ctx_t ctx)
int gr_mat_div_fmpq(gr_mat_t res, const gr_mat_t mat, const fmpq_t x, gr_ctx_t ctx)

  Perform the matrix-scalar operation \(A / x\).

int gr_mat_addmul_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr c, gr_ctx_t ctx)
int gr_mat_submul_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr c, gr_ctx_t ctx)

## Matrix multiplication

int gr_mat_mul(gr_mat_t res, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)

  Compute \(AB\) using the default algorithm chosen by the element ring. If the element ring does not overload matrix multiplication, this will fall back to `gr_mat_mul_generic()` by default.

int gr_mat_mul_generic(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)

  Multiply matrices using a generic algorithm choice. Currently this always chooses classical multiplication, but may implement other strategies in the future.

int gr_mat_mul_classical(gr_mat_t res, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)

  Computes the \(m \times n \times p\) matrix product using the classical algorithm, performing \(mp\) dot products of length \(n\).

int gr_mat_mul_strassen(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)

  Uses Strassen’s algorithm to evaluate \(AB\) using 7 recursive matrix multiplications of roughly half the size. This function calls `gr_mat_mul()` for the recursive multiplications; to use Strassen recursively, the base ring must overload `gr_mat_mul()` to choose Strassen multiplication above some cutoff. This results in \(O(n^{2.81})\) asymptotic complexity in the case of an \(n \times n \times n\) product.

  This function does not implement Strassen’s original evaluation sequence but that of Bodrato [Bodrato2010] which uses fewer additions and offers some further time savings when squaring.

int gr_mat_mul_waksman(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)

  Compute the product \(AB\) using roughly half the number of multiplications of the classical algorithm, using Waksman’s algorithm [Waksman1970]. The base ring must be commutative and must support exact division by two.

int gr_mat_mul_rosowski(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)

  Compute the product \(AB\) using roughly half the number of multiplications of the classical algorithm, using Rosowski’s algorithm [Rosowski2023]. The base ring must be commutative. This uses the same number of multiplications as Waksman’s algorithm when the inner dimension n is even, but uses fewer multiplications when n is odd, and does not require division by two.

## Powering

See the section on matrix functions for non-integer exponents.

int gr_mat_sqr(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)

int gr_mat_pow_ui(gr_mat_t res, const gr_mat_t mat, ulong e, gr_ctx_t ctx)
int gr_mat_pow_si(gr_mat_t res, const gr_mat_t mat, slong e, gr_ctx_t ctx)
int gr_mat_pow_fmpz(gr_mat_t res, const gr_mat_t mat, const fmpz_t e, gr_ctx_t ctx)

## Polynomial evaluation

int _gr_mat_gr_poly_evaluate(gr_mat_t res, gr_srcptr poly, slong len, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_gr_poly_evaluate(gr_mat_t res, const gr_poly_t poly, const gr_mat_t mat, gr_ctx_t ctx)

  Sets res to the matrix obtained by evaluating the scalar polynomial poly with matrix argument mat.

## Diagonal and triangular matrices

truth_t gr_mat_is_upper_triangular(const gr_mat_t mat, gr_ctx_t ctx)
truth_t gr_mat_is_lower_triangular(const gr_mat_t mat, gr_ctx_t ctx)

  Returns whether mat is upper (respectively lower) triangular, having zeros everywhere below (respectively above) the main diagonal. The matrix need not be square.

truth_t gr_mat_is_diagonal(const gr_mat_t mat, gr_ctx_t ctx)

  Returns whether mat is a diagonal matrix, having zeros everywhere except on the main diagonal. The matrix need not be square.

int gr_mat_mul_diag(gr_mat_t res, const gr_mat_t A, const gr_vec_t D, gr_ctx_t ctx)
int gr_mat_diag_mul(gr_mat_t res, const gr_vec_t D, const gr_mat_t A, gr_ctx_t ctx)

  Set res to the product \(AD\) or \(DA\) respectively, where \(D\) is a diagonal matrix represented as a vector of entries.

## Gaussian elimination

int gr_mat_find_nonzero_pivot_large_abs(slong *pivot_row, gr_mat_t mat, slong start_row, slong end_row, slong column, gr_ctx_t ctx)
int gr_mat_find_nonzero_pivot_generic(slong *pivot_row, gr_mat_t mat, slong start_row, slong end_row, slong column, gr_ctx_t ctx)
int gr_mat_find_nonzero_pivot(slong *pivot_row, gr_mat_t mat, slong start_row, slong end_row, slong column, gr_ctx_t ctx)

  Attempts to find a nonzero element in column number column of the matrix mat in a row between start_row (inclusive) and end_row (exclusive). On success, sets `pivot_row` to the row index and returns `GR_SUCCESS`. If no nonzero pivot element exists, returns `GR_DOMAIN`. If no nonzero pivot element exists and zero-testing fails for some element, returns the flag `GR_UNABLE`.

  This function may be destructive: any elements that are nontrivially zero but can be certified zero may be overwritten by exact zeros.

int gr_mat_lu_classical(slong *rank, slong *P, gr_mat_t LU, const gr_mat_t A, int rank_check, gr_ctx_t ctx)
int gr_mat_lu_recursive(slong *rank, slong *P, gr_mat_t LU, const gr_mat_t A, int rank_check, gr_ctx_t ctx)
int gr_mat_lu_generic(slong *rank, slong *P, gr_mat_t LU, const gr_mat_t A, int rank_check, gr_ctx_t ctx)
int gr_mat_lu(slong *rank, slong *P, gr_mat_t LU, const gr_mat_t A, int rank_check, gr_ctx_t ctx)

  Computes a generalized LU decomposition \(A = PLU\) of a given matrix A, writing the rank of A to rank.

  If A is a nonsingular square matrix, LU will be set to a unit diagonal lower triangular matrix L and an upper triangular matrix U (the diagonal of L will not be stored explicitly).

  If A is an arbitrary matrix of rank r, U will be in row echelon form having r nonzero rows, and L will be lower triangular but truncated to r columns, having implicit ones on the r first entries of the main diagonal. All other entries will be zero.

  If a nonzero value for `rank_check` is passed, the function will abandon the output matrix in an undefined state and set the rank to 0 if A is detected to be rank-deficient. This currently only works as expected for square matrices.

  The algorithm can fail if it fails to certify that a pivot element is zero or nonzero, in which case the correct rank cannot be determined. It can also fail if a pivot element is not invertible. In these cases the `GR_UNABLE` and/or `GR_DOMAIN` flags will be returned. On failure, the data in the output variables `rank`, `P` and `LU` will be meaningless.

  The classical version uses iterative Gaussian elimination. The recursive version uses a block recursive algorithm to take advantage of fast matrix multiplication. The generic version calls the recursive algorithm with a default cutoff.

int gr_mat_fflu(slong *rank, slong *P, gr_mat_t LU, gr_ptr den, const gr_mat_t A, int rank_check, gr_ctx_t ctx)

  Similar to `gr_mat_lu()`, but computes a fraction-free LU decomposition using the Bareiss algorithm. The denominator is written to den.

## Solving

int gr_mat_nonsingular_solve_tril_classical(gr_mat_t X, const gr_mat_t L, const gr_mat_t B, int unit, gr_ctx_t ctx)
int gr_mat_nonsingular_solve_tril_recursive(gr_mat_t X, const gr_mat_t L, const gr_mat_t B, int unit, gr_ctx_t ctx)
int gr_mat_nonsingular_solve_tril_generic(gr_mat_t X, const gr_mat_t L, const gr_mat_t B, int unit, gr_ctx_t ctx)
int gr_mat_nonsingular_solve_tril(gr_mat_t X, const gr_mat_t L, const gr_mat_t B, int unit, gr_ctx_t ctx)
int gr_mat_nonsingular_solve_triu_classical(gr_mat_t X, const gr_mat_t U, const gr_mat_t B, int unit, gr_ctx_t ctx)
int gr_mat_nonsingular_solve_triu_recursive(gr_mat_t X, const gr_mat_t U, const gr_mat_t B, int unit, gr_ctx_t ctx)
int gr_mat_nonsingular_solve_triu_generic(gr_mat_t X, const gr_mat_t U, const gr_mat_t B, int unit, gr_ctx_t ctx)
int gr_mat_nonsingular_solve_triu(gr_mat_t X, const gr_mat_t U, const gr_mat_t B, int unit, gr_ctx_t ctx)

  Solves the lower triangular system \(LX = B\) or the upper triangular system \(UX = B\), respectively. Division by the the diagonal entries must be possible; if not a division fails, `GR_DOMAIN` is returned even if the system is solvable. If unit is set, the main diagonal of L or U is taken to consist of all ones, and in that case the actual entries on the diagonal are not read at all and can contain other data.

  The classical versions perform the computations iteratively while the recursive versions perform the computations in a block recursive way to benefit from fast matrix multiplication. The default versions choose an algorithm automatically.

int gr_mat_nonsingular_solve_fflu(gr_mat_t X, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
int gr_mat_nonsingular_solve_lu(gr_mat_t X, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
int gr_mat_nonsingular_solve(gr_mat_t X, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)

  Solves \(AX = B\). If A is not invertible, returns `GR_DOMAIN` even if the system has a solution.

int gr_mat_nonsingular_solve_fflu_precomp(gr_mat_t X, const slong *perm, const gr_mat_t LU, const gr_mat_t B, gr_ctx_t ctx)
int gr_mat_nonsingular_solve_lu_precomp(gr_mat_t X, const slong *perm, const gr_mat_t LU, const gr_mat_t B, gr_ctx_t ctx)

  Solves \(AX = B\) given a precomputed FFLU or LU factorization of A.

int gr_mat_nonsingular_solve_den_fflu(gr_mat_t X, gr_ptr den, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
int gr_mat_nonsingular_solve_den(gr_mat_t X, gr_ptr den, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)

  Solves \(AX = B\) over the fraction field of the present ring (assumed to be an integral domain), returning \(X\) with an implied denominator den. If A is not invertible over the fraction field, returns `GR_DOMAIN` even if the system has a solution.

int gr_mat_solve_field(gr_mat_t X, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)

  Solves \(AX = B\) where A is not necessarily square and not necessarily invertible. Assuming that the ring is a field, a return value of `GR_DOMAIN` indicates that the system has no solution. If there are multiple solutions, an arbitrary solution is returned.

## Determinant and trace

int gr_mat_det_fflu(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_det_berkowitz(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_det_lu(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_det_cofactor(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_det_generic_field(gr_ptr res, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_det_generic_integral_domain(gr_ptr res, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_det_generic(gr_ptr res, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_det(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)

  Sets res to the determinant of the square matrix mat. Various algorithms are available:

  - The berkowitz version uses the division-free Berkowitz algorithm performing \(O(n^4)\) operations. Since no zero tests are required, it is guaranteed to succeed if the ring arithmetic succeeds.

  - The cofactor version performs cofactor expansion. This is currently only supported for matrices up to size 4, and for larger matrices returns the `GR_UNABLE` flag.

  - The lu and fflu versions use rational LU decomposition and fraction-free LU decomposition (Bareiss algorithm) respectively, requiring \(O(n^3)\) operations. These algorithms can fail if zero certification or inversion fails, in which case the `GR_UNABLE` flag is returned.

  - The generic, generic_field and generic_integral_domain versions choose an appropriate algorithm for a generic ring depending on the availability of division.

  - The default method can be overloaded.

  If the matrix is not square, `GR_DOMAIN` is returned.

int gr_mat_trace(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)

  Sets res to the trace (sum of entries on the main diagonal) of the square matrix mat. If the matrix is not square, `GR_DOMAIN` is returned.

## Permanent

int gr_mat_permanent_cofactor(gr_ptr res, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_permanent_ryser(gr_ptr res, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_permanent_glynn(gr_ptr res, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_permanent_glynn_threaded(gr_ptr res, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_generic(gr_ptr res, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_permanent(gr_ptr res, const gr_mat_t A, gr_ctx_t ctx)

  Sets res to the permanent of the square matrix A, The permanent \(\operatorname{perm}(A)\) is a polynomial expression in the entries of \(A\) equivalent to that for the determinant but with all subtractions replaced by additions. Several algorithms are implemented:

  - The cofactor version uses recursive cofactor expansion, requiring \(O(n!)\) additions and multiplications.

  - The ryser version uses Ryser’s formula [Rys1963] with Gray code traversal, requiring \(O(2^n n)\) additions, subtractions and multiplications.

  - The glynn version uses Glynn’s formula [Gly2010] (actually equivalent to a method described earlier by Nijenhuis and Wilf [NW1978]). This requires about half as many operations as ryser but requires (exact) division by 2.

  - The glynn_threaded version is a multithreaded implementation of Glynn’s formula.

  The generic method chooses cofactor expansion for \(n \le 4\) and otherwise chooses ryser, glynn or glynn_threaded depending on whether the ring supports division by 2. The default method can be overloaded.

## Rank

int gr_mat_rank_fflu(slong *rank, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_rank_lu(slong *rank, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_rank(slong *rank, const gr_mat_t mat, gr_ctx_t ctx)

  Sets res to the rank of mat. The default method returns `GR_DOMAIN` if the element ring is not an integral domain, in which case the usual rank is not well-defined. The fflu and lu variants currently do not check the element domain, and simply return this flag if they encounter an impossible inverse in the execution of the respective algorithms.

## Row echelon form

int gr_mat_rref_lu(slong *rank, gr_mat_t R, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_rref_fflu(slong *rank, gr_mat_t R, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_rref(slong *rank, gr_mat_t R, const gr_mat_t A, gr_ctx_t ctx)

  Sets R to the reduced row echelon form of A, also setting rank to its rank.

int gr_mat_rref_den_fflu(slong *rank, gr_mat_t R, gr_ptr den, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_rref_den(slong *rank, gr_mat_t R, gr_ptr den, const gr_mat_t A, gr_ctx_t ctx)

  Like rref, but computes the reduced row echelon multiplied by a common (not necessarily minimal) denominator which is written to den. This can be used to compute the rref over an integral domain which is not a field.

## Nullspace

int gr_mat_nullspace(gr_mat_t X, const gr_mat_t A, gr_ctx_t ctx)

  Sets X to a basis for the (right) nullspace of A. On success, the output matrix will be resized to the correct number of columns.

  The basis is not guaranteed to be presented in a canonical or minimal form.

  If the ring is not a field, this is implied to compute a nullspace basis over the fraction field. The result may be meaningless if the ring is not an integral domain.

int gr_mat_nullspace_from_rref(gr_mat_t X, const gr_mat_t A, gr_srcptr Aden, slong rank, gr_ctx_t ctx)

  Computes nullspace given the precomputed reduced row echelon form matrix A with rank rank. If Aden is not NULL, assume tha A has been multiplied by this common denominator as in the output of `gr_mat_rref_den()`.

int gr_mat_nullspace_no_resize(slong *nullity, gr_mat_t X, const gr_mat_t A, gr_ctx_t ctx)

  Similar to `gr_mat_nullspace()`, but does not resize the matrix X, instead zero-padding if needed and returning the nullity (number of basis columns) in a separate variable. The user must supply a output matrix X with at least as many columns as the nullity.

## Inverse and adjugate

int gr_mat_inv(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)

  Sets res to the inverse of mat, computed by solving \(A A^{-1} = I\).

  Returns `GR_DOMAIN` if it can be determined that mat is not invertible over the present ring (warning: this may not work over non-integral domains). If invertibility cannot be proved, returns `GR_UNABLE`.

  To compute the inverse over the fraction field, one may use `gr_mat_nonsingular_solve_den()` or `gr_mat_adjugate()`.

int gr_mat_adjugate_charpoly(gr_mat_t adj, gr_ptr det, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_adjugate_cofactor(gr_mat_t adj, gr_ptr det, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_adjugate(gr_mat_t adj, gr_ptr det, const gr_mat_t mat, gr_ctx_t ctx)

  Sets adj to the adjugate matrix of mat, simultaneously setting det to the determinant of mat. We have \(\operatorname{adj}(A) A = A \operatorname{adj}(A) = \det(A) I\), and \(A^{-1} = \operatorname{adj}(A) / \det(A)\) when A is invertible.

  The cofactor version uses cofactor expansion, requiring the evaluation of \(n^2\) determinants. The charpoly version computes and then evaluates the characteristic polynomial, requiring \(O(n^{1/2})\) matrix multiplications plus \(O(n^3)\) or \(O(n^4)\) operations for the characteristic polynomial itself depending on the algorithm used.

## Characteristic polynomial

int _gr_mat_charpoly(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_charpoly(gr_poly_t res, const gr_mat_t mat, gr_ctx_t ctx)

  Computes the characteristic polynomial using an algorithm choice which defaults to `_gr_mat_charpoly_generic()` but may be overridden by specific rings for performance. The underscore method assumes that res is a preallocated array of \(n + 1\) coefficients.

int _gr_mat_charpoly_generic(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_charpoly_generic(gr_poly_t res, const gr_mat_t mat, gr_ctx_t ctx)

  Computes the characteristic polynomial using a generic algorithm choice.

int _gr_mat_charpoly_berkowitz(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_charpoly_berkowitz(gr_poly_t res, const gr_mat_t mat, gr_ctx_t ctx)

  Sets res to the characteristic polynomial of the square matrix mat, computed using the division-free Berkowitz algorithm. The number of operations is \(O(n^4)\) where n is the size of the matrix.

int _gr_mat_charpoly_danilevsky_inplace(gr_ptr res, gr_mat_t mat, gr_ctx_t ctx)
int _gr_mat_charpoly_danilevsky(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_charpoly_danilevsky(gr_poly_t res, const gr_mat_t mat, gr_ctx_t ctx)
int _gr_mat_charpoly_gauss(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_charpoly_gauss(gr_poly_t res, const gr_mat_t mat, gr_ctx_t ctx)
int _gr_mat_charpoly_householder(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_charpoly_householder(gr_poly_t res, const gr_mat_t mat, gr_ctx_t ctx)

  Sets res to the characteristic polynomial of the square matrix mat, computed using the Danilevsky algorithm, Hessenberg reduction using Gaussian elimination, and Hessenberg reduction using Householder reflections. The number of operations of each method is \(O(n^3)\) where n is the size of the matrix. The inplace version overwrites the input matrix.

  These methods require divisions and can therefore fail when the ring is not a field. They also require zero tests. The householder version also requires square roots. The flags `GR_UNABLE` or `GR_DOMAIN` are returned when an impossible division or square root is encountered or when a comparison cannot be performed.

int _gr_mat_charpoly_faddeev(gr_ptr res, gr_mat_t adj, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_charpoly_faddeev(gr_poly_t res, gr_mat_t adj, const gr_mat_t mat, gr_ctx_t ctx)
int _gr_mat_charpoly_faddeev_bsgs(gr_ptr res, gr_mat_t adj, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_charpoly_faddeev_bsgs(gr_poly_t res, gr_mat_t adj, const gr_mat_t mat, gr_ctx_t ctx)

  Sets res to the characteristic polynomial of the square matrix mat, computed using the Faddeev-LeVerrier algorithm. If the optional output argument adj is not NULL, it is set to the adjugate matrix, which is computed free of charge.

  The bsgs version uses a baby-step giant-step strategy, also known as the Preparata-Sarwate algorithm. This reduces the complexity from \(O(n^4)\) to \(O(n^{3.5})\) operations at the cost of requiring \(n^{0.5}\) temporary matrices to be stored.

  This method requires divisions by small integers and can therefore fail (returning the `GR_UNABLE` or `GR_DOMAIN` flags) in finite characteristic or when the underlying ring does not implement a division algorithm.

int _gr_mat_charpoly_from_hessenberg(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_charpoly_from_hessenberg(gr_poly_t res, const gr_mat_t mat, gr_ctx_t ctx)

  Sets res to the characteristic polynomial of the square matrix mat, which is assumed to be in Hessenberg form (this is currently not checked).

## Minimal polynomial

int gr_mat_minpoly_field(gr_poly_t res, const gr_mat_t mat, gr_ctx_t ctx)

  Compute the minimal polynomial of the matrix mat. The algorithm assumes that the coefficient ring is a field.

## Companion matrix

int _gr_mat_companion(gr_mat_t res, gr_srcptr poly, gr_ctx_t ctx)
int gr_mat_companion(gr_mat_t res, const gr_poly_t poly, gr_ctx_t ctx)

  Sets the n by n matrix res to the companion matrix of the polynomial poly which must have degree n. The underscore method reads \(n + 1\) input coefficients. The algorithm assumes that the leading coefficient of poly is invertible.

int _gr_mat_companion_fraction(gr_mat_t res_num, gr_ptr res_den, gr_srcptr poly, gr_ctx_t ctx)
int gr_mat_companion_fraction(gr_mat_t res_num, gr_ptr res_den, const gr_poly_t poly, gr_ctx_t ctx)

  Sets the n by n matrix res_num and the polynomial res_den so that the fraction is the companion matrix of the polynomial poly which must have degree n. The underscore method reads \(n + 1\) input coefficients.

## Similarity transformations

int gr_mat_apply_row_similarity(gr_mat_t M, slong r, gr_ptr d, gr_ctx_t ctx)

  Applies an elementary similarity transform to the \(n\times n\) matrix \(M\) in-place.

  If \(P\) is the \(n\times n\) identity matrix the zero entries of whose row \(r\) (\(0\)-indexed) have been replaced by \(d\), this transform is equivalent to \(M = P^{-1}MP\).

  Similarity transforms preserve the determinant, characteristic polynomial and minimal polynomial.

## Eigenvalues

int gr_mat_eigenvalues(gr_vec_t lambda, gr_vec_t mult, const gr_mat_t mat, int flags, gr_ctx_t ctx)
int gr_mat_eigenvalues_other(gr_vec_t lambda, gr_vec_t mult, const gr_mat_t mat, gr_ctx_t mat_ctx, int flags, gr_ctx_t ctx)

  Finds all eigenvalues of the given matrix in the ring defined by ctx, storing the eigenvalues without duplication in lambda (a vector with elements of type `ctx`) and the corresponding multiplicities in mult (a vector with elements of type `fmpz`).

  The interface is essentially the same as that of `gr_poly_roots()`; see its documentation for details.

int gr_mat_diagonalization_precomp(gr_vec_t D, gr_mat_t L, gr_mat_t R, const gr_mat_t A, const gr_vec_t eigenvalues, const gr_vec_t mult, gr_ctx_t ctx)
int gr_mat_diagonalization_generic(gr_vec_t D, gr_mat_t L, gr_mat_t R, const gr_mat_t A, int flags, gr_ctx_t ctx)
int gr_mat_diagonalization(gr_vec_t D, gr_mat_t L, gr_mat_t R, const gr_mat_t A, int flags, gr_ctx_t ctx)

  Computes a diagonalization \(LAR = D\) given a square matrix \(A\), where \(D\) is a diagonal matrix (returned as a vector) of the eigenvalues repeated according to their multiplicities, \(L\) is a matrix of left eigenvectors, and \(R\) is a matrix of right eigenvectors, normalized such that \(L = R^{-1}\). This implies that \(A = RDL = RDR^{-1}\). Either \(L\) or \(R\) (or both) can be set to `NULL` to omit computing the respective matrix.

  If the matrix has entries in a field then a return flag of `GR_DOMAIN` indicates that the matrix is non-diagonalizable over this field.

  The precomp version requires as input a precomputed set of eigenvalues with corresponding multiplicities, which can be computed with `gr_mat_eigenvalues()`.

## Jordan decomposition

int gr_mat_set_jordan_blocks(gr_mat_t mat, const gr_vec_t lambda, slong num_blocks, slong *block_lambda, slong *block_size, gr_ctx_t ctx)
int gr_mat_jordan_blocks(gr_vec_t lambda, slong *num_blocks, slong *block_lambda, slong *block_size, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_jordan_transformation(gr_mat_t mat, const gr_vec_t lambda, slong num_blocks, slong *block_lambda, slong *block_size, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_jordan_form(gr_mat_t J, gr_mat_t P, const gr_mat_t A, gr_ctx_t ctx)

## Matrix functions

int gr_mat_func_jordan(gr_mat_t res, const gr_mat_t A, gr_method_vec_op jet_func, gr_ctx_t ctx)
int gr_mat_func_param_jordan(gr_mat_t res, const gr_mat_t A, gr_method_vec_scalar_op jet_func, gr_srcptr c, gr_ctx_t ctx)

  Computes the matrix function \(f(A)\) using Jordan decomposition. The user supplies `int jet_func(gr_ptr r, gr_srcptr x, slong n, gr_ctx_t ctx)` which given a scalar \(x\) writes the jet \(f(x), f'(x), \ldots, f^{(n-1)} / (n-1)!\) to the array \(r\).

  The param version takes as input a function with signature `int jet_func(gr_ptr r, gr_srcptr x, slong n, gr_srcptr c, gr_ctx_t ctx)` for evaluating a function \(f(x, c)\) depending on an extra parameter \(c\). Although `c` is nominally passed as a `gr_srcptr`, it can be a void pointer to arbitrary data that `jet_func` knows how to handle.

int gr_mat_exp_jordan(gr_mat_t res, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_exp(gr_mat_t res, const gr_mat_t A, gr_ctx_t ctx)

int gr_mat_log_jordan(gr_mat_t res, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_log(gr_mat_t res, const gr_mat_t A, gr_ctx_t ctx)

int gr_mat_pow_scalar_jordan(gr_mat_t res, const gr_mat_t A, gr_srcptr c, gr_ctx_t ctx)
int gr_mat_pow_scalar(gr_mat_t res, const gr_mat_t A, gr_srcptr c, gr_ctx_t ctx)
int gr_mat_pow_fmpq_jordan(gr_mat_t res, const gr_mat_t mat, const fmpq_t exp, gr_ctx_t ctx)
int gr_mat_pow_fmpq(gr_mat_t res, const gr_mat_t mat, const fmpq_t exp, gr_ctx_t ctx)

  Compute \(A^c\) using Jordan decomposition. The non-Jordan methods also check for small integer exponents and delegate those to the standard powering method.

int gr_mat_sqrt(gr_mat_t res, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_rsqrt(gr_mat_t res, const gr_mat_t A, gr_ctx_t ctx)

  Compute a square root \(A^{1/2}\) or a reciprocal square root \(A^{-1/2}\). Currently the only implemented algorithm is the Jordan decomposition. Warning: this will often fail and return `GR_UNABLE` when the scalar type is not an algebraically closed field, even if the matrix is a perfect square.

## Hessenberg form

truth_t gr_mat_is_hessenberg(const gr_mat_t mat, gr_ctx_t ctx)

  Returns whether mat is in upper Hessenberg form.

int gr_mat_hessenberg_gauss(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_hessenberg_householder(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)
int gr_mat_hessenberg(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)

  Sets res to an upper Hessenberg form of mat. The gauss version uses Gaussian elimination. The householder version uses Householder reflections.

  These methods require divisions and zero testing and can therefore fail (returning `GR_UNABLE` or `GR_DOMAIN`) when the ring is not a field. The householder version additionally requires complex conjugation and the ability to compute square roots.

## Random matrices

int gr_mat_randtest(gr_mat_t res, flint_rand_t state, gr_ctx_t ctx)

  Sets res to a random matrix. The distribution is nonuniform.

int gr_mat_randops(gr_mat_t mat, flint_rand_t state, slong count, gr_ctx_t ctx)

  Randomises mat in-place by performing elementary row or column operations. More precisely, at most count random additions or subtractions of distinct rows and columns will be performed.

int gr_mat_randpermdiag(int *parity, gr_mat_t mat, flint_rand_t state, gr_ptr diag, slong n, gr_ctx_t ctx)

  Sets mat to a random permutation of the diagonal matrix with n leading entries given by the vector `diag`. Returns `GR_DOMAIN` if the main diagonal of `mat` does not have room for at least n entries. The parity (0 or 1) of the permutation is written to `parity`.

int gr_mat_randrank(gr_mat_t mat, flint_rand_t state, slong rank, gr_ctx_t ctx)

  Sets `mat` to a random sparse matrix with the given rank, having exactly as many non-zero elements as the rank. The matrix can be transformed into a dense matrix with unchanged rank by subsequently calling `gr_mat_randops()`.

  This operation only makes sense over integral domains (currently not checked).

int gr_mat_randsimilar(gr_mat_t mat, flint_rand_t state, slong opcount, gr_ctx_t ctx)

  Randomises mat in-place by conjugating by elementary row/column operations. More precisely, at most opcount conjugations by random elementary row/column operations will be performed.

## Orthogonal matrices

truth_t gr_mat_is_orthogonal(const gr_mat_t A, gr_ctx_t ctx)

  Returns whether A is an orthogonal matrix (orthonormal matrix), i.e. a square matrix satisfying \(A A^T = A^T A = I\). It is assumed (not checked) that the scalar ring is commutative.

truth_t gr_mat_is_row_orthogonal(const gr_mat_t A, gr_ctx_t ctx)
truth_t gr_mat_is_col_orthogonal(const gr_mat_t A, gr_ctx_t ctx)
truth_t gr_mat_is_row_orthonormal(const gr_mat_t A, gr_ctx_t ctx)
truth_t gr_mat_is_col_orthonormal(const gr_mat_t A, gr_ctx_t ctx)

  Returns whether A is of the following type:

  - Row-orthogonal: \(A A^T = D\) for some diagonal matrix D.

  - Column-orthogonal: \(A^T A = D\) for some diagonal matrix D.

  - Row-orthonormal: \(A A^T = I\).

  - Column-orthonormal: \(A^T A = I\).

int gr_mat_randtest_orthogonal(gr_mat_t A, flint_rand_t state, gr_ctx_t ctx)

  Generates a random orthogonal matrix. Uses Cayley’s construction, with a permutation matrix as a fallback. It is assumed (not checked) that the scalar ring is commutative. Fails with `GR_DOMAIN` if A is not square.

## QR decomposition

int gr_mat_lq_gso(gr_mat_t L, gr_mat_t Q, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_lq_recursive(gr_mat_t L, gr_mat_t Q, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_lq_generic(gr_mat_t L, gr_mat_t Q, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_lq(gr_mat_t L, gr_mat_t Q, const gr_mat_t A, gr_ctx_t ctx)
int gr_mat_qr(gr_mat_t Q, gr_mat_t R, const gr_mat_t A, gr_ctx_t ctx)

  Computes a QR or LQ decomposition. The \(A = QR\) decomposition orthogonalizes the columns of \(A\):

  - \(A\) is \(m \times n\) with \(m \ge n\) and full rank (\(A\) has \(n\) linearly independent columns)

  - \(Q\) is \(m \times n\) and column-orthonormal

  - \(R\) is \(n \times n\) and upper triangular

  The \(A = LQ\) decomposition is the transposed operation, which orthogonalizes the rows of \(A\):

  - \(A\) is \(m \times n\) with \(m \le n\) and full rank (\(A\) has \(m\) linearly independent rows)

  - \(L\) is \(m \times m\) and lower triangular

  - \(Q\) is \(m \times n\) and row-orthonormal

  For rectangular matrices, these definitions correspond to the “reduced”, “compact” or “economy-sized” QR or LQ decomposition. To construct a “full” QR decomposition where \(Q\) is a square, orthogonal matrix, one can compute the reduced QR decomposition, extend \(Q\) with its orthogonal complement, and pad \(R\) or \(L\) with zeros.

  The input matrix \(A\) should have elements in a subfield of \(\mathbb{R}\) or real floating-point entries. These functions may succeed for other types, but the Euclidean norm normalizations will not necessarily be meaningful.

  The gso algorithm uses Gram-Schmidt orthogonalization. The recursive algorithm uses block recursion. The QR decomposition is a simple wrapper around the LQ decomposition. We use the LQ decomposition internally as row operations are more efficient than column operations in the row-major format of `gr_mat_t`.

  Aliasing of \(A\) and \(Q\) is handled efficiently in-place. Aliasing of \(A\) and \(R\) or \(L\) is allowed, but requires a temporary internal copy.

  These methods return `GR_DOMAIN` if the matrix dimensions are not compatible, if \(A\) does not have full rank, or if the base ring does not support the necessary divisions or square roots to normalize vectors.

## Special matrices

For the following functions, the user supplies an output matrix with the intended number of rows and columns.

int gr_mat_ones(gr_mat_t res, gr_ctx_t ctx)

  Sets all entries in res to one.

int gr_mat_pascal(gr_mat_t res, int triangular, gr_ctx_t ctx)

  Sets res to a Pascal matrix, whose entries are binomial coefficients. If triangular is 0, constructs a full symmetric matrix with the rows of Pascal’s triangle as successive antidiagonals. If triangular is 1, constructs the upper triangular matrix with the rows of Pascal’s triangle as columns, and if triangular is -1, constructs the lower triangular matrix with the rows of Pascal’s triangle as rows.

int gr_mat_stirling(gr_mat_t res, int kind, gr_ctx_t ctx)

  Sets res to a Stirling matrix, whose entries are Stirling numbers. If kind is 0, the entries are set to the unsigned Stirling numbers of the first kind. If kind is 1, the entries are set to the signed Stirling numbers of the first kind. If kind is 2, the entries are set to the Stirling numbers of the second kind.

int gr_mat_hilbert(gr_mat_t res, gr_ctx_t ctx)

  Sets res to the Hilbert matrix, which has entries \(1/(i+j+1)\) for \(i, j \ge 0\).

int gr_mat_hadamard(gr_mat_t res, gr_ctx_t ctx)

  If possible, sets res to a Hadamard matrix of the provided size and returns `GR_SUCCESS`. Returns `GR_DOMAIN` if no Hadamard matrix of the given size exists, and `GR_UNABLE` if the implementation does not know how to construct a Hadamard matrix of the given size.

  A Hadamard matrix of size n can only exist if n is 0, 1, 2, or a multiple of 4. It is not known whether a Hadamard matrix exists for every size that is a multiple of 4. This function uses the Paley construction, which succeeds for all n of the form \(n = 2^e\) or \(n = 2^e (q + 1)\) where q is an odd prime power. Orders n for which Hadamard matrices are known to exist but for which this construction fails are 92, 116, 156, … (OEIS A046116).

## Helper functions for reduction

int gr_mat_reduce_row_generic(slong *column, gr_mat_t A, slong *P, slong *L, slong m, gr_ctx_t ctx)
int gr_mat_reduce_row(slong *column, gr_mat_t A, slong *P, slong *L, slong m, gr_ctx_t ctx)

  Reduce row n of the matrix \(A\), assuming the prior rows are in Gauss form. However those rows may not be in order. The entry \(i\) of the array \(P\) is the row of \(A\) which has a pivot in the \(i\)-th column. If no such row exists, the entry of \(P\) will be \(-1\). The function sets column to the column in which the \(n\)-th row has a pivot after reduction. This will always be chosen to be the first available column for a pivot from the left. This information is also updated in \(P\). Entry \(i\) of the array \(L\) contains the number of possibly nonzero columns of \(A\) row \(i\). This speeds up reduction in the case that \(A\) is chambered on the right. Otherwise the entries of \(L\) can all be set to the number of columns of \(A\). We require the entries of \(L\) to be monotonic increasing.

  By default the generic version is called; specific rings can overload this (typically to implement delayed canonicalisation).

## LLL

Let \(A = (a_0, \ldots, a_{n-1})\) be a set of linearly independent vectors over \(\mathbb{R}\) with Gram-Schmidt orthogonalization \((b_0, \ldots, b_{n-1})\) and Gram-Schmidt coefficients \(\mu_{i,j} = \langle a_i, b_i \rangle / \| b_j \|^2\). The basis \(A\) is said to be LLL-reduced with parameter (\(\delta\), \(\eta\)) where \(0.25 < \delta \le 1\) and \(0.5 \le \eta < \sqrt{\delta}\) if they satisfy the size reduction condition

\[|\mu_{i,j}| \le \eta, \quad 0 \le j < i < n\]

and the Lovász condition

\[(\delta - \mu_{i,i-1}^2) \| b_{i-1} \|_2^2 \le \| b_i \|_2^2, \quad 1 \le i \le n - 1.\]

truth_t gr_mat_is_row_lll_reduced_naive(const gr_mat_t A, gr_srcptr delta, gr_srcptr eta, gr_ctx_t ctx)
truth_t gr_mat_is_row_lll_reduced_with_removal_naive(const gr_mat_t A, gr_srcptr delta, gr_srcptr eta, gr_srcptr gs_B, slong newd, gr_ctx_t ctx)

  Check if the rows of A are LLL-reduced by naively performing the Gram-Schmidt orthogonalization and checking the conditions one row at a time.

  In interval arithmetic, these functions terminate eagerly: `T_UNKNOWN` is returned if the tests are inconclusive for one row, even if there is a possibility that a later row could prove that the result should be `T_FALSE`.

## Linear ODEs

int _gr_mat_gr_poly_solve_lode_newton_start(gr_mat_t Y, gr_mat_t Z, gr_poly_t A_denominator_inv, const gr_mat_t A_numerator, const gr_poly_t A_denominator, const gr_mat_t Y0, gr_ctx_t sol_poly_ctx)
int _gr_mat_gr_poly_solve_lode_newton_step(gr_mat_t Y, gr_mat_t Z, gr_poly_t A_denominator_inv, slong *len, const gr_mat_t A_numerator, const gr_poly_t A_denominator, int A_is_companion, gr_ctx_t sol_poly_ctx)
int _gr_mat_gr_poly_solve_lode_newton(gr_mat_t Y, gr_mat_t Z, const gr_mat_t A_numerator, const gr_poly_t A_denominator, const gr_mat_t Y0, slong len, gr_ctx_t A_poly_ctx, gr_ctx_t sol_poly_ctx)
int gr_mat_gr_poly_solve_lode_newton(gr_mat_t Y, const gr_mat_t A_numerator, const gr_poly_t A_denominator, const gr_mat_t Y0, slong len, gr_ctx_t A_poly_ctx, gr_ctx_t sol_poly_ctx)

  Solves the system of linear ordinary differential equations \(Y'(t) = A(t)Y(t)\) with rational function coefficients \(A(t)\) and initial condition \(Y(0) = Y_0\). Sets Y to the power series solution truncated to length len. The algorithm is a Newton iteration as described in [BCOSSS2007], but using an iteration for the power series inverse of A_denominator, and using a specialized implementation of multiplication by \(A(t)\) mod \(t^m\) in the case where \(A(t)\) is a companion matrix. An optimal sequence of precisions is used to reach the target len, even when it is not a power of two.

## Test functions

The following functions run iters test iterations, generating matrices up to size maxn. If ctx is set to `NULL`, a random ring is generated on each test iteration, otherwise the given ring is tested.

void gr_mat_test_mul(gr_method_mat_binary_op mul_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)

  Tests the given function `mul_impl` for correctness as an implementation of `gr_mat_mul()`.

void gr_mat_test_lu(gr_method_mat_lu_op lu_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)

  Tests the given function `mul_impl` for correctness as an implementation of `gr_mat_lu()`.

void gr_mat_test_det(gr_method_mat_unary_op_get_scalar det_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)

  Tests the given function `det_impl` for correctness as an implementation of `gr_mat_det()`.

void gr_mat_test_charpoly(gr_method_mat_unary_op_get_scalar charpoly_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)

  Tests the given function `charpoly_impl` for correctness as an implementation of `_gr_mat_charpoly()`.

void gr_mat_test_nonsingular_solve_tril(gr_method_mat_binary_op_with_flag solve_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
void gr_mat_test_nonsingular_solve_triu(gr_method_mat_binary_op_with_flag solve_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)

  Tests the given function `solve_impl` for correctness as an implementation of `gr_mat_nonsingular_solve_tril()` / `gr_mat_nonsingular_solve_triu()`.

void gr_mat_test_approx_mul_max_norm(gr_method_mat_binary_op mul_impl, gr_srcptr rel_tol, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)

  Tests the given implementation of matrix multiplication for accuracy over an approximate numerical ring by checking that \(|C-AB| \le |A||B| rel\_tol\) holds in the max norm, using classical multiplication for reference.

void gr_mat_test_approx_mul_pos_entrywise_accurate(gr_method_mat_binary_op mul_impl, gr_srcptr rel_tol, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)

  Tests the given implementation of matrix multiplication for accuracy over an approximate numerical ring by generating nonnegative matrices and checking that the entrywise relative error compared to classical multiplication does not exceed rel_tol.