Add method to compute integer roots by GiacomoPope · Pull Request #72 · flintlib/python-flint · GitHub

This is a work in progress. Currently both fmpz_poly and fmpq_poly have .roots() methods which return complex roots of the polynomial. The hope is to address issue #62.

This first set of commits renames these methods to complex_roots() and also introduces a new method to compute integer roots on fmpz_poly() which is named .roots()

I have updated and included a few more tests and added doctests to the new function.

If this is something the flint team like, I can add similar methods to fmpq_poly for rational roots and also look at the other polynomial classes.

It would be much easier to implement this as just Python code:

In [26]: x = flint.fmpz_poly([0, 1])

In [27]: p = x*(x - 1)*(x**2 + 2)

In [28]: [(-fac[0], mul) for fac, mul in p.factor()[1] if fac.degree() == 1]
Out[28]: [(0, 1), (1, 1)]

Then it could be a generic method defined in the superclass.

I kept this as cython in case it helped with performance but agree that we could instead write it in python this way.

If the code was in the base superclass is the best thing then to have import statements in the doctest for each flint type?

For the arb and acb we can simply have roots write over the superclass for the specific impl which currently exist?

Performance of anything else is unlikely to be a significant concern if the base routine is factor. Other algorithms could be much more efficient in common cases e.g. using the rational root theorem (and derivative for multiplicity):
https://en.wikipedia.org/wiki/Rational_root_theorem

Note also that Python code in a .pyx file will run faster than Python code in a .py file: Cython already translates the code to C to some extent even if it just looks like Python code.

If Flint had dedicated algorithms for this then it would be worth calling those from cython. Otherwise I don't think it is worth implementing anything nontrivial here and we should just try to keep the code here as simple as possible: it is better to spend time expanding the coverage of Flint's existing functionality than trying to optimise things that are not already optimised in Flint.

For the arb and acb we can simply have roots write over the superclass for the specific impl which currently exist?

Yes, and if any specialised algorithms are implemented in future for other cases like fmpz_poly they can override the generic method as well.

New commit adds a pythonic way to find roots.

