Computing roots of polynomials with real algebraic coefficients

[mmaaz-git]

References to other Issues or PRs

Addresses a challenge brought up in the discussion in #26785. #23077 was also previously opened to do the same thing but was not merged. This was initially brought up in #22943 I believe.

Brief description of what is fixed or changed

We can now compute roots of polynomials with algebraic coefficients. This would previously fail. It computes the "rational lift" of the polynomial, i.e., a polynomial with rational coefficients whose roots are a superset of the original polynomial's, and then uses numerical comparison to find which roots are also roots of the original polynomial. This leverages the recent fixes in #26813 and #26816. The numerical comparison code is taken nearly verbatim from #22943 (comment).

Examples:

>>> from sympy import *
>>> from sympy.abc import x,y,z
>>> Poly(x**3 + sqrt(2)*x**2 - 1, x).real_roots()
[CRootOf(x**6 - 2*x**4 - 2*x**3 + 1, 0)]
>>> Poly(x**3 + sqrt(2)*x**2 - 1, x).all_roots()
[CRootOf(x**6 - 2*x**4 - 2*x**3 + 1, 0), CRootOf(x**6 - 2*x**4 - 2*x**3 + 1, 2), CRootOf(x**6 - 2*x**4 - 2*x**3 + 1, 3)]
>>> Poly(x**3 + CRootOf(y**5+y**3-1,0)*x**2 - 1, x).real_roots()
[CRootOf(x**15 + x**13 - 5*x**12 - 2*x**10 + 10*x**9 + 3*x**7 - 10*x**6 - x**4 + 5*x**3 - 1, 0)]

The fix is in ComplexRootOf._get_roots(), which is ultimately what Poly.real_roots() or Poly.all_roots() call (and hence ultimately what functions like solve() call I think). I think this makes logical sense in terms of integrating well into how the root functions already work. But, happy to hear feedback on this.

Other comments

This change breaks a test in test_quintics_3() in solvers/tests/test_solvers.py. The test asserted:

y = x**5 + x**3 - 2**Rational(1, 3)
assert solve(y) == solve(-y) == []

But now that this polynomial can be solved the roots should not be an empty list, but solve(y) == solve(-y) should still be correct. I changed the test accordingly. As well, the computation for complex roots can be somewhat slow, eg this test takes 210 seconds on my computer, hence I added a slow decorator.

Release Notes

• polys
  • It is now possible to compute the roots of polynomials with algebraic coefficients as RootOf using real_roots(p, extension=True) or all_roots(p, extension=True).

[sympy-bot]

✅

Hi, I am the SymPy bot. 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:

