Add roots, complex_roots and real_roots

[GiacomoPope]

The flint_poly is a child class of flint_elem and is used to build other polynomial types on top of, however, flint_poly has a method roots() which converts the self to acb type, which I assume is not the intension.

Suggestion is to remove this method and include it to acb if needed.

python-flint/src/flint/pyflint.pyx

Lines 215 to 220 in 71bf110

	def roots(self, **kwargs):
	"""
	Isolates the complex roots of *self*. See :meth:`.acb_poly.roots`
	for details.
	"""
	return acb_poly(self).roots(**kwargs)

[oscarbenjamin]

Is this a problem because of refactoring (like a cyclic import)?

I guess acb is the only type that would always be able to represent the roots of any fmpz_poly or fmpq_poly although it probably shouldn't be used for nmod_poly.

Perhaps e.g. fmpz_poly should have some sort of roots method that only finds integer roots? It doesn't look like Flint has an immediate function for that.

Asking the user to convert to acb_poly explicitly does not seem unreasonable to me although it would be a breaking change.

[fredrik-johansson]

By default, roots ought to return roots in the coefficient ring, e.g. integer roots for fmpz_poly and rational roots for fmpq_poly.

This is supported explicitly by the generics system (though not all conversions have been implemented).

In flint_ctypes, the roots method for gr_poly currently looks like this:

    def roots(self, domain=None):
        """
        Computes the roots in the coefficient ring of this polynomial,
        returning a tuple (``roots``, ``multiplicities``).
        If the ring is not algebraically closed, the sum of multiplicities
        can be smaller than the degree of the polynomial.
        If ``domain`` is given, returns roots in that ring instead.

            >>> (ZZx([3,2]) * ZZx([15,1])**2 * ZZx([-10,1])).roots()
            ([10, -15], [1, 2])
            >>> ZZx([1]).roots()
            ([], [])

        We consider roots of the zero polynomial to be ill-defined:

            >>> ZZx([]).roots()
            Traceback (most recent call last):
              ...
            ValueError

        We construct an integer polynomial with rational, real algebraic
        and complex algebraic roots and extract its roots over
        different domains:

            >>> f = ZZx([-2,0,1]) * ZZx([1, 0, 1]) * ZZx([3, 2])**2
            >>> f.roots()   # integer roots (there are none)
            ([], [])
            >>> f.roots(domain=QQ)    # rational roots
            ([-3/2], [2])
            >>> f.roots(domain=AA)     # real algebraic roots
            ([Root a = 1.41421 of a^2-2, Root a = -1.41421 of a^2-2, -3/2], [1, 1, 2])
            >>> f.roots(domain=QQbar)     # complex algebraic roots
            ([Root a = 1.00000*I of a^2+1, Root a = -1.00000*I of a^2+1, Root a = 1.41421 of a^2-2, Root a = -1.41421 of a^2-2, -3/2], [1, 1, 1, 1, 2])
            >>> f.roots(domain=RR)      # real ball roots
            ([[-1.414213562373095 +/- 4.89e-17], [1.414213562373095 +/- 4.89e-17], -1.500000000000000], [1, 1, 2])
            >>> f.roots(domain=CC)      # complex ball roots
            ([[-1.414213562373095 +/- 4.89e-17], [1.414213562373095 +/- 4.89e-17], 1.000000000000000*I, -1.000000000000000*I, -1.500000000000000], [1, 1, 1, 1, 2])
            >>> f.roots(RF)     # real floating-point roots
            ([-1.414213562373095, 1.414213562373095, -1.500000000000000], [1, 1, 2])
            >>> f.roots(CF)     # complex floating-point roots
            ([-1.414213562373095, 1.414213562373095, 1.000000000000000*I, -1.000000000000000*I, -1.500000000000000], [1, 1, 1, 1, 2])

        Calcium examples/tests:

            >>> PolynomialRing(CC_ca)([2,11,20,12]).roots()
            ([-0.666667 {-2/3}, -0.500000 {-1/2}], [1, 2])
            >>> PolynomialRing(RR_ca)([1,-1,0,1]).roots()
            ([-1.32472 {a where a = -1.32472 [a^3-a+1=0]}], [1])
            >>> PolynomialRing(CC_ca)([1,-1,0,1]).roots()
            ([-1.32472 {a where a = -1.32472 [a^3-a+1=0]}, 0.662359 + 0.562280*I {a where a = 0.662359 + 0.562280*I [a^3-a+1=0]}, 0.662359 - 0.562280*I {a where a = 0.662359 - 0.562280*I [a^3-a+1=0]}], [1, 1, 1])

        """
        Rx = self.parent()
        R = Rx._coefficient_ring
        mult = VecZZ()
        if domain is None:
            roots = Vec(R)()
            status = libgr.gr_poly_roots(roots._ref, mult._ref, self._ref, 0, R._ref)
        else:
            C = domain
            roots = Vec(C)()
            status = libgr.gr_poly_roots_other(roots._ref, mult._ref, self._ref, R._ref, 0, C._ref)
        if status:
            if status & GR_UNABLE: raise NotImplementedError
            if status & GR_DOMAIN: raise ValueError
        return (roots, mult)

