Add Flint matrices for ZZ and QQ

[oscarbenjamin]

References to other Issues or PRs

Followup from gh-25474

Built on top of gh-25492 (which could be merged first or otherwise this could be merged instead).

See also gh-25410.

Brief description of what is fixed or changed

Uses python-flint's fmpz and fmpq if SYMPY_GROUND_TYPES=flint for dense matrices over ZZ and QQ.

Other comments

Plenty more work can be done to make this faster but first step is to get it working.

Release Notes

• matrices
  • The dense implementation of DomainMatrix will now use python-flint matrices for matrices over ZZ and QQ if SYMPY_GROUND_TYPES=flint and python-flint is installed. A new low-level matrix type DFM is added to be a wrapper for dense flint matrices. DomainMatrix will use this internally by default for dense matrices when the ground types are flint. This makes many operations with dense DomainMatrix much faster although for now to test this you need to explicitly choose to use the dense implementation by doing e.g. M = Matrix([[1, 2], [3, 4]]).to_DM().to_dense() to get a matrix that has a flint matrix as its internal representation.

[sympy-bot]

✅

Hi, I am the SymPy bot (v170). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

• matrices
  • The dense implementation of DomainMatrix will now use python-flint matrices for matrices over ZZ and QQ if SYMPY_GROUND_TYPES=flint and python-flint is installed. A new low-level matrix type DFM is added to be a wrapper for dense flint matrices. DomainMatrix will use this internally by default for dense matrices when the ground types are flint. This makes many operations with dense DomainMatrix much faster although for now to test this you need to explicitly choose to use the dense implementation by doing e.g. M = Matrix([[1, 2], [3, 4]]).to_DM().to_dense() to get a matrix that has a flint matrix as its internal representation. (#25495 by @oscarbenjamin)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.13.

Click here to see the pull request description that was parsed.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->

Followup from gh-25474

Built on top of gh-25492 (which could be merged first or otherwise this could be merged instead).

See also gh-25410.

#### Brief description of what is fixed or changed

Uses python-flint's fmpz and fmpq if `SYMPY_GROUND_TYPES=flint` for dense matrices over ZZ and QQ.

#### Other comments

Plenty more work can be done to make this faster but first step is to get it working.

#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* matrices
    * The dense implementation of DomainMatrix will now use python-flint matrices for matrices over ZZ and QQ if `SYMPY_GROUND_TYPES=flint` and python-flint is installed. A new low-level matrix type `DFM` is added to be a wrapper for dense flint matrices. `DomainMatrix` will use this internally by default for dense matrices when the ground types are `flint`. This makes many operations with dense DomainMatrix much faster although for now to test this you need to explicitly choose to use the dense implementation by doing e.g. `M = Matrix([[1, 2], [3, 4]]).to_DM().to_dense()` to get a matrix that has a flint matrix as its internal representation.
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[sympy-bot]

🟠

Hi, I am the SymPy bot (v170). I've noticed that some of your commits add or delete files. Since this is sometimes done unintentionally, I wanted to alert you about it.

This is an experimental feature of SymPy Bot. If you have any feedback on it, please comment at sympy/sympy-bot#75.

The following commits add new files:

• ea3d527:
  • requirements-test.txt
• 2ca8816:
  • sympy/polys/matrices/_dfm.py
  • sympy/polys/matrices/dfm.py
• 86e585c:
  • sympy/polys/matrices/tests/test_xxm.py

If these files were added/deleted on purpose, you can ignore this message.

[oscarbenjamin]

At this point the existing tests should pass although I still need to add docstrings and new tests for DFM. Also for now at least the DFM methods mostly convert to the existing pure Python implementations so they do not use the flint routines (I will change that once tests are added) so you should call methods on the underlying flint matrix directly to test timings.

To test the effect of this PR you can do:

# works on OSX/Windows, otherwise see docs for building
pip install python-flint

Then you need to set SYMPY_GROUND_TYPES=flint before importing sympy e.g.:

SYMPY_GROUND_TYPES=flint isympy

Or otherwise from Python:

import os
os.environ['SYMPY_GROUND_TYPES'] = 'flint'
import sympy

Now if you make a dense DomainMatrix over ZZ or QQ (integers or rationals) then its internal representation will be DFM which is a wrapper around a python-flint fmpz_mat or fmpq_mat:

In [1]: M = randMatrix(5)

In [2]: M
Out[2]: 
⎡98  81  49  15  8 ⎤
⎢                  ⎥
⎢60  26  73  51  60⎥
⎢                  ⎥
⎢52  8   3   48  22⎥
⎢                  ⎥
⎢9   45  37  3   76⎥
⎢                  ⎥
⎣25  49  53  8   39⎦

In [3]: M.to_DM()
Out[3]: DomainMatrix({0: {0: 98, 1: 81, 2: 49, 3: 15, 4: 8}, 1: {0: 60, 1: 26, 2: 73, 3: 51, 4: 60}, 2: {0: 52, 1: 8, 2: 3, 3: 48, 4: 22}, 3: {0: 9, 1: 45, 2: 37, 3: 3, 4: 76}, 4: {0: 25, 1: 49, 2: 53, 3: 8, 4: 39}}, (5, 5), ZZ)

In [4]: M.to_DM().to_dense()
Out[4]: DomainMatrix([[98, 81, 49, 15, 8], [60, 26, 73, 51, 60], [52, 8, 3, 48, 22], [9, 45, 37, 3, 76], [25, 49, 53, 8, 39]], (5, 5), ZZ)

In [5]: M.to_DM().to_dense().rep
Out[5]: DFM([[98, 81, 49, 15, 8], [60, 26, 73, 51, 60], [52, 8, 3, 48, 22], [9, 45, 37, 3, 76], [25, 49, 53, 8, 39]], (5, 5), ZZ)

In [6]: M.to_DM().to_dense().rep.rep
Out[6]: 
[98, 81, 49, 15,  8]
[60, 26, 73, 51, 60]
[52,  8,  3, 48, 22]
[ 9, 45, 37,  3, 76]
[25, 49, 53,  8, 39]

In [7]: type(M.to_DM().to_dense().rep.rep)
Out[7]: flint._flint.fmpz_mat

In [8]: M.to_DM().to_dense().rep.rep.charpoly()
Out[8]: x^5 + (-169)*x^4 + (-8741)*x^3 + 115920*x^2 + 22323182*x + 250829643

In [9]: type(M.to_DM().to_dense().rep.rep.charpoly())
Out[9]: flint._flint.fmpz_poly

Here are some comparative timings:

In [13]: for n in [10, 50, 100]:
    ...:     M = randMatrix(n)
    ...:     dM = M.to_DM()
    ...:     fM = dM.to_dense().rep.rep
    ...:     print(f'\ntime for {n}x{n} charpoly with matrix of small integers')
    ...:     print('Matrix:')
    ...:     %time M.charpoly()
    ...:     print('DomainMatrix:')
    ...:     %time dM.charpoly()
    ...:     print('Flint:')
    ...:     %time fM.charpoly()
    ...: 

time for 10x10 charpoly with matrix of small integers
Matrix:
CPU times: user 30 ms, sys: 117 µs, total: 30.1 ms
Wall time: 30.2 ms
DomainMatrix:
CPU times: user 0 ns, sys: 1.84 ms, total: 1.84 ms
Wall time: 1.67 ms
Flint:
CPU times: user 0 ns, sys: 75 µs, total: 75 µs
Wall time: 78.4 µs

time for 50x50 charpoly with matrix of small integers
Matrix:
CPU times: user 2.56 s, sys: 8 µs, total: 2.56 s
Wall time: 2.56 s
DomainMatrix:
CPU times: user 841 ms, sys: 0 ns, total: 841 ms
Wall time: 841 ms
Flint:
CPU times: user 6.13 ms, sys: 0 ns, total: 6.13 ms
Wall time: 6.13 ms

time for 100x100 charpoly with matrix of small integers
Matrix:
CPU times: user 37.1 s, sys: 316 ms, total: 37.4 s
Wall time: 37.4 s
DomainMatrix:
CPU times: user 18.9 s, sys: 35.9 ms, total: 19 s
Wall time: 19 s
Flint:
CPU times: user 144 ms, sys: 0 ns, total: 144 ms
Wall time: 144 ms

So DomainMatrix is a lot faster than Matrix although for ZZ the difference is not that great. See #25458 (comment) for more comparisons of Matrix vs DomainMatrix. That PR will speed up Matrix.charpoly() by using DomainMatrix which is very significant for QQ and for other domains. In the case of ZZ or QQ though we can see that flint is a lot faster for charpoly. I believe that the big-O is not the same here but I don't know exactly what the algorithmic differences are. Scaling up with flint we have:

In [14]: M = randMatrix(200)

In [15]: %time ok = M.to_DM().to_dense().rep.rep.charpoly()
CPU times: user 1.33 s, sys: 4.03 ms, total: 1.33 s
Wall time: 1.33 s

In [18]: M = randMatrix(300)

In [19]: %time ok = M.to_DM().to_dense().rep.rep.charpoly()
CPU times: user 6.87 s, sys: 35.8 ms, total: 6.91 s
Wall time: 6.89 s

In [20]: M = randMatrix(400)

In [21]: %time ok = M.to_DM().to_dense().rep.rep.charpoly()
CPU times: user 22.2 s, sys: 4.04 ms, total: 22.3 s
Wall time: 22.3 s

[github-actions]

Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with + are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with - are speedups.

Significantly changed benchmark results (PR vs master)

Significantly changed benchmark results (master vs previous release)

       before           after         ratio
     [8059df73]       [fdce9ca4]
     <sympy-1.12^0>                 
-      87.8±0.8ms       56.8±0.5ms     0.65  integrate.TimeIntegrationRisch02.time_doit(10)
-        86.3±2ms       55.3±0.4ms     0.64  integrate.TimeIntegrationRisch02.time_doit_risch(10)
-        2.11±0ms        658±0.8μs     0.31  polys.TimePREM_LinearDenseQuadraticGCD.time_op(3, 'sparse')
-     10.4±0.04ms      1.96±0.01ms     0.19  polys.TimePREM_LinearDenseQuadraticGCD.time_op(5, 'sparse')
-       366±0.7μs       82.1±0.2μs     0.22  polys.TimePREM_QuadraticNonMonicGCD.time_op(1, 'sparse')
-     4.86±0.01ms        364±0.8μs     0.07  polys.TimePREM_QuadraticNonMonicGCD.time_op(3, 'sparse')
-     10.7±0.04ms         1.10±0ms     0.10  polys.TimePREM_QuadraticNonMonicGCD.time_op(5, 'sparse')
-     6.30±0.01ms      3.96±0.01ms     0.63  polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(2, 'sparse')
-     27.3±0.08ms      12.1±0.03ms     0.44  polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(3, 'sparse')
-     6.97±0.03ms         1.17±0ms     0.17  polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(1, 'sparse')
-     16.2±0.02ms      9.23±0.03ms     0.57  polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(2, 'sparse')
-       212±0.3ms      70.8±0.09ms     0.33  polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(3, 'sparse')
-     6.49±0.02ms          534±1μs     0.08  polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(3, 'sparse')
-     27.9±0.03ms          857±3μs     0.03  polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(5, 'sparse')
-         616±2μs          198±2μs     0.32  polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(1, 'sparse')
-     6.56±0.02ms        203±0.9μs     0.03  polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(3, 'sparse')
-     17.2±0.03ms          206±1μs     0.01  polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(5, 'sparse')
-       164±0.4μs       88.5±0.2μs     0.54  solve.TimeMatrixOperations.time_rref(3, 0)
-       313±0.5μs          110±4μs     0.35  solve.TimeMatrixOperations.time_rref(4, 0)
-     31.7±0.05ms      13.5±0.05ms     0.43  solve.TimeSolveLinSys189x49.time_solve_lin_sys

Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).

