Skip to content

Add method to compute integer roots #72

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 11 commits into from
Sep 7, 2023
14 changes: 10 additions & 4 deletions src/flint/test/test.py
Original file line number Diff line number Diff line change
@@ -438,11 +438,17 @@ def test_fmpz_poly():
assert Z([1,2,2]).sqrt() is None
assert Z([1,0,2,0,3]).deflation() == (Z([1,2,3]), 2)
assert Z([1,1]).deflation() == (Z([1,1]), 1)
[(r,m)] = Z([1,1]).roots()
[(r,m)] = Z([1,1]).complex_roots()
assert m == 1
assert r.overlaps(-1)
assert Z([]).complex_roots() == []
assert Z([1]).complex_roots() == []
[(r,m)] = Z([1,1]).roots()
assert m == 1
assert r == -1
assert Z([]).roots() == []
assert Z([1]).roots() == []
assert Z([1, 2]).roots() == []

def test_fmpz_poly_factor():
Z = flint.fmpz_poly
@@ -985,9 +991,9 @@ def set_bad():
assert Q.bernoulli_poly(3) == Q([0,1,-3,2],2)
assert Q.euler_poly(3) == Q([1,0,-6,4],4)
assert Q.legendre_p(3) == Q([0,-3,0,5],2)
assert Q([]).roots() == []
assert Q([1]).roots() == []
[(r,m)] = Q([1,1]).roots()
assert Q([]).complex_roots() == []
assert Q([1]).complex_roots() == []
[(r,m)] = Q([1,1]).complex_roots()
assert m == 1
assert r.overlaps(-1)

6 changes: 3 additions & 3 deletions src/flint/types/fmpq_poly.pyx
Original file line number Diff line number Diff line change
@@ -384,16 +384,16 @@ cdef class fmpq_poly(flint_poly):
fac[i] = (base, exp)
return c / self.denom(), fac

def roots(self, **kwargs):
def complex_roots(self, **kwargs):
"""
Computes the complex roots of this polynomial. See
:meth:`.fmpz_poly.roots`.

>>> from flint import fmpq
>>> fmpq_poly([fmpq(2,3),1]).roots()
>>> fmpq_poly([fmpq(2,3),1]).complex_roots()
[([-0.666666666666667 +/- 3.34e-16], 1)]
"""
return self.numer().roots(**kwargs)
return self.numer().complex_roots(**kwargs)

@staticmethod
def bernoulli_poly(n):
52 changes: 48 additions & 4 deletions src/flint/types/fmpz_poly.pyx
Original file line number Diff line number Diff line change
@@ -344,18 +344,62 @@ cdef class fmpz_poly(flint_poly):
return c, res

def roots(self, bint verbose=False):
"""
Computes all the integer roots of this polynomial.
Returns a list of all pairs (*v*, *m*) where *v* is the
integer root and *m* is the multiplicity of the root

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

"""
cdef fmpz_poly_factor_t fac
cdef fmpz constant_coeff, linear_coeff
cdef long deg, i
cdef int exp

if not self:
return []

roots = []
constant_coeff = fmpz()
linear_coeff = fmpz()

fmpz_poly_factor_init(fac)
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)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

if deg == 1 and linear_coeff == 1:
fmpz_poly_get_coeff_fmpz(constant_coeff.val, &fac.p[i], 0)
exp = fac.exp[i]
roots.append((-constant_coeff, exp))
fmpz_poly_factor_clear(fac)
return roots

def complex_roots(self, bint verbose=False):
"""
Computes all the complex roots of this polynomial.
Returns a list of pairs (*c*, *m*) where *c* is the root
as an *acb* and *m* is the multiplicity of the root.

>>> fmpz_poly([]).roots()
>>> fmpz_poly([]).complex_roots()
[]
>>> fmpz_poly([1]).roots()
>>> fmpz_poly([1]).complex_roots()
[]
>>> fmpz_poly([2,0,1]).roots()
>>> fmpz_poly([2,0,1]).complex_roots()
[([1.41421356237310 +/- 4.96e-15]j, 1), ([-1.41421356237310 +/- 4.96e-15]j, 1)]
>>> for c, m in (fmpz_poly([2,3,4]) * fmpz_poly([5,6,7,11])**3).roots():
>>> for c, m in (fmpz_poly([2,3,4]) * fmpz_poly([5,6,7,11])**3).complex_roots():
... print((c,m))
...
([-0.375000000000000 +/- 1.0e-19] + [0.599478940414090 +/- 5.75e-17]j, 1)