[GiacomoPope]

The generics system in FLINT3 is really beautiful. Potentially the reworking into submodules in #61 ties in really nicely with the generics from FLINT3, as we can have these as parent classes of specific types.

For now what I've marked ready should be "fine" but it's really only baby steps to something which hopefully is more modular and easier to develop.

[oscarbenjamin]

roots ought to return roots in the coefficient ring, e.g. integer roots for fmpz_poly and rational roots for fmpq_poly

Does Flint already have functions for computing these?

I didn't see anything in the docs for fmpz_poly and fmpq_poly when I looked earlier.

Is there a more efficient method than just calling factor?

I would have thought that linear factors could be found more efficiently somehow than by computing a full factorisation.

[fredrik-johansson]

Does Flint already have functions for computing these?

Within the generics system gr_poly_roots_other offers such an interface. For example to support computing the arb roots of an fmpz poly, the arb domain provides this code: https://github.com/flintlib/flint2/blob/dcd25c0eb86defc67b5c42b270da8a9b68aa004b/src/gr/arb.c#L1634

Such code ought to be factored out to dedicated functions at least for the most common cases (e.g. fmpz poly -> arb roots); this just hasn't been done yet.

I would have thought that linear factors could be found more efficiently somehow than by computing a full factorisation.

Yes, this is definitely possible. There is such root-finding code for fmpz_mod polynomials, but for many other cases (e.g. fmpz) an optimized algorithm hasn't been implemented yet.

[oscarbenjamin]

If it is causing refactoring problems then I think that it is fine to just remove the fmpz_poly etc roots methods and leave roots as a method on acb_poly. If longer term we don't want to have fmpz_poly.roots return acb then at some point there will be a compatibility break so we might as well just remove it now.

I suggest adding a CHANGELOG section in the README file with a section for version 0.5.0 and saying that the roots method is removed but acb_poly(p).roots() can be used instead.

[GiacomoPope]

Yeah. We can even have a depreciated function for flint_poly.roots() which points to the right way to use acb.roots()

[GiacomoPope]

Just to add to this, the method of roots() currently on fmpz_poly doesn't convert then call roots, but instead computes complex roots as its own method. Which feels like the wrong thing to do for me, but obviously someone else disagrees

python-flint/src/flint/fmpz_poly.pyx