[oscarbenjamin]

In the last commit I have changed it so that DFM.det calls through to the FLINT implementation of either fmpz_mat_det or fmpq_mat_det.

There is no sparse implementation for SDM.det which just converts to dense and calls DDM.det. I have changed that here to convert to DFM if possible so that now sparse domain matrices will end up using the Flint det method for matrices over ZZ or QQ. That makes computing determinants much faster:

In [3]: M = randMatrix(1000) # 1000x1000 matrix

In [4]: %time ok = M.to_DM().det() # Using ground types = flint
CPU times: user 18.9 s, sys: 74.5 ms, total: 19 s
Wall time: 19 s

[oscarbenjamin]

The last commit uses python-flint for the inverse of matrices over QQ in DFM. As for det the SDM method just calls the DDM method and so is changed here to call the DFM method if possible.

This makes computing the inverse of DomainMatrix over QQ much faster:

In [19]: M = randMatrix(200) # 200x200 matrix

In [20]: %time ok = M.to_DM().to_dense().to_field().inv()
CPU times: user 3.41 s, sys: 15.9 ms, total: 3.42 s
Wall time: 3.41 s

Without flint (but with gmpy2) the existing DomainMatrix implementation gives

In [7]: M = randMatrix(200) # 200x200 matrix

In [8]: %time ok = M.to_DM().to_dense().to_field().inv()
CPU times: user 1min 57s, sys: 51.1 ms, total: 1min 57s
Wall time: 1min 57s