To ensure this works with nmod_poly I first needed to allow a comparison between nmod and int types to check the linear term was one
roots now works nicely with nmod_poly, fmpz_poly and fmpq_poly
roots works differently on abc_poly (this doesn't show multiplicities?)
roots does not currently work on arb_poly as this class has no factor() method, or special method for computing roots.

>>> from flint import *
>>> fmpz_poly([12,7,1]).roots()
[(-3, 1), (-4, 1)]
>>> fmpq_poly([12,7,1], 163).roots()
[(-3, 1), (-4, 1)]
>>> nmod_poly([12,7,1], 11).roots()
[(8, 1), (7, 1)]
>>> acb_poly([12,7,1]).roots()
[[-4.0000000 +/- 1.96e-8] + [+/- 1.50e-8]j, [-3.00000000 +/- 7.46e-9] + [+/- 7.46e-9]j]
>>> arb_poly([12,7,1]).roots()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "src/flint/flint_base/flint_base.pyx", line 104, in flint.flint_base.flint_base.flint_poly.roots
    raise NotImplementedError("Polynomial has no factor method, roots cannot be determined")
NotImplementedError: Polynomial has no factor method, roots cannot be determined

roots does not currently work on arb_poly as this class has no factor() method, or special method for computing roots.

That is a shame because this would be useful. Actually for SymPy what would be even more useful is a function that can compute non-intersecting arb real roots of a (square free or irreducible) fmpq_poly. I wonder if that exists anywhere...

For real roots, we could take the complex roots and then simply filter those which are purely real? Not pretty, but it would at least "work".

I'm not familiar with what the state of the art, or best idea would be for this, but I'm happy to do the "stupid" idea following the above.

It would be better to add most of these examples just to the tests rather than the docstring. The primary purpose of the docstring should be for users to read as documentation rather than to function as a test suite. In that sense it is probably also good to mention how to find complex roots here because I expect some users would come to this function expecting it to do that.

For real roots, we could take the complex roots and then simply filter those which are purely real? Not pretty, but it would at least "work".

What I want is a function that provably separates the real and non-real roots and then bounds each real root into an interval [a, b] (open or closed?) on the real axis and each complex root into a complex square/disc that does not overlap the real axis. I need to be fully confident though that the real roots are definitely real and the complex roots are definitely not. Given an arb_poly it is impossible to do that because only exact arithmetic can definitely distinguish real and complex roots and the problem can be ill-posed at any level of approximate precision:
https://en.wikipedia.org/wiki/Wilkinson%27s_polynomial

However given an fmpq_poly it is possible and the results could be represented as arb. It is possible to compute only the real roots much more efficiently than also computing the complex roots and this is actually something that I would use a lot.

Current state of polynomial roots:

>>> from flint import *
>>> fmpz_poly([1,2]).roots()
[]
>>> fmpz_poly([2,1]).roots()
[(-2, 1)]
>>> fmpz_poly([1,2]).complex_roots()
[(-0.500000000000000, 1)]
>>> fmpz_poly([2,1]).complex_roots()
[(-2.00000000000000, 1)]
>>> 
>>> fmpq_poly([2,1]).roots()
[(-2, 1)]
>>> fmpq_poly([1,2]).roots()
[(-1/2, 1)]
>>> fmpq_poly([1,2]).complex_roots()
[(-0.500000000000000, 1)]
>>> fmpq_poly([2,1]).complex_roots()
[(-2.00000000000000, 1)]
>>> 
>>> nmod_poly([1,2],11).roots()
[(5, 1)]
>>> nmod_poly([1,2],11).complex_roots()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'flint.types.nmod_poly.nmod_poly' object has no attribute 'complex_roots'
>>> 
>>> acb_poly([1,2]).roots()
[-0.500000000000000]
>>> acb_poly([1,2]).complex_roots()
[-0.500000000000000]
>>> 
>>> arb_poly([1,2]).roots()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "src/flint/flint_base/flint_base.pyx", line 88, in flint.flint_base.flint_base.flint_poly.roots
    raise NotImplementedError("Polynomial has no factor method, roots cannot be determined")
NotImplementedError: Polynomial has no factor method, roots cannot be determined
>>> arb_poly([1,2]).complex_roots()
[-0.500000000000000]

I think the complex method should include multiplicities, but there's no support for squareful input atm?

>>> acb_poly([1,2,1]).roots()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "src/flint/types/acb_poly.pyx", line 402, in flint.types.acb_poly.acb_poly.roots
    raise ValueError("roots() failed to converge: insufficient precision, or squareful input")
ValueError: roots() failed to converge: insufficient precision, or squareful input

Perhaps nmod_poly.complex_roots could give a better error message.

For nmod_poly.complex_roots I had no real idea of the best practice. Here you'd want the closure of Z/NZ, maybe the best thing to do is add a NotImplementedError for complex_roots() on the superclass and include nothing within nmod_poly here?

I think the complex method should include multiplicities, but there's no support for squareful input atm?

It is not possible to distinguish close roots vs multiple roots if the coefficients are not exact so multiplicity is not well defined there.

The situation in practice would be that you

Start with some exact polynomials in e.g. $\mathbb{Q}(i,\sqrt{2})[x]$ or something.
Compute an exact factorisation into square free and coprime (e.g. irreducible) factors.
Then numerically evaluate the coefficients.
Finally ask acb_poly to find the roots of the square free factors.
Go back to step 3 if not enough precision was used.

maybe the best thing to do is add a NotImplementedError for complex_roots() on the superclass and include nothing within nmod_poly here?

There are different ways to arrange the code. Possibly if the other poly types all have this method then nmod_poly should have it too even if just to raise an error but in that case the error should not be AttributeError and the error message should make clear that it is intentional that nmod_poly does not implement this method.

It is not possible to distinguish close roots vs multiple roots if the coefficients are not exact so multiplicity is not well defined there.

I guess my issue is that we have a divergence of behaviour here:

>>> from flint import *
>>> fmpz_poly([1,2,1]).roots()
[(-1, 2)]
>>> fmpz_poly([1,2,1]).complex_roots()
[(-1.00000000000000, 2)]
>>> acb_poly([1,2,1]).complex_roots()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "src/flint/types/acb_poly.pyx", line 402, in flint.types.acb_poly.acb_poly.roots
    raise ValueError("roots() failed to converge: insufficient precision, or squareful input")
ValueError: roots() failed to converge: insufficient precision, or squareful input

Added an AttributeError to the base class, which means nmod_poly fails more gracefully:

>>> from flint import *
>>> 
>>> fmpz_poly([1,2,1]).roots()
[(-1, 2)]
>>> fmpz_poly([1,2,1]).complex_roots()
[(-1.00000000000000, 2)]
>>> fmpq_poly([1,2,1]).roots()
[(-1, 2)]
>>> fmpq_poly([1,2,1]).complex_roots()
[(-1.00000000000000, 2)]
>>> nmod_poly([1,2,1], 13).roots()
[(12, 2)]
>>> nmod_poly([1,2,1], 13).complex_roots()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "src/flint/flint_base/flint_base.pyx", line 100, in flint.flint_base.flint_base.flint_poly.complex_roots
    raise AttributeError("Complex roots are not supported for this polynomial")
AttributeError: Complex roots are not supported for this polynomial
>>> arb_poly([1,2,1]).roots()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "src/flint/flint_base/flint_base.pyx", line 89, in flint.flint_base.flint_base.flint_poly.roots
    raise NotImplementedError("Polynomial has no factor method, roots cannot be determined")
NotImplementedError: Polynomial has no factor method, roots cannot be determined
>>> arb_poly([1,2,1]).complex_roots()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "src/flint/types/arb_poly.pyx", line 121, in flint.types.arb_poly.arb_poly.complex_roots
    return acb_poly(self).roots(**kwargs)
  File "src/flint/types/acb_poly.pyx", line 402, in flint.types.acb_poly.acb_poly.roots
    raise ValueError("roots() failed to converge: insufficient precision, or squareful input")
ValueError: roots() failed to converge: insufficient precision, or squareful input
>>> arb_poly([1,2]).complex_roots()
[-0.500000000000000]
>>> acb_poly([1,2,1]).complex_roots()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "src/flint/types/acb_poly.pyx", line 402, in flint.types.acb_poly.acb_poly.roots
    raise ValueError("roots() failed to converge: insufficient precision, or squareful input")
ValueError: roots() failed to converge: insufficient precision, or squareful input
>>> acb_poly([1,2,1]).roots()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "src/flint/types/acb_poly.pyx", line 402, in flint.types.acb_poly.acb_poly.roots
    raise ValueError("roots() failed to converge: insufficient precision, or squareful input")
ValueError: roots() failed to converge: insufficient precision, or squareful input
>>> acb_poly([1,2]).roots()
[-0.500000000000000]

I'm marking this ready for review as it does what I thought needed to be done. There's a question about the consistency of the complex roots for various methods, but this is a different problem I suppose.

AttributeError in particular makes no sense here rather than e.g. NotImplementedError

Ideally it would be like DomainError or something at least for nmod_poly.

Perhaps not for this PR but exception handling needs a broad cleanup.

    def complex_roots(self):
        raise NotImplementedError("The method complex_roots is only defined for certain polynomial rings")

Might be better but I agree something like DomainError would be nice.

Okay, let's get this in.

It is not possible to distinguish close roots vs multiple roots if the coefficients are not exact so multiplicity is not well defined there.

I guess my issue is that we have a divergence of behaviour here:

>>> from flint import *
>>> fmpz_poly([1,2,1]).roots()
[(-1, 2)]
>>> fmpz_poly([1,2,1]).complex_roots()
[(-1.00000000000000, 2)]
>>> acb_poly([1,2,1]).complex_roots()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "src/flint/types/acb_poly.pyx", line 402, in flint.types.acb_poly.acb_poly.roots
    raise ValueError("roots() failed to converge: insufficient precision, or squareful input")
ValueError: roots() failed to converge: insufficient precision, or squareful input

I'm not sure what in general can be done about this. Essentially computing roots of polynomials with approximate coefficients is not something that can be completely well defined. I know that in context I would always use this the way that fmpz_poly.complex_roots does: compute the square free factorisation exactly first and then compute the approximate roots of the square free polynomials after.

I guess fundementally the difference is for fmpz_poly we start with full accuracy coefficients which allow us to make more assumptions than in acb_poly? Happy to do additional work, but I don't have a clear idea.

GiacomoPope added 3 commits September 7, 2023 12:10

update tests for complex and integer roots

1ddfb2f

rename method to show roots are complex

c4d7769

Rename complex roots and add integer roots method

288c7af

oscarbenjamin reviewed Sep 7, 2023

View reviewed changes

GiacomoPope added 3 commits September 7, 2023 13:18

allow comparison of nmod values to integers

221ac0b

Add generic method to compute roots from factorisation

29797bf

Simplify and remove specific method

9b6e697

oscarbenjamin reviewed Sep 7, 2023

View reviewed changes

GiacomoPope added 4 commits September 7, 2023 14:12

simplify doctests

ce720ab

simplify doctests

783ff05

include additional tests from unittest

eab8bdd

Add complex roots method

0a1c264

Add attribute error for base class for complex_roots

ce633dc

GiacomoPope marked this pull request as ready for review September 7, 2023 15:44

oscarbenjamin reviewed Sep 7, 2023

View reviewed changes

oscarbenjamin merged commit 4b7bab3 into flintlib:master Sep 7, 2023

GiacomoPope mentioned this pull request Sep 16, 2023

Bump version to 0.5.0 #84

Merged

oscarbenjamin mentioned this pull request Aug 26, 2024

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

Open

+                      fmpz_poly_factor(fac, self.val)
+                      for 0 <= i < fac.num:
+                          deg = fmpz_poly_degree(&fac.p[i])
+                          fmpz_poly_get_coeff_fmpz(linear_coeff.val, &fac.p[i], 1)

+                          >>> from flint import fmpz_poly
+                          >>> fmpz_poly([]).roots()
+                          []
+                          >>> fmpz_poly([1]).roots()
+                          []
+                          >>> fmpz_poly([2, 1]).roots()
+                          [(-2, 1)]
+                          >>> fmpz_poly([1, 2]).roots()
+                          []
+                          >>> fmpz_poly([12, 7, 1]).roots()
+                          [(-3, 1), (-4, 1)]
+                          >>> (fmpz_poly([1, 2]) *  fmpz_poly([-3,1]) *  fmpz_poly([1, 2, 3]) *  fmpz_poly([12, 7, 1])).roots()
+                          [(-3, 1), (-4, 1), (3, 1)]
+                          >>> from flint import nmod_poly
+                          >>> nmod_poly([1], 11).roots()
+                          []
+                          >>> nmod_poly([1, 2, 3], 11).roots()
+                          [(8, 1), (6, 1)]
+                          >>> nmod_poly([1, 6, 1, 8], 11).roots()
+                          [(5, 3)]
+                          >>> from flint import fmpq_poly
+                          >>> fmpq_poly([1,2]).roots()
+                          [(-1/2, 1)]
+                          >>> fmpq_poly([2,1]).roots()
+                          [(-2, 1)]
+                          >>> f = fmpq_poly([fmpq(1,3), fmpq(3,5)]) * fmpq_poly([fmpq(4,11), fmpq(9)])
+                          >>> f.roots()
+                          [(-4/99, 1), (-5/9, 1)]

Add method to compute integer roots by GiacomoPope · Pull Request #72 · flintlib/python-flint · GitHub

Repository navigation

Add method to compute integer roots - #72

Add method to compute integer roots#72oscarbenjamin merged 11 commits intoflintlib:masterflintlib/python-flint:masterfrom GiacomoPope:new_rootsGiacomoPope/python-flint:new_rootsCopy head branch name to clipboard

Conversation

GiacomoPope commented Sep 7, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

oscarbenjamin Sep 7, 2023

Choose a reason for hiding this comment

Uh oh!

GiacomoPope Sep 7, 2023

Choose a reason for hiding this comment

Uh oh!

oscarbenjamin Sep 7, 2023

Choose a reason for hiding this comment

Uh oh!

GiacomoPope commented Sep 7, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

oscarbenjamin commented Sep 7, 2023

Uh oh!

GiacomoPope commented Sep 7, 2023

Uh oh!

oscarbenjamin Sep 7, 2023

Choose a reason for hiding this comment

Uh oh!

oscarbenjamin commented Sep 7, 2023

Uh oh!

GiacomoPope commented Sep 7, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

oscarbenjamin commented Sep 7, 2023

Uh oh!

GiacomoPope commented Sep 7, 2023

Uh oh!

oscarbenjamin commented Sep 7, 2023

Uh oh!

GiacomoPope commented Sep 7, 2023

Uh oh!

GiacomoPope commented Sep 7, 2023

Uh oh!

GiacomoPope commented Sep 7, 2023

Uh oh!

oscarbenjamin Sep 7, 2023

Choose a reason for hiding this comment

Uh oh!

oscarbenjamin Sep 7, 2023

Choose a reason for hiding this comment

Uh oh!

GiacomoPope Sep 7, 2023

Choose a reason for hiding this comment

Uh oh!

oscarbenjamin commented Sep 7, 2023

Uh oh!

oscarbenjamin commented Sep 7, 2023

Uh oh!

GiacomoPope commented Sep 7, 2023

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Add method to compute integer roots#72
oscarbenjamin merged 11 commits into
flintlib:masterfrom
GiacomoPope:new_roots

GiacomoPope commented Sep 7, 2023 •
edited
Loading

GiacomoPope commented Sep 7, 2023 •
edited
Loading

GiacomoPope commented Sep 7, 2023 •
edited
Loading