updates to unrad, continuous_domain, roots_(cubic, quartic, quintic)

[smichr]

References to other Issues or PRs

fixes #19869
fixes #21263
fixes #21268
fixes #21287
closes #21049 and closes #21002 as an alternate
closes #21295

Brief description of what is fixed or changed

unrad was changed in #18324 and #21032 to return an expression even when the original expression didn't have persisting radicals; it should only return an expression when the expression is going to give a superset of the solutions to an expression, not when the expression can be simply expanded to remove radicals. Those changes have been reverted and the function modified to handle the case that prompted that change.

Other comments

Release Notes

• polys
  • roots_quintic better detects non-rational coefficients
  • roots_cubic results are simpler for cubics that can be written (or rewritten) as x**3 + y = 0
  • corner-case bugs fixed in to_rational_coeffs
• calculus
  • continuous_domain returns results for any power with even denominator, not just 2
• solvers
  • unrad better detects radicals to be removed even when they are hidden as bases of rational powers
  • unrad no longer makes the highest order term positive when making sign canonical
  • _has_rational_power (a private function) is no longer used and is removed
  • _solve_as_rational now fully expands the expression to be solved
  • the signature of _solve_radical was changed to accept the unrad form of f
  • _solveset will now use unrad to see if an equation needs to use _solve_radical before doing so
  • solveset now uses factor as a last-resort

[sympy-bot]

✅

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