The difference is firstly algorithmic because fraction-free methods are better in this case so DomainMatrix inv_den is a lot faster:

In [10]: %time ok = M.to_DM().inv_den()
CPU times: user 14.2 s, sys: 0 ns, total: 14.2 s
Wall time: 14.2 s

Compared to Flint that is about 4x slower which is consistent with other comparisons of Flint vs DomainMatrix fraction-free methods. I believe that 4x difference is primarily due to some algorithmic difference in solving linear systems like maybe Flint uses LU solve whereas DomainMatrix always uses RREF for solving.

[oscarbenjamin]

The last commit enables using flint for charpoly if possible. Same story as before: there is no sparse implementation and so SDM uses DDM.charpoly. Here that is changed to use DFM.charpoly if possible which in turn calls into the Flint routine.

Note that I have a PR gh-25458 which adds a sparse implementation of charpoly. The two PRs are compatible but the method for selecting an algorithm in gh-25458 should probably always choose Flint when possible. The Flint algorithm is modular with CRT reconstruction and is significantly better than any algorithms currently in SymPy.

Timing:

In [1]: M = randMatrix(200) # 200x200 matrix

In [2]: %time ok = M.to_DM().charpoly()
CPU times: user 1.49 s, sys: 81 µs, total: 1.49 s
Wall time: 1.49 s

