flintlib / python-flint

Python bindings for Flint and Arb
MIT License
132 stars 27 forks source link

Make `fmpz_mat` preserve the domain in `.inv()` and `.rref()`. #198

Open oscarbenjamin opened 3 months ago

oscarbenjamin commented 3 months ago

This follows on from gh-109 which makes division (a / b) exact and gh-72 which ensures that poly roots methods return roots in the ground domain. As far as I know there are now only two places in python-flint where the domain is implicitly changed during an operation which are fmpz_mat.inv() and fmpz_mat.solve():

In [2]: M = fmpz_mat([[1, 2], [3, 4]])

In [3]: M.inv()
Out[3]: 
[ -2,    1]
[3/2, -1/2]

In [4]: type(_)
Out[4]: flint.types.fmpq_mat.fmpq_mat

In [16]: M2 = fmpz_mat([[3, 5], [6, 7]])

In [17]: M.solve(M2)
Out[17]: 
[  0, -3]
[3/2,  4]

Another inconsistency is found with fmpz_mat.rref() where most types have an .rref() method that returns a 2-tuple (RREF, rank) the fmpz_mat.rref() method returns a 3-tuple (RREF, denominator, rank) from fraction-free RREF. This preserves the domain but not the interface of the .rref() method:

In [5]: M.rref()
Out[5]: 
([-2,  0]
 [ 0, -2],
 -2,
 2)

In [6]: fmpq_mat(M).rref()
Out[6]: 
([1, 0]
 [0, 1],
 2)

Currently fmpz_mat is the only matrix type where the ground domain is an integral-domain but not a field (all other types are fields except nmod/fmpz_mod with non-prime modulus). That means that right now fmpz_mat is the only type for which it makes sense to have fraction-free methods but ideally in future there will be more e.g. matrices with polynomial entries, gr_mat etc.

We need a better system that ensures that all types can have consistent interfaces and preserve the domain. The way I handled this with SymPy's DomainMatrix is very similar to the way it is handled in Flint's gr_mat generic type and uses the exact same naming convention (I don't think I had seen the gr_mat function names at the time...):

Here the first set like inv, rref etc are defined for fields and the second set are defined for all integral domains. Both sets always preserve the domain.

We should add methods inv_den, rref_den, etc to all matrix types (they could still raise for non-prime modulus). The current fmpz_mat.rref() method should just be renamed to rref_den().

Then fmpz_mat.rref() should be a function that either returns a 2-tuple (fmpz_mat, rank) or raises an error. The implementation is just:

def rref(M):
    (Mrref, den, rank) = M.rref_den()
    # Mrref / den raises DomainError if division is not exact.
    return (Mrref / den, rank)

Given python-flint's use of exact division I think that what makes most sense is that fmpz_mat.inv() should return an fmpz_mat if possible and raise an error otherwise (if the matrix is not unimodular). Again the implementation is simple:

def inv(M):
    (Minv, den) = M.inv_den()
    return Minv / den  # raises if not exact

The gr_mat functions for some of these operations are:

If I understand the docs correctly these all do exactly what I described above for the corresponding method names. There is no gr_mat_inv_den but the docs suggest using gr_mat_nonsingular_solve_den and I think it would be worth having a method called inv_den in python-flint since the implementation is trivial:

def inv_den(M):
    return M.solve_den(eye(M.nrows()))

I'm not sure about the name nonsingular_solve vs solve. The convention in python-flint is already solve for fmpz_mat, fmpq_mat, arb_mat and acb_mat. This is like the Flint functions fmpz_mat_solve and fmpq_mat_solve etc whose meaning corresponds either to solve or solve_den as described above. The gr_mat methods are all called nonsingular_solve or nonsingular_solve_den with the exception of gr_solve_field which does not require an invertible matrix.

The other method names relevant here but not used yet in python-flint are lu and fflu which is analogous to e.g. rref vs rref_den. I think the same convention could be used in python-flint for LU decompositions.

What this all means for now is that we need to change fmpz_mat.inv, fmpz_mat.solve and fmpz_mat.rref in an incompatible way for consistency with current and future matrix types. All three methods need to be changed to return fmpz_mat always and to raise an error if the result cannot be represented as an fmpz_mat.

I don't think anything special needs to be done to ease the impact of compatibility breakage because the old behaviour is available trivially by converting to fmpq_mat e.g.

M.inv() -> fmpq_mat(M).inv()
M.solve(b) -> fmpq_mat(M).solve(fmpq_mat(b))

or by just calling the new method name:

M.rref() -> M.rref_den()

If someone needs to straddle both new and old versions of python-flint then converting to fmpq_mat already works in older versions. In the RREF case the can do something like:

rref_den = getattr(fmpz_mat, "rref_den", fmpz_mat.rref)
M.rref() -> rref_den(M)