✅

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
  • Implemented power in CRootOf. (#27154 by @haru-44)

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

I have implemented the power function for `CRootOf`. For example, `r = CRootOf(x**3 + x + 1, 0)` is one of the roots of the equation $x^3 + x + 1 = 0$. Since $x^3 = -x - 1$, `r**3` can be calculated as `-r - 1`. Similarly, the powers of `CRootOf` can be transformed into a polynomial of degree less than the original polynomial.

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


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


#### Other comments


#### 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
  * Implemented power in `CRootOf`.
<!-- END RELEASE NOTES -->

See also gh-26743, gh-26740, gh-8552, gh-21343, gh-15753.

I think there was previous discussion that maybe this simplification should not happen automatically. When you look at something like gh-26816 it seems obvious that the simplification is desired. The case whee it may not be desired is just where a simple power expands to a larger polynomial.

On balance I think that this should simplify automatically just like integer powers of rational powers do.

@oscarbenjamin @sylee957
Thank you for reviewing.
I have made comprehensive revisions, including handling negative numbers, based on your suggestions.

We have set a condition that prevents expansion when abs(other) < self.poly.degree(), so depending on the degree, CRootOf**(-1) will not be expanded. I believe it is better for this not to expand automatically, but there may be cases where we intentionally want to expand it. Therefore, how about implementing CRootOf.invert() for that purpose?

    def invert(self, g=None, *gens, **args):
        if g is not None:
            raise NotImplementedError()
        try:
            inv = invert(self.poly.gen, self.poly)
        except NotInvertible:
            raise NotImplementedError()
        return inv.subs(self.poly.gen, self)

how about implementing CRootOf.invert()

That seems reasonable to me.

Also it would be good to have something that can rationalise the denominator to bring all RootOf into the numerator.

so depending on the degree, CRootOf**(-1) will not be expanded. I believe it is better for this not to expand automatically

I'm not sure about if it better not to expand that. I think that the only concern is that the number of terms can increase, for example, CRootOf(x**5 - x + 1, x, 0)**-1 would have 2 terms 1 - CRootOf**4, which was previously one fraction 1/CRootOf.

However, that is none issue because that should not be different than expanding to any positive power, let's say, CRootOf**6 and the number of term can increase CRootOf**2 - CRootOf which was previously one power term.

I'm not sure about if it better not to expand that. I think that the only concern is that the number of terms can increase, for example, CRootOf(x**5 - x + 1, x, 0)**-1 would have 2 terms 1 - CRootOf**4, which was previously one fraction 1/CRootOf.

I don't have any strong evidence either.
I was just concerned that users might be surprised if 1/x is automatically rewritten. However, it also seems unreasonable to expand only within a specific range of negative numbers.

I was just concerned that users might be surprised if 1/x is automatically rewritten.

Other kinds of division such as 1/(r + r**2) won't reduce to a polynomial function of r without something more to rationalise the denominator. If we have 1/r**n though with n >= deg we could either have that expand in the denominator like 1/p(r) or into the numerator like (r**n).invert(q(r)).

It would be worth checking how other systems handle their equivalents of RootOf.

The inverse of r + r**2 can be found as gcdex(r + r**2, r**3 + r + 1)[0] when r is CRootOf(x**3 + x + 1, 0). Similarly, I think the inverse of a polynomial in CRootOf can be converted into a polynomial. However, I am not sure if this should be handled in _eval_power.

The inverse of r + r**2 can be found as gcdex(r + r**2, r**3 + r + 1)[0] when r is CRootOf(x**3 + x + 1, 0)

Isn't that just what invert does?

However, I am not sure if this should be handled in _eval_power.

No, I don't think it should be handled in _eval_power. If that should be handled somewhere else (radsimp?) then the question is: what should be handled in _eval_power?

Here is how it works with surds:

In [17]: cbrt(2)
Out[17]: 
3 ___
╲╱ 2 

In [18]: cbrt(2)**2
Out[18]: 
 2/3
2   

In [19]: cbrt(2)**3
Out[19]: 2

In [20]: cbrt(2)**-1
Out[20]: 
 2/3
2   
────
 2  

In [21]: 1/(sqrt(2) + 1)
Out[21]: 
  1   
──────
1 + √2

In [22]: radsimp(1/(sqrt(2) + 1))
Out[22]: -1 + √2

RootOf of degree 2 expands to surds so we should perhaps want manipulation of RootOf and manipulation of surds to be similar to some extent.

If that should be handled somewhere else (radsimp?) then the question is: what should be handled in _eval_power?

How about organizing it so that _eval_power handles CRootOf**n and radsimp handles everything else? The expected role of _eval_power is to deal with powers of itself, and it probably isn’t intended to do anything beyond that.

I might not be following the conversation correctly. Apologies if I’m saying something irrelevant.

How about organizing it so that _eval_power handles CRootOf**n and radsimp handles everything else?

Yes, but there is the question about what to do when n is negative. Let's make this more concrete:

In [1]: p = x**5 - x + 1

In [2]: r = RootOf(p, 0)

In [3]: r
Out[3]: 
       ⎛ 5           ⎞
CRootOf⎝x  - x + 1, 0⎠

For 2 <= n < 5 we will have that r**n is just a Pow but for n >= 5 it will reduce e.g.:

r**5 -> r - 1
r**100 -> 627401*r**4 - 735723*r**3 + 864339*r**2 - 1006897*r + 540536

This is analogous to what happens with surds e.g.:

In [12]: s = 3**(S(1)/5)

In [13]: s
Out[13]: 
5 ___
╲╱ 3 

In [14]: s**5
Out[14]: 3

In [15]: s**100
Out[15]: 3486784401

For negative n there are two cases. For n < -5 we would have r**-n = 1/r**n where r**n would usually reduce. If we want to reduce r**5 in the numerator then we should also reduce it in the denominator in which case we can have e.g.:

r**-5 -> 1/r**5 -> 1/(r - 1)

We can instead though invert r**5 mod p to get

r**-5 -> -r**4 - r**3 - r**2 - r

This means that we do not introduce an Add in the denominator but we also don't allow r**n with abs(n) >= 5 to exist in either the numerator or the denominator.

The other negative n case is -5 < n < 0 in which case we could keep the negative power r**n or we could also transform that into an Add of positive powers in the numerator like:

1/r -> 1 - r**4

This is like what happens with surds except in the surd case there is no Add e.g.:

In [21]: s
Out[21]: 
5 ___
╲╱ 3 

In [22]: 1/s
Out[22]: 
 4/5
3   
────
 3 

The question now is of these should be chosen for automatic simplification. If we working in the algebraic field QQ[r] then everything would always canonicalise to a polynomial in positive powers:

In [26]: a = AlgebraicNumber(r, alias='alpha')

In [27]: K = QQ.algebraic_field(a)

In [28]: aa = K.from_sympy(a)

In [29]: K.to_sympy(1 + aa/(1 - aa))
Out[29]: 
 3        2    4
α  + α + α  + α 

I think that for Expr we don't expect this canonical form if there is an Add in the denominator unless some function is used to rationalise the denominator. We could make it so that all integer powers of r canonicalise to an Add of positive powers. What I am not sure about is if there is some situation in which keeping a power is better than introducing an Add somewhere.

I think we should investigate how other computer algebra systems like Mathematica, Maple, etc handle their versions of RootOf.

Thank you for your explanation. This is a difficult issue, and I believe it won't be resolved overnight (which is why there have been so many discussions on it so far).

I think I lack the ability to conduct further research or reach a conclusion here, so it might be best to close this PR.

Just want to point out that, when it comes time to reduce x**n mod p, SymPy does have an implementation of the linear recurrence relations for computing this reduction without actually doing polynomial division.

I don't know if it would be faster or not. Presumably the reason the linear recurrence approach is ever used is that it can be faster, but whether it actually is in this case might depend on implementation details.

In python-flint some poly types have a pow_mod method. It hasn't been added to fmpq_poly yet but e.g. for polynomials mod 11:

In [14]: x = flint.nmod_poly([0, 1], 11)

In [15]: p = x**5 + 2

In [16]: %time (1 + x).pow_mod(100000000000, p)
CPU times: user 98 μs, sys: 10 μs, total: 108 μs
Wall time: 282 μs
Out[16]: 4*x^4 + 6*x^3 + 8*x^2 + 4*x + 5

We don't use this yet in SymPy but it could be used in DUP_Flint for the domains where it is available. It could be added for fmpq_poly in python-flint as well but would have to be implemented there because FLINT does not provide a ready made function.

In the latest prerelease of python-flint (pip install -U --pre python-flint) it is possible to access number fields via the experimental generics system:

In [31]: import flint.types._gr as gr

In [32]: x = flint.fmpq_poly([0, 1])

In [33]: p = x**5 + 2

In [34]: K = gr.gr_nf_ctx.new(p)

In [35]: K
Out[35]: gr_nf_ctx(x^5 + 2)

In [36]: q = 1 + K.gen()

In [37]: q
Out[37]: a+1

In [38]: q**2
Out[38]: a^2 + 2*a + 1

In [39]: %time q**100
CPU times: user 79 μs, sys: 8 μs, total: 87 μs
Wall time: 109 μs
Out[39]: 2482843830166169701737607986025*a^4 + 2762435876506655586044308684300*a^3 + 1858222353424857283143655476550*a^2 - 191302639563165531348043321100*a - 2807500490603043882196052199999

In principle we can make computing powers of RootOf very fast. I suspect though that actually computing the power will likely not be the dominant cost. What matters more is the downstream impact on other symbolic manipulation routines that need to manipulate the expressions after reduction. Consider something like:

In [3]: rs = all_roots(Poly([1, 2, 3, 4, 5, 6], x))

In [4]: rs
Out[4]: 
⎡       ⎛ 5      4      3      2             ⎞         ⎛ 5      4      3      2             ⎞       ↪
⎣CRootOf⎝x  + 2⋅x  + 3⋅x  + 4⋅x  + 5⋅x + 6, 0⎠, CRootOf⎝x  + 2⋅x  + 3⋅x  + 4⋅x  + 5⋅x + 6, 1⎠, CRoo ↪

↪    ⎛ 5      4      3      2             ⎞         ⎛ 5      4      3      2             ⎞          ↪
↪ tOf⎝x  + 2⋅x  + 3⋅x  + 4⋅x  + 5⋅x + 6, 2⎠, CRootOf⎝x  + 2⋅x  + 3⋅x  + 4⋅x  + 5⋅x + 6, 3⎠, CRootOf ↪

↪ ⎛ 5      4      3      2             ⎞⎤
↪ ⎝x  + 2⋅x  + 3⋅x  + 4⋅x  + 5⋅x + 6, 4⎠⎦

In [5]: prod(rs)**5
Out[5]: 
                                             5                                              5       ↪
       ⎛ 5      4      3      2             ⎞         ⎛ 5      4      3      2             ⎞        ↪
CRootOf⎝x  + 2⋅x  + 3⋅x  + 4⋅x  + 5⋅x + 6, 0⎠ ⋅CRootOf⎝x  + 2⋅x  + 3⋅x  + 4⋅x  + 5⋅x + 6, 1⎠ ⋅CRoot ↪

↪                                         5                                              5          ↪
↪   ⎛ 5      4      3      2             ⎞         ⎛ 5      4      3      2             ⎞         ⎛ ↪
↪ Of⎝x  + 2⋅x  + 3⋅x  + 4⋅x  + 5⋅x + 6, 2⎠ ⋅CRootOf⎝x  + 2⋅x  + 3⋅x  + 4⋅x  + 5⋅x + 6, 3⎠ ⋅CRootOf⎝ ↪

↪                                      5
↪  5      4      3      2             ⎞ 
↪ x  + 2⋅x  + 3⋅x  + 4⋅x  + 5⋅x + 6, 4⎠ 

Now reduce each of these powers and the expression sort of explodes. Is that a likely problem? Perhaps not but the downstream manipulation of the expressions is likely a bigger performance concern than the reduction itself.

Also on performance python-flint has routines for isolating polynomial roots that are much faster than those in SymPy. Using those would make the construction of RootOf much faster in the first place.