Without Flint that example is about 60x slower:

In [11]: M = randMatrix(200) # 200x200 matrix

In [12]: %time ok = M.to_DM().charpoly()
CPU times: user 1min 26s, sys: 43.5 ms, total: 1min 26s
Wall time: 1min 26s

The difference is algorithmic though so I expect the ratio to grow for larger matrices.

[oscarbenjamin]

This is ready for review now.

I was originally thinking that the much simpler PR gh-25492 might be merged first just because it is simpler and sort of independent of the other changes here. It doesn't matter too much but I have already needed to fix merge conflicts in gh-25492 so it would be good to get that in rather than leaving it open.

This PR adds a DFM class to parallel the existing DDM and SDM classes that are the dense and sparse internal implementations of DomainMatrix. The DFM class wraps a Flint matrix which for now is either an fmpz_mat (a dense integer matrix) or an fmpq_mat (a dense rational matrix). Since Flint is implemented in C and has good algorithms for all of the operations used here this can bring some significant speedups.

This is still just a start in terms of making use of python-flint in SymPy. I could open many more PRs like this but I have a lot of open PRs right now so I'm going to wait for some to get reviewed or merged before opening any more.

To test the PR you should install python-flint:

pip install python-flint

and then set the SYMPY_GROUND_TYPES=flint environment variable as described above. Here are some timings.