• polys
  • It is now possible to compute the roots of polynomials with algebraic coefficients as RootOf using real_roots(p, extension=True) or all_roots(p, extension=True). (#26945 by @mmaaz-git)

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

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

Addresses a challenge brought up in the discussion in #26785. #23077 was also previously opened to do the same thing but was not merged. This was initially brought up in #22943 I believe.

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

We can now compute roots of polynomials with algebraic coefficients. This would previously fail. It computes the "rational lift" of the polynomial, i.e., a polynomial with rational coefficients whose roots are a superset of the original polynomial's, and then uses numerical comparison to find which roots are also roots of the original polynomial. This leverages the recent fixes in #26813 and #26816. The numerical comparison code is taken nearly verbatim from https://github.com/sympy/sympy/issues/22943#issuecomment-2220562152.

Examples:
```
>>> from sympy import *
>>> from sympy.abc import x,y,z
>>> Poly(x**3 + sqrt(2)*x**2 - 1, x).real_roots()
[CRootOf(x**6 - 2*x**4 - 2*x**3 + 1, 0)]
>>> Poly(x**3 + sqrt(2)*x**2 - 1, x).all_roots()
[CRootOf(x**6 - 2*x**4 - 2*x**3 + 1, 0), CRootOf(x**6 - 2*x**4 - 2*x**3 + 1, 2), CRootOf(x**6 - 2*x**4 - 2*x**3 + 1, 3)]
>>> Poly(x**3 + CRootOf(y**5+y**3-1,0)*x**2 - 1, x).real_roots()
[CRootOf(x**15 + x**13 - 5*x**12 - 2*x**10 + 10*x**9 + 3*x**7 - 10*x**6 - x**4 + 5*x**3 - 1, 0)]
```

The fix is in ComplexRootOf._get_roots(), which is ultimately what Poly.real_roots() or Poly.all_roots() call (and hence ultimately what functions like solve() call I think). I think this makes logical sense in terms of integrating well into how the root functions already work. But, happy to hear feedback on this.

#### Other comments

This change breaks a test in `test_quintics_3()` in solvers/tests/test_solvers.py. The test asserted:

```
y = x**5 + x**3 - 2**Rational(1, 3)
assert solve(y) == solve(-y) == []
```

But now that this polynomial can be solved the roots should not be an empty list, but `solve(y) == solve(-y)` should still be correct. I changed the test accordingly. As well, the computation for complex roots can be somewhat slow, eg this test takes 210 seconds on my computer, hence I added a slow decorator.


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

* polys
  * It is now possible to compute the roots of polynomials with algebraic coefficients as RootOf using `real_roots(p, extension=True)` or `all_roots(p, extension=True)`.

<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[mmaaz-git]

Had to change the real_roots() and all_roots() doctests in polys/polytools.py as they previously expected a not implemented exception for algebraic coefficients

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

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

[oscarbenjamin]

There seem to be bugs somewhere because these cases should succeed EDIT: they do succeed on master:

In [12]: p = Poly(x**5 + x**3 - 2**Rational(1, 3), x, extension=True)

In [13]: p
Out[13]: Poly(x**5 + x**3 - 2**(1/3), x, domain='QQ<2**(1/3)>')

In [14]: p.lift()
---------------------------------------------------------------------------
CoercionFailed 
...
CoercionFailed: Cannot convert ANP([-1, 0, 0], [1, 0, 0, -2], QQ) of type <class 'sympy.polys.polyclasses.ANP'> from QQ<2**(1/3)> to QQ

Also

In [15]: Poly(x**2 + I).lift()
---------------------------------------------------------------------------
DomainError
...
DomainError: computation can be done only in an algebraic domain

[oscarbenjamin]

There seem to be bugs somewhere because these cases should succeed:

Forget that. I was testing the wrong branch.

[mmaaz-git]

I've changed the implementation. There are now two functions _get_roots_qq() and _get_roots_alg(). It performs a check and dispatches as suggested. As well, the extension is not constructed within CRootOf anymore, the responsibility is punted "higher up", so that it can be turned off. Lastly, I've modified the loop for _numerically_filter_roots so that it keeps track of if it's already processed the root, in the case of duplicates.

Since the extension has to be set before calling the root, I changed the real_roots() and all_roots() functions to set extension=True. This conflicted with the current greedy setting so I removed it for now. I'm not sure if this breaks anything else. Also, the unit tests now specify extension=True.

Basically, this won't work:

Poly(sqrt(2)*x**2 - 1, x).real_roots()

But this will:

Poly(sqrt(2)*x**2 - 1, x, extension=True).real_roots()

Not sure if you think this is good. I'm not sure what the best way is for users to interact with this.

[oscarbenjamin]

I don't think that this fallback should be used. It cannot reliably order the roots. If we get here then it is almost certainly because two or more roots have identical real part but calling re(root.eval()) will make them seem not equal leading to an arbitrary order.

Is it not possible to continue using the existing code that sorts the roots?

I realise everything here is that the code here is quite convoluted so making it difficult to do these things.

Here is one way to do this:

roots = p.lift().real_roots() # Existing QQ code can sort the roots

subroots = set()
for f, m in p.sqf_list()[1]:
    subroots.update(p.which_roots(roots, real=True))

roots_filtered = [r for r in roots if r not in subroots]

We only need to make this square-free for the real_roots case and only because otherwise count_roots does not give us the multiplicities. For the all_roots case it is just:

roots_filtered = p.which_roots(p.lift().all_roots(), real=False)

Given that functions like real_roots and all_roots are separate functions rather than being separated by a keyword argument I think it would make sense to have which_real_roots_sqf and which_roots as separate methods. The which_real_roots_sqf function requires the polynomial to be square free. The which_roots function does not. A function which_real_roots() could do the dance with sqf_list that I showed above before calling which_real_roots_sqf().

[mmaaz-git]

This is great! Thanks for the help.

After some experimenting, I think that even in the all_roots case, we need to get square-free factors. Eg if f= Poly((x-1) * (x-sqrt(2))**2 * (x-I) * (x+I), x, extension=True) (which has roots [1,sqrt(2),sqrt(2), I, -I]).

Bu it's lift has roots [-sqrt(2), -sqrt(2), 1, 1, sqrt(2), sqrt(2), -I, -I, I, I]. The -sqrt(2) are filtered out easily. However, in which_roots, if we just compare the sum of counts to the degree (as in the all_roots setting), this will loop forever. Thus, I think it makes sense to always return unique roots in which_roots(). And then handle multiplicity at the level of the square-free factorization.

In light of that, I think both cases are similar enough that we can keep which_roots one function, separated by the real argument? The only difference is the root counting. Otherwise, both should return unique elements.

[oscarbenjamin]

I think that even in the all_roots case, we need to get square-free factors.

Yes, you are right. The multiplicities might not be equal so we can't count the roots with multiplicities except in the square free case where the multiplicity is always 1.

[oscarbenjamin]

We can write this as just subroots[r] = m.

Is it possible for distinct square-free factors to have any roots in common?

I don't think it is. Let's suppose that fi and fj are factors with multiplicity mi and mj. If some number r is a root of both fi and fj then it must be a root of some fij that is a factor of both fi and fj and then the multiplicity of fij would be at least mi + mj.

To put it another way the distinct square-free factorisation is defined so that the factors are pairwise coprime (the algorithm literally divides out the gcds). That means that they cannot have any roots in common.

[mmaaz-git]

Ah great point, I put that in case, but from your explanation I can see it wouldn't happen.

[oscarbenjamin]

We can make this simpler:

roots_flat = []
for r in roots:
    if r in subroots:
        m = subroots[r]
        roots_flat.extend([r] * m)

[mmaaz-git]

Sounds good, although we have to add a way to not add the same root again to roots_flat. Eg with this code we get

In [5]: Poly((x-1) * (x-sqrt(2))**2, x, extension=True).real_roots()
Out[5]: [1, 1, sqrt(2), sqrt(2), sqrt(2), sqrt(2)]

Bc the lift has roots [-sqrt(2), -sqrt(2), 1, 1, sqrt(2), sqrt(2)], so the sqrt(2) gets looped over twice (and mult is 2, so we get 4 of them in the final output).

I've added a check for this to your code.

        roots_seen = set()
        roots_flat = []
        for r in roots:
            if r in subroots and r not in roots_seen:
                m = subroots[r]
                roots_flat.extend([r] * m)
                roots_seen.add(r)

[oscarbenjamin]

There should be a simpler way of writing this. One problem with the above is that we call Poly(f.expr, extension=True) if we have a Poly with domain QQ which is unnecessary.

It seems like something is wrong if we can't do all of this in one call. But I can't think what the call should be...

[mmaaz-git]

One reason I have all these checks is because if we just did something like:

if extension:
            F = Poly(f, extension=True)
else:
            F = Poly(f, greedy=False)

the issue is that if specify extension=True if the Poly is already made, it won't change it. E.g.

In [3]: f = x - sqrt(2)

In [4]: f
Out[4]: x - sqrt(2)

In [5]: Poly(f, x)
Out[5]: Poly(x - sqrt(2), x, domain='EX')

In [6]: Poly(f, x, extension=True)
Out[6]: Poly(x - sqrt(2), x, domain='QQ<sqrt(2)>')

In [7]: Poly(Poly(f, x), extension=True)
Out[7]: Poly(x - sqrt(2), x, domain='EX')

In [8]: Poly(Poly(f, x).expr, extension=True)
Out[8]: Poly(x - sqrt(2), x, domain='QQ<sqrt(2)>')

But as the last output shows, if you use the expr then you can construct the extension. It seems important to me to handle the case that extension=True and the user passes a Poly. I'm not sure if there is a better way to construct the domain in this case.

[oscarbenjamin]

Maybe if the input is a Poly it is reasonable to just ignore the other options like extension=True.

Either way I was more complaining that there isn't a simple way to do what is being done here that does not require handling 4 separate cases.

[mmaaz-git]

Hmm looks like the PR is failing bc of Ruff (mentioned in #26964)

[oscarbenjamin]

A lot of the code in which_real_roots and which_all_roots is duplicated. It would be better factor it out into a separate method that could be called by both public methods.

[oscarbenjamin]

Apart from this it all looks good though so let's fix this and get it merged.

I suggest having a private method like

def _which_roots(self, num_roots, candidates):
   ...

Then both which_all_roots and which_real_roots can call that method and it can house all the otherwise repeated code.

[mmaaz-git]

Done!

[mmaaz-git]

I changed the release notes entry
"Added ability to compute roots of polynomials with real algebraic coefficients"
to
"Added ability to compute roots of polynomials with algebraic coefficients"

since we also handle gaussian domains now

[oscarbenjamin]

Okay, looks good. Thanks!

[oscarbenjamin]

This is a nice addition and I think it was good that we went through the full review process on a smaller PR so that there is a clear idea of what it is that is needed to the final product into the codebase.

Looking forward to seeing the next steps for CAD!

[mmaaz-git]

Thanks! I learned a lot during the process :)