Lines 323 to 368 in da48588

	def 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([1]).roots()
	[]
	>>> fmpz_poly([2,0,1]).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():
	... print((c,m))
	...
	([-0.375000000000000 +/- 1.0e-19] + [0.599478940414090 +/- 5.75e-17]j, 1)
	([-0.375000000000000 +/- 1.0e-19] + [-0.599478940414090 +/- 5.75e-17]j, 1)
	([-0.735284727404843 +/- 4.11e-16], 3)
	([0.0494605455206031 +/- 1.33e-17] + [0.784693167647185 +/- 2.85e-16]j, 3)
	([0.0494605455206031 +/- 1.33e-17] + [-0.784693167647185 +/- 2.85e-16]j, 3)
	"""
	cdef fmpz_poly_factor_t fac
	cdef long deg, i, j
	cdef int exp, flags
	cdef acb_ptr croots
	if not self:
	return []
	flags = 0
	if verbose:
	flags = 1
	roots = []
	fmpz_poly_factor_init(fac)
	fmpz_poly_factor_squarefree(fac, self.val)
	for 0 <= i < fac.num:
	deg = fmpz_poly_degree(&fac.p[i]);
	exp = fac.exp[i]
	croots = _acb_vec_init(deg)
	arb_fmpz_poly_complex_roots(croots, &fac.p[i], flags, getprec())
	for 0 <= j < deg:
	v = acb()
	acb_set(v.val, &croots[j])
	roots.append((v, exp))
	_acb_vec_clear(croots, deg)
	fmpz_poly_factor_clear(fac)
	return roots

[fredrik-johansson]

Yes, that code ought to be changed.

[GiacomoPope]

I was looking at the docs and found that, for example, fmpz_poly has methods to find complex and real roots but not integer roots?

https://flintlib.org/doc/fmpz_poly.html

Are integer roots a Flint3 thing? Not sure how to progress best with this issue

[oscarbenjamin]

fmpz_poly has methods to find complex and real roots but not integer roots?

A complete factorisation algorithm can always find roots in the ground domain of any polynomial. Both fmpz_poly and fmpq_poly can use .factor() and check for linear factors. As mentioned above faster algorithms are not yet implemented.

[GiacomoPope]

Haha I just returned to this issue to say I realised this. I'll work on a naive implementation which returns a list of the (negative of) the constant coefficient of each linear factor and then we can also expose real_roots() and complex_roots()

[oscarbenjamin]

Looks like we haveroots and complex_roots now but not real_roots.

Types were these methods do not make sense should raise something like DomainError rather than AttributeError:

In [3]: nmod_poly([1, 1], 10).complex_roots()
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[3], line 1
----> 1 nmod_poly([1, 1], 10).complex_roots()

File flint/flint_base/flint_base.pyx:262, in flint.flint_base.flint_base.flint_poly.complex_roots()

AttributeError: Complex roots are not supported for this polynomial

[GiacomoPope]

Maybe we can go through in a new PR and add failures for the real roots and make sure the error raised is DomainError then?

[oscarbenjamin]

Yes, that sounds right.

[oscarbenjamin]

I left some comments about real_roots: #199 (comment).

There are two things to consider:

• Should fmpz_poly.real_roots() use arb_fmpz_poly_complex_roots or acb_poly_validate_real_roots (or something else)?
• Should there be a uniform interface for roots methods?

In [12]: acb_poly([1, 0, -1]).roots()
Out[12]: [-1.00000000000000, 1.00000000000000]

In [13]: fmpz_poly([1, 0, -1]).roots()
Out[13]: [(1, 1), (-1, 1)]

It isn't possible to handle multiplicity correctly for arb_poly or acb_poly because exact square-free factorisation is needed for that and if the poly is square-free then the multiplicity would always be 1. Maybe that just means that they have to be different...

[oscarbenjamin]

Another issue is about overlapping roots:

In [3]: r1 = 10**100

In [4]: r2 = 10**100 + 1

In [5]: x = fmpz_poly([0, 1])

In [6]: p = (x - r1)*(x - r2)**2

In [7]: p.complex_roots()
Out[7]: 
[([1.00000000000000e+100 +/- 3e+81], 1),
 ([1.00000000000000e+100 +/- 3e+81], 2)]

In [8]: [(rr1, _), (rr2, _)] = p.complex_roots()

In [9]: rr1.overlaps(rr2)
Out[9]: True

Ideally these roots would be refined until they no longer overlap but since they come from distinct square free factors that doesn't happen. If the multiplicities are equal then it should be possible to distinguish them:

In [10]: p = (x - r1)*(x - r2)

In [11]: [(rr1, _), (rr2, _)] = p.complex_roots()

In [12]: rr1.overlaps(rr2)
Out[12]: False

I guess that is why we want to use arb_fmpz_poly_complex_roots because it guarantees this isolation but it only does it one square-free factor at a time.

Maybe we need to take the square-free part (radical) of the polynomial and use that to compute the roots but then identify each root with a square-free factor of known multiplicity somehow (e.g. by evaluating each factor at each root, but then might need to increase precision...).

[oscarbenjamin]

This approach can be used to handle multiplicity of different roots while ensuring that they don't overlap. This is for fmpz_poly:

from flint import *

def complex_roots(p, prec=None):
    return _roots(p, prec, real=False)

def real_roots(p, prec=None):
    return _roots(p, prec, real=True)

def _roots(p, prec, real):
    """Real roots of an fmpz_poly."""
    if prec is None:
        prec = ctx.prec
    _, fac = p.factor_squarefree()
    rad = p / p.gcd(p.derivative())
    while True:
        roots = _get_roots(fac, rad, prec, real)
        if roots is not None:
            return roots
        prec *= 2

def _get_roots(fac, rad, prec, real):
    """Get real roots from squarefree factors and radical.

    If prec is not high enough returns None.
    """
    oldprec = ctx.prec
    ctx.prec = prec

    # Current implementation except should be simplified
    # to avoid factorisation.
    complex_roots = rad.complex_roots()

    # prec affects rad.complex_roots() and f(r)
    # May as well use a higher precision for f(r)
    ctx.prec *= 2

    roots = []
    for r, m_ in complex_roots:
        assert m_ == 1
        if real and not r.imag.overlaps(0):
            break
        elif real:
            r = r.real
        roots.append(r)

    roots_mult = []

    for r in roots:

        factor_already_found = False

        for f, m in fac:

            if f(r).overlaps(0):
                if factor_already_found:
                    # Could be a root of two factors...
                    return None
                else:
                    factor_already_found = True
                    roots_mult.append((r, m))

    ctx.prec = oldprec
    return roots_mult

I don't know if there is some better way than this.

Of course there should be a shortcut for the case where the original polynomial has only one square-free factor.