Matrix multiplication:

In [25]: M = randMatrix(100) # 100x100 matrix of smallish integers

In [26]: dM = M.to_DM().to_dense() # Convert to a dense DomainMatrix

In [27]: type(dM.rep) # The DFM type added in this PR
Out[27]: sympy.polys.matrices._dfm.DFM

In [28]: type(dM.rep.rep) # The underlying Flint type
Out[28]: flint._flint.fmpz_mat

In [29]: %time ok = dM * dM # time to multiply two 100x100 Flint matrices
CPU times: user 4.64 ms, sys: 4 µs, total: 4.64 ms
Wall time: 4.79 ms

In [30]: %time ok = M * M # time to multiply with SymPy pure Python implementation
CPU times: user 458 ms, sys: 4.01 ms, total: 462 ms
Wall time: 461 ms

So the difference here is that Flint can multiply two dense integer matrices about 100x faster than SymPy's existing implementation. Note that for this operation DomainMatrix is used internally so this really comparing DomainMatrix against Flint. The factor of 100 difference is typical of optimised Python code vs optimised C code. For larger matrices the difference is larger and then the difference is algorithmic:

In [31]: M = randMatrix(200) # 200x200 matrix of smallish integers

In [32]: dM = M.to_DM().to_dense() # Convert to a dense DomainMatrix

In [33]: %time ok = dM * dM # time to multiply two 200x200 Flint matrices
CPU times: user 11.6 ms, sys: 0 ns, total: 11.6 ms
Wall time: 11.5 ms

In [34]: %time ok = M * M # time to multiply with SymPy pure Python implementation
CPU times: user 3.38 s, sys: 8 µs, total: 3.38 s
Wall time: 3.38 s

