Add method to compute integer roots #72

GiacomoPope · 2023-09-07T11:13:51Z

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.

oscarbenjamin · 2023-09-07T11:29:28Z

src/flint/types/fmpz_poly.pyx

+        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)


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.

GiacomoPope · 2023-09-07T12:21:15Z

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

oscarbenjamin · 2023-09-07T12:36:58Z

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...

GiacomoPope · 2023-09-07T12:43:10Z

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.

oscarbenjamin · 2023-09-07T12:46:10Z

src/flint/flint_base/flint_base.pyx

+            >>> 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)]


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.

oscarbenjamin · 2023-09-07T12:54:21Z

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.

GiacomoPope · 2023-09-07T13:25:07Z

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

oscarbenjamin · 2023-09-07T13:33:42Z

Perhaps nmod_poly.complex_roots could give a better error message.

GiacomoPope · 2023-09-07T13:36:53Z

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?

oscarbenjamin · 2023-09-07T13:50:36Z

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.

GiacomoPope · 2023-09-07T14:41:18Z

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

GiacomoPope · 2023-09-07T14:46:14Z

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]

GiacomoPope · 2023-09-07T15:44:50Z

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.

oscarbenjamin · 2023-09-07T15:53:33Z

src/flint/flint_base/flint_base.pyx

+        return roots
+
+    def complex_roots(self):
+        raise AttributeError("Complex roots are not supported for this polynomial")


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.

oscarbenjamin · 2023-09-07T15:54:39Z

Okay, let's get this in.

oscarbenjamin · 2023-09-07T16:10:10Z

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.

GiacomoPope · 2023-09-07T16:17:24Z

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
16 checks passed

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

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add method to compute integer roots #72

Add method to compute integer roots #72

GiacomoPope commented Sep 7, 2023 •

edited

Loading

oscarbenjamin Sep 7, 2023

GiacomoPope Sep 7, 2023

oscarbenjamin Sep 7, 2023

GiacomoPope commented Sep 7, 2023 •

edited

Loading

oscarbenjamin commented Sep 7, 2023

GiacomoPope commented Sep 7, 2023

oscarbenjamin Sep 7, 2023

oscarbenjamin commented Sep 7, 2023

GiacomoPope commented Sep 7, 2023 •

edited

Loading

oscarbenjamin commented Sep 7, 2023

GiacomoPope commented Sep 7, 2023

oscarbenjamin commented Sep 7, 2023

GiacomoPope commented Sep 7, 2023

GiacomoPope commented Sep 7, 2023

GiacomoPope commented Sep 7, 2023

oscarbenjamin Sep 7, 2023

oscarbenjamin Sep 7, 2023

GiacomoPope Sep 7, 2023

oscarbenjamin commented Sep 7, 2023

oscarbenjamin commented Sep 7, 2023

GiacomoPope commented Sep 7, 2023

Add method to compute integer roots #72

Add method to compute integer roots #72

Conversation

GiacomoPope commented Sep 7, 2023 • edited Loading

oscarbenjamin Sep 7, 2023

Choose a reason for hiding this comment

GiacomoPope Sep 7, 2023

Choose a reason for hiding this comment

oscarbenjamin Sep 7, 2023

Choose a reason for hiding this comment

GiacomoPope commented Sep 7, 2023 • edited Loading

oscarbenjamin commented Sep 7, 2023

GiacomoPope commented Sep 7, 2023

oscarbenjamin Sep 7, 2023

Choose a reason for hiding this comment

oscarbenjamin commented Sep 7, 2023

GiacomoPope commented Sep 7, 2023 • edited Loading

oscarbenjamin commented Sep 7, 2023

GiacomoPope commented Sep 7, 2023

oscarbenjamin commented Sep 7, 2023

GiacomoPope commented Sep 7, 2023

GiacomoPope commented Sep 7, 2023

GiacomoPope commented Sep 7, 2023

oscarbenjamin Sep 7, 2023

Choose a reason for hiding this comment

oscarbenjamin Sep 7, 2023

Choose a reason for hiding this comment

GiacomoPope Sep 7, 2023

Choose a reason for hiding this comment

oscarbenjamin commented Sep 7, 2023

oscarbenjamin commented Sep 7, 2023

GiacomoPope commented Sep 7, 2023

GiacomoPope commented Sep 7, 2023 •

edited

Loading

GiacomoPope commented Sep 7, 2023 •

edited

Loading

GiacomoPope commented Sep 7, 2023 •

edited

Loading