• calculus

  • continuous_domain returns results for any power with even denominator, not just 2 (#21276 by @smichr)

• polys

  • roots_quintic better detects non-rational coefficients (#21276 by @smichr)

  • roots_cubic results are simpler for cubics that can be written (or rewritten) as x**3 + y = 0 (#21276 by @smichr)

  • corner-case bugs fixed in to_rational_coeffs (#21276 by @smichr)

• solvers

  • unrad better detects radicals to be removed even when they are hidden as bases of rational powers (#21276 by @smichr)

  • unrad no longer makes the highest order term positive when making sign canonical (#21276 by @smichr)

  • _has_rational_power (a private function) is no longer used and is removed (#21276 by @smichr)

  • _solve_as_rational now fully expands the expression to be solved (#21276 by @smichr)

  • the signature of _solve_radical was changed to accept the unrad form of f (#21276 by @smichr)

  • _solveset will now use unrad to see if an equation needs to use _solve_radical before doing so (#21276 by @smichr)

  • solveset now uses factor as a last-resort (#21276 by @smichr)

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

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. -->
fixes #19869 
fixes #21263 
fixes #21268
fixes #21287
closes #21049 and closes #21002 as an alternate
closes #21295
#### Brief description of what is fixed or changed

`unrad` was changed in #18324 and #21032 to return an expression even when the original expression didn't have persisting radicals; it should only return an expression when the expression is going to give a superset of the solutions to an expression, not when the expression can be simply expanded to remove radicals. Those changes have been reverted and the function modified to handle the case that prompted that change.

#### Other comments


#### Release Notes

<!-- BEGIN RELEASE NOTES -->
* polys
  * `roots_quintic` better detects non-rational coefficients
  * `roots_cubic` results are simpler for cubics that can be written (or rewritten) as `x**3 + y = 0`
  * corner-case bugs fixed in `to_rational_coeffs`
* calculus
  * `continuous_domain` returns results for any power with even denominator, not just 2
* solvers
  * `unrad` better detects radicals to be removed even when they are hidden as bases of rational powers
  * `unrad` no longer makes the highest order term positive when making sign canonical
  * `_has_rational_power` (a private function) is no longer used and is removed
  * `_solve_as_rational` now fully expands the expression to be solved
  * the signature of `_solve_radical` was changed to accept the unrad form of f
  * `_solveset` will now use `unrad` to see if an equation needs to use `_solve_radical` before doing so
  * `solveset` now uses `factor` as a last-resort
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[smichr]

The case from issue #21263 now gives the simpler result:

>>> roots_cubic(Poly(x**3+3*x**2 + 3*x + a*b*c + 1, x))
[(-a*b*c)**(1/3) - 1, (-a*b*c)**(1/3)*(-1/2 + sqrt(3)*I/2) - 1, (-a*b*c)**(1/3)*(-1/2 - sqrt(3)*I/2) - 1]
>>> [(x**3+3*x**2 + 3*x + a*b*c + 1).subs(x,i).simplify() for i in _]
[0, 0, 0]

And it can be demonstrated that the alternate forms calculated internally are correct:

>>> from sympy.polys.polyroots import roots_cubic
>>> var('y',positive=True)
>>> s=roots_cubic(Poly(x**3 +y, x))
>>> [(x**3 +y).subs(x,i).simplify() for i in s]
[0, 0, 0]
>>> s=roots_cubic(Poly(x**3 -y, x))
>>> [(x**3 -y).subs(x,i).simplify() for i in s]
[0, 0, 0]
>>>

[smichr]

@oscarbenjamin , unrad has not been optimized to reduce the number of symbols added but some of the behaviors that were causing problems have been eliminated.

[oscarbenjamin]

Why has the answer here become more complicated?

[smichr]

Because solveset bails out if it can't solve the expression as written:

>>> e
-_x*erf(sqrt(2)*(_x - mu)/(2*sigma)) - _x + mu*erf(sqrt(2)*(_x - mu)/(2*sigma))
+ mu + 2*x*(_x - mu)
>>> solveset(_, _x)
ConditionSet(_x, Eq(-_x*erf(sqrt(2)*(_x - mu)/(2*sigma)) - _x + mu*erf(sqrt(2)*(
_x - mu)/(2*sigma)) + mu + 2*x*(_x - mu), 0), Complexes)
>>> solveset(e.factor(), _x)
FiniteSet(mu, sqrt(2)*sigma*(sqrt(2)*mu/sigma + 2*erfinv(2*x - 1))/2)

[smichr]

This has been fixed. The equation represents an interesting case where signsimp helps a lot by doing a simple thing. And then factor tidies up the simplified expression and gives a product of two factors which can be easily solved:

>>> fprint(e)

-w*exp(-w**2 + w*y - y**2)*exp(w**2 - w*y + y**2)*erf(w - y) +
w*exp(-w**2 + w*y - y**2)*exp(w**2 - w*y + y**2) + 2*x*(-w + y) +
y*exp(-w**2 + w*y - y**2)*exp(w**2 - w*y + y**2)*erf(w - y) -
y*exp(-w**2 + w*y - y**2)*exp(w**2 - w*y + y**2)
>>> signsimp(e)
-w*erf(w - y) + w - 2*x*(w - y) + y*erf(w - y) - y
>>> factor(_)
(-w + y)*(2*x + erf(w - y) - 1)

[oscarbenjamin]

Here p comes from p = -e**2/12 - g above. The reason this simplify is added is because if e is an Add like 1 + sqrt(2) then we need to expand the power to see whether it is zero. I think it is better to use _mexpand here since it is clear why that might be needed in some cases.

There might be some cases that get left unresolved where simplify could handle them but in those cases we return a Piecewise which is still correct. The caller can use simplify if they want. Ideally refine would be better at handling this sort of thing or there should be a simplify_piecewise function whose job is precisely to simplify the conditional parts of a Piecewise.

[oscarbenjamin]

I think this looks good although I had a couple of questions

[smichr]

I think this looks good although I had a couple of questions

thanks for taking the time to look this over; I'll make a few mods and look it over again and commit.

[smichr]

@oscarbenjamin , what do you think about using factor as the fallback in solveset like we do in solve?

[oscarbenjamin]

I'm not sure... When is it needed?

[smichr]

I'm not sure... When is it needed?

see the test added in the last commit; also, the different solve path used by the changes here apparently put the expression in a state that wasn't recognized without factoring; you will see that the original stats tests (except for the additional exclusion of mu from the solution) has been restored.

[smichr]

Another way to think about such equations is to think about them in terms of equations that can be easily solved when written in terms of "separation of generators'. Such equations will have solutions for generators that are each independent of the others:

>>> eq = 2*x*(y - z) - y*erf(y - z) - y + z*erf(y - z) + z
>>> Poly(2*x*(y - z) - y*erf(y - z) - y + z*erf(y - z) + z).gens
(x, y, z, erf(y - z))

there are two generators involving y

>>> solve(eq, erf(y-z))
[2*x - 1]

This solution does not depend on y so eq has a factor of (erf(y-z)-2x+1). To eliminate it, replace that factor with 1

>>> eq.subs(erf(y - z), (2*x - 1) + 1).expand()
-y + z

And now we can solve that easily for the second generator

>>> solve(_, y)
[z]

So we know how eq can be factored as a result but did not need to factor to find the roots:

>>> cancel(eq/((y-z)*(erf(y - z) - (2*x - 1)))).has(y)
False

An equation that can be solved like this will have these properties:

• there are multiple generators in the symbol of interest
• each generator will appear in front of factors such that when collecting on each generator all symbols in the factors will be in the part of the collected expression that is independent of the generator (or else the expression could have been trivially had a term removed or it cannot be factored because it will have extra terms):

>>> collect(eq,u)
u*(-y + z) + 2*x*y - 2*x*z - y + z  <-- y and z are in the non-u part of the expression
>>> collect(eq,y)
u*z - 2*x*z + y*(-u + 2*x - 1) + z  <-- u and x are in the non-y pare of the expression

• when solving for its generators it will have solutions independent of the other generators
• the generators can be inverted to find the base symbol of interest

Interesting to note that if you replace erf(y - z) with sqrt(y) then you end up with an expression that SymPy does not recognize as being factorable...

>>> uu = eq.subs(erf(y - z), sqrt(y))
2*x*y - 2*x*z - y**(3/2) + sqrt(y)*z - y + z
>>> factor(uu) == uu
True

So currently, unrad would be used to obtain a cubic which would then be solved;

>>> unrad(uu)
(4*x**2*y**2 - 8*x**2*y*z + 4*x**2*z**2 - 4*x*y**2 + 8*x*y*z - 4*x*z**2 - y**3 +
 2*y**2*z + y**2 - y*z**2 - 2*y*z + z**2, [])
>>> [factor(i) for i in solve(uu,y)]
[z, (2*x - 1)**2]

To do the separation we would have to carefully separate the generators symbolically (e.g. getting sqrt(y)*y as u*y)

>>> sep = uu.xreplace({y**(S(3)/2):u*y}).xreplace({sqrt(y):u}); sep
-u*y + u*z + 2*x*y - 2*x*z - y + z

Then solve for each as

>>> solve(sep,u)
[2*x - 1]
>>> solve(sqrt(y)-_[0],y)  # invert the solution for u to give the solution for y
[(2*x - 1)**2]

and

>>> solve(sep,y)
[z]

I suspect this gets a little trickier to identify as factors increase? But maybe not -- if we had 3 generators, gi and the expression could be written as (g1 - a)*(g2 - b)*(g3-c) then we would have to have all of the permutations of terms involving g in the unfactored expression. Let's see:

>>> [i for i in Poly((sqrt(x)-a)*(root(x,3)-b)*(x-c)).gens if i.has(x)]
[x, sqrt(x), x**(1/3), x**(1/6)]
>>> [i for i in Poly((sqrt(x)-a)*(sqrt(x)-b)*(x-c)).gens if i.has(x)]
[x, sqrt(x)]

I'm not sure how you would identify that the last expression can be factored into 3 terms instead of 2. Looking at the highest power and thinking how it can be written as the sum of the generator powers? This is not a trivial exercise for me; it may be so for others (@jksuom ?). So my preference is to leave this at the present use of unrad and the new use of factor in solveset for now.

[oscarbenjamin]

Poly doesn't handle square roots of symbols very well:

In [53]: Poly(sqrt(x) + x)
Out[53]: Poly(x + (sqrt(x)), x, sqrt(x), domain='ZZ')

In [54]: Poly(sqrt(x) + x, sqrt(x))
Out[54]: Poly((sqrt(x)) + x, sqrt(x), domain='ZZ[x]')

There are two possible approaches to handling that in Poly:

1. Identify sqrt(y) as a generator and recognise y as sqrt(y)**2. This is a relatively simple approach that reduces such expressions to multivariate polynomials that have well used routines. This works even if the coefficients are EX domain (in so far as anything does).
2. Construct an extension field so the domain is QQ[y,t]/(t**2-y). This can be done but there is less code for implementing such a domain and it doesn't necessarily combine well with other more complicated expressions.

Your last example can be factored by manually implementing 1.:

In [49]: a, b, c, x = symbols('a, b, c, x')

In [50]: eq = expand((sqrt(x)-a)*(sqrt(x)-b)*(x-c))

In [51]: eq
Out[51]: 
                             3/2               3/2          2
-a⋅b⋅c + a⋅b⋅x + a⋅c⋅√x - a⋅x    + b⋅c⋅√x - b⋅x    - c⋅x + x 

In [52]: eq.subs(sqrt(x), t).factor().subs(t, sqrt(x))
Out[52]: (-a + √x)⋅(-b + √x)⋅(-c + x)

The extension approach is needed for more complicated examples where it isn't possible to just substitute out one symbol of interest e.g. if you have sqrt(x-1) and sqrt(x-2) etc.

[smichr]

OK, I'm going to commit what is here. I will open an issue about Poly identifying sqrt(x) as a generator instead of both x and sqrt(x).