Here for 200x200 matrices we see that Flint is about 300x faster. The difference in ratio here is because the underlying algorithm is different (I'm not sure exactly what the algorithmic difference is yet). We can expect that for large matrices that difference in timing will grow.

This PR means that in principle SymPy could use Flint for operations like multiplication although the PR does not actually enable that. Deciding whether or not to use Flint for a given operation is complicated because the existing SymPy implementation for things like this uses a sparse DomainMatrix representation which is optimised for sparse matrices. If we make the matrices more and more sparse then at some point the existing sparse DomainMatrix implementation will be faster than Flint's dense representation:

In [35]: M = randMatrix(200, percent=1) # 200x200 matrix with 1% entries nonzero

In [36]: dM = M.to_DM().to_dense() # Convert to a dense DomainMatrix

In [37]: %time ok = dM * dM # time to multiply two 200x200 sparse Flint matrices
CPU times: user 24.3 ms, sys: 2 µs, total: 24.3 ms
Wall time: 24.1 ms

In [38]: %time ok = M * M # time to multiply with SymPy pure Python implementation
CPU times: user 2.67 ms, sys: 0 ns, total: 2.67 ms
Wall time: 2.69 ms

The situation in master right now is that SymPy's current sparse representation is always faster than its dense representation so we do not need to have any code that would decide which representation to use: just always use the sparse representation. After this PR the situation would be more complicated and we would need to have code that checks exactly how sparse the matrix is to decide whether to use the sparse implementation or to use Flint if it is available. I have not added that here so as of this PR the sparse representation is still used by default except in cases that already convert to the dense representation (e.g. because a sparse implementation of an algorithm is not implemented such as the operations mentioned below).

In future the approach here could be used to add support for many more types of matrices but for now this only applies to matrices with integer or rational entries. If SYMPY_GROUND_TYPES=flint then this PR immediately speeds up a number of operations for DomainMatrix:

• det
• charpoly
• inv
• lu_solve
• lll

The differences can be very significant because Flint is faster both as a result of being implemented in C rather than Python and also because it generally uses the best algorithms (at least better than many of SymPy's existing algorithms). In some cases the SymPy algorithms could be improved to give performance that more closely resembles the Flint performance. In other cases like the matrix multiplication example above the Python overhead is dominant and a pure Python implementation could never become close in speed to a C implementation. In any case Flint is the state of the art for these kinds of exact/symbolic calculations so being able to use them makes SymPy as fast as anything else. Note that other computer algebra systems (Mathematica etc) already use Flint for exactly this reason.

There are more operations that are implemented in python-flint that I have not immediately added in this PR for various reasons. The PR is already large and for some of the operations that it would be good to add the best thing would be to make some nontrivial changes to the existing DomainMatrix classes. I would prefer to do that in separate PRs that can be reviewed independently more easily then the changes here.

Some examples of things that I have not implemented and reasons are:

• Reduced row echelon form (RREF). The flint matrices do support RREF and it is fast. However we now have a more complicated method for selecting RREF algorithms and so if we add support for Flint's RREF algorithms then it needs to be incorporated into that. Having worked on this PR I am now less happy with the compromises I agreed to in perf(matrices): Use DomainMatrix for Matrix.rref over ZZ/QQ #25443 but I would prefer to discuss that in a separate PR from this one.
• Smith normal form (SNF) and Hermite normal form (HNF). The flint integer matrices have this and it is fast. However the existing DDM and SDM types do not have an snf method and the existing code works only with DomainMatrix as input rather than DDM. This is not difficult to change but it requires a small refactor that I want to avoid in this PR.
• Nullspace. There are some problems with the signature of the DomainMatrix, DDM and SDM nullspace methods that make it difficult to use Flint here. That is not hard to fix but again requires a small refactor that I want to avoid in this PR. Also in python-flint only matrices over ZZ have a nullspace method and not matrices over QQ which is the reverse of the way that DomainMatrix currently works. The nullspace method should be added for QQ in python-flint. Also the DDM and SDM nullspace methods should be changed to work over rings rather than just fields.
• Certain DomainMatrix (and DDM and SDM) methods require a field and will raise an error otherwise e.g. inv, solve, rref, lu, lu_solve. The corresponding python-flint matrix methods will just convert to a field so e.g. fmpz_mat.inv will just return a fmpq_mat. The fmpz_mat.inv method has an integer=True option that will compute the inverse of an integer matrix as an integer matrix if possible (if the matrix is unimodular) and raise otherwise. I think that I prefer the python-flint approach here and DomainMatrix should generally move to that instead. Again this is a change that I don't want to make in this PR so that it can be discussed separately.

I have also added a test module that uses pytest's parametrised test functionality to run the same test suite across all of SDM, DDM and the new DFM introduced here. Running the same tests on all implementations reveals some inconsistencies. Since these types are supposed to be interchangeable I have made some small changes to the existing types to improve compatibility but I have deliberately left out cases where larger changes would be needed or new methods would need to be added to the existing types (e.g. hnf or snf).

In future work a lot can be done to expose more of Flint's functionality in python-flint and expose more of python-flint's functionality in DomainMatrix. The final step in each case is exposing DomainMatrix's functionality and performance in Matrix which I am not doing in this PR. I think that none of the changes here will impact anything with Matrix because the relevant methods are not used by Matrix (also the new class introduced here is only use if SYMPY_GROUND_TYPES=flint).

I have separate PRs to make Matrix use DomainMatrix e.g. gh-25458 (charpoly), gh-25452 (inv) but the only one that is merged is gh-25443 (RREF) and RREF is one of the operations that I am not changing here. I could easily open more PRs like that but I am struggling to get any significant review on the existing open PRs so I won't push further on that right now.

This PR should not adversely affect anything that is directly used by users. It does however add the capability to SymPy to leverage state of the art implementations for many key matrix operations. More work would be needed to translate that into tangible benefits for end users (e.g. like gh-25458 and gh-25452) but this is an important step. Also there is plenty more in python-flint and in Flint itself that could be leveraged so this should be seen as just a first step in making use of the features there. The biggest benefit to SymPy from Flint would come from Flint's sparse polynomials which are not yet exposed in python-flint but someone else is working on that as I write this and I expect that the next release of python-flint will have something that SymPy can use.

[oscarbenjamin]

Thanks @smichr!

There was some discussion elsewhere about the requirements file but I'll follow that up in gh-25468