real nth root function

[burcin]

See sage-devel threads at:

http://groups.google.com/group/sage-devel/t/cea9b562ea49c9c1

and

https://groups.google.com/forum/#!topic/sage-devel/Q8VLKBypcJk

CC: @orlitzky @kcrisman [email protected] @eviatarbach @slel

Component: symbolics

Author: Burcin Erocal, Kwankyu Lee

Branch: f8cb7a0

Reviewer: Karl-Dieter Crisman, Nils Bruin, Kwankyu Lee

Issue created by migration from https://trac.sagemath.org/ticket/12074

[burcin]

Work Issues: needs tests and documentation

[jdemeyer]

comment:2

Is there any chance this could be made into a function that we can do calculus with, like computing derivatives, integrals, solving equations... (I'm afraid the answer will be no though because we need maxima).

[burcin]

comment:3

Attachment: trac_12074-nth_root.patch.gz

Replying to @jdemeyer:

Is there any chance this could be made into a function that we can do calculus with, like computing derivatives, integrals, solving equations... (I'm afraid the answer will be no though because we need maxima).

This is already a symbolic function, so it plays well with symbolics generally (as opposed to piecewise functions for instance):

sage: v = nth_root(x,3)
sage: v*sin(x) + x^2
x^2 + real_nth_root(x, 3)*sin(x)

I updated the patch to add custom exponentiation and derivative methods as well:

sage: v^2
real_nth_root(x, 3/2)
sage: v*v
real_nth_root(x, 3/2)
sage: v.diff(x)
1/3*real_nth_root(x, -3/2)

This all needs a lot of work of course.

For integration and solving equations we call out to maxima. One way to get sensible results from these calls would be to convert this function to a regular (base)^(exp) representation when passing it to maxima. I don't think there is any way to read it back from the maxima result though.

[slel]

comment:10

See also sympy's real_root

• https://docs.sympy.org/latest/modules/functions/elementary.html#real-root

[kcrisman]

comment:11

There is some interest in this again, due to an upcoming book's 2nd edition.

Sympy's solution may be more elegant in some ways, because it uses Abs(x)**(1/3)*sign(x) for real_root(x,3), apparently. If that solution (in Sage, of course, with translations to Sympy) has better support overall in Sage we could just do that instead (for odd roots).

Here's the full Sympy code. They have better handling of complicated piecewise things, I guess.

    from sympy.functions.elementary.complexes import Abs, im, sign
    from sympy.functions.elementary.piecewise import Piecewise
    if n is not None:
        return Piecewise(
            (root(arg, n, evaluate=evaluate), Or(Eq(n, S.One), Eq(n, S.NegativeOne))),
            (Mul(sign(arg), root(Abs(arg), n, evaluate=evaluate), evaluate=evaluate),
            And(Eq(im(arg), S.Zero), Eq(Mod(n, 2), S.One))),
            (root(arg, n, evaluate=evaluate), True))
    rv = sympify(arg)
    n1pow = Transform(lambda x: -(-x.base)**x.exp,
                      lambda x:
                      x.is_Pow and
                      x.base.is_negative and
                      x.exp.is_Rational and
                      x.exp.p == 1 and x.exp.q % 2)
    return rv.xreplace(n1pow)

[kwankyu]

Changed author from Burcin Erocal to Burcin Erocal; Kwankyu Lee

[kwankyu]

Branch: public/12074

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

5028e22

Add real_nth_root symbolic function

[sagetrac-git]

Commit: 5028e22

[kwankyu]

comment:14

Revived Burcin's implementation of real_nth_root symbolic function.

It is in public. Feel free to improve.

[kwankyu]

Changed work issues from needs tests and documentation to none

[kwankyu]

Changed keywords from nth_root to none

[kwankyu]

Changed author from Burcin Erocal; Kwankyu Lee to Burcin Erocal, Kwankyu Lee

[kwankyu]

comment:36

Replying to @kcrisman:

by injecting some code into expression.pyx. The question is whether we want this.

I'm gonna say no on that - not without some healthy sage-devel discussion. I do think some people use the assume command for precisely this purpose, though.

Injecting the code incur many doctest failures in other places. So I give up this idea.

[kwankyu]

comment:37

Replying to @kcrisman:

• Do we want a latex method? E.g. latex_name=r"\mathrm{abs}" used elsewhere in this file. This is obviously not ideal:

sage: latex(F)
{\rm real}_{{\rm nth}_{{\rm root}}}\left(x, 5\right)

But I'm not sure whether just using the usual root syntax is "correct" either since we have a different function now. Using abs(x)^(1/3)*sgn(x) doesn't work for even n, else I'd recommend that for latex.

I see real_nth_root as representing the real function x^1/n with variable x taking real numbers (possibly negative for odd n). The latex print is for humans. So why not just print x^1/n?

• Integration, perhaps surprisingly, works, as you mention above:

sage: f = real_nth_root(x,3)
sage: f.diff()
1/3*real_nth_root(x^(-2), 3)
sage: f.integrate(x)
integrate(abs(x)^(1/3)*sgn(x), x)
sage: _.diff()
abs(x)^(1/3)*sgn(x)

Apparently integration now tries Maxima, Giac, and then Sympy, so that explains it. We should probably document this behavior.

Done.

• Upon further consideration, maybe the derivative can be a little better. Because aren't we sort of assuming this is real input? Maybe we aren't. Here is Sympy.

>>> from sympy import root, real_root, Rational
>>> from sympy.abc import x, n
>>> real_root(x,5)
Piecewise((Abs(x)**(1/5)*sign(x), Eq(im(x), 0)), (x**(1/5), True))
>>> diff(real_root(x,5),x)
Piecewise((Abs(x)**(1/5)*Derivative(sign(x), x) + (re(x)*Derivative(re(x), x) + im(x)*Derivative(im(x), x))*sign(x)**2/(5*x*Abs(x)**(4/5)), Eq(im(x), 0)), (1/(5*x**(4/5)), True))
>>> x = Symbol('x', real=True)
>>> diff(real_root(x,5),x)
2*Abs(x)**(1/5)*DiracDelta(x) + sign(x)**2/(5*Abs(x)**(4/5))

and here is Sage

sage: G = abs(x)^(1/5)*sgn(x)
sage: G
abs(x)^(1/5)*sgn(x)
sage: G.diff(x)
2*abs(x)^(1/5)*dirac_delta(x) + 1/10*(x + conjugate(x))*sgn(x)/abs(x)^(9/5)
sage: G.diff(x).simplify()
2*abs(x)^(1/5)*dirac_delta(x) + 1/5*x*sgn(x)/abs(x)^(9/5)

at least

sage: assume(x,'real')
sage: F.diff(x)
1/5*real_nth_root(x^(-4), 5)

should mimic this, if not perhaps always. Any thoughts on this?

I don't get your concern. I think there is nothing wrong with

sage: f = real_nth_root(x,3)
sage: f.diff()
1/3*real_nth_root(x^(-2), 3)

• To elaborate, Sympy's function somehow allows complex input, but do we want to?

sage: real_nth_root(3.*I,5)
1.18476052767182 + 0.384952030759866*I # should this happen?

That's of course related to Nils question regarding negative input too.

sage: H = real_nth_root(x,2)
sage: (H.subs(x=-1))^2
---------------------------------------------------------------------------
...
ValueError: no real nth root of negative real number with even n
sage: ((H)^2).subs(x=-1)
1

I'm not sure the latter should be allowed or not.

No. real_nth_root takes real input and outputs real. For the symbolic function 1/x, we assume x takes nonzero value. Likewise, we assume real value x for real_nth_root(x, 3).

[sagetrac-git]

Changed commit from d12ab31 to 00380b0

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

00380b0

Fixes for reviewer comments

[orlitzky]

comment:39

Replying to @kwankyu:

Injecting the code incur many doctest failures in other places. So I give up this idea.

I was thinking of something simpler, like

sage: nth_root(RR(-1), 3)
-1.00000000000000
sage: nth_root(CC(-1), 3)
0.500000000000000 + 0.866025403784439*I

[kcrisman]

comment:40

Thanks for all these fixes - we are converging :-) Couple more questions/comments:

---

• I don't think you needed to make this one not tested; it should have passed fine. Only the subsequent one would have thrown an error.

-    sage: plot(real_nth_root(x, 3), (x, -1, 1))
-    Graphics object consisting of 1 graphics primitive
+    sage: plot(real_nth_root(x, 3), (x, -1, 1))  # not tested

• latex:

I see real_nth_root as representing the real function x^1/n with variable x taking real numbers (possibly negative for odd n). The latex print is for humans. So why not just print x^1/n?

It's a little tricky, because the humans it is intended for might not realize it's not quite the same. mjo or Nils, any thoughts? I can go with it for now.

• calculus:

Done.

Thanks. However, that brings up the related issue that the online documentation is only going to show the _init_ method doc (and apparently not the tests?).

If you are up for it, I'd recommend making a more "narrative" structure which includes examples of calculus, plotting, etc., more like real_part does. If you'd prefer me to do that, I can - might just take a little longer as I recall how to add to a public branch :-)

• diff output:

sage: assume(x,'real')
sage: F.diff(x)
1/5*real_nth_root(x^(-4), 5)

should mimic this, if not perhaps always. Any thoughts on this?

I don't get your concern. I think there is nothing wrong with

sage: f = real_nth_root(x,3)
sage: f.diff()
1/3*real_nth_root(x^(-2), 3)

If n is odd, though, 2*abs(x)^(1/5)*dirac_delta(x) + 1/5*x*sgn(x)/abs(x)^(9/5) is somehow "better", and available. Just thinking out loud; could be for a followup ticket.

• input:

sage: H = real_nth_root(x,2)
sage: (H.subs(x=-1))^2
---------------------------------------------------------------------------
...
ValueError: no real nth root of negative real number with even n
sage: ((H)^2).subs(x=-1)
1

I'm not sure the latter should be allowed or not.

No. real_nth_root takes real input and outputs real. For the symbolic function 1/x, we assume x takes nonzero value. Likewise, we assume real value x for real_nth_root(x, 3).

Wait, so are you saying that we should disallow that function with x=-1 to be evaluated? Because 1/x blows up if we inject x=0 - though then again Sage simplifies 1/x*x =1 which isn't technically correct. Anyway, I think this is Nils' point.

[kcrisman]

comment:41

I was thinking of something simpler, like

sage: nth_root(RR(-1), 3)
-1.00000000000000
sage: nth_root(CC(-1), 3)
0.500000000000000 + 0.866025403784439*I

I feel like this would be another ticket, and need sage-devel discussion. This ticket is about adding something to aid in plotting, and then making sure it has enough bells and whistles that it doesn't end up causing confusion. I'm a little uncomfortable with this since

sage: RR(-1) == CC(-1)
True

and the point of the function is strictly for pedagogical purposes (which should also be made clear in the documentation).

[orlitzky]

comment:42

Replying to @kcrisman:

I feel like this would be another ticket, and need sage-devel discussion. This ticket is about adding something to aid in plotting, and then making sure it has enough bells and whistles that it doesn't end up causing confusion. I'm a little uncomfortable with this since

sage: RR(-1) == CC(-1)
True

and the point of the function is strictly for pedagogical purposes (which should also be made clear in the documentation).

We run into that same problem whenever we have a subclass method that does something different from its superclass method. I was only trying to clear up my previous suggestion. If the authors/reviewers don't think this is a good idea, feel free to ignore it.

[kcrisman]

comment:43

I was only trying to clear up my previous suggestion.

Fair enough.

If the authors/reviewers don't think this is a good idea, feel free to ignore it.

If this gets merged, why not open a new ticket for an nth_root function if you think it's more generally useful? SymPy has some explicit warnings about these things for its function but we could follow that model (its function is intended for more than the current ticket, I believe).

[orlitzky]

comment:44

Replying to @kcrisman:

If this gets merged, why not open a new ticket for an nth_root function if you think it's more generally useful? SymPy has some explicit warnings about these things for its function but we could follow that model (its function is intended for more than the current ticket, I believe).

I'm not 100% convinced it's the right thing to do either. Especially if we already have real_nth_root, then we have ^(1/n) for the complex root, and... presumably you know which one you want? I dunno. But until someone finds the existing methods lacking, adding them proactively just increases the chances that we'll have gotten it wrong when the need does arise.

[sagetrac-git]

Changed commit from 00380b0 to f8cb7a0

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

f8cb7a0

Refactored docs

[kwankyu]

comment:46

Replying to @kcrisman:

However, that brings up the related issue that the online documentation is only going to show the _init_ method doc (and apparently not the tests?).

If you are up for it, I'd recommend making a more "narrative" structure which includes examples of calculus, plotting, etc., more like real_part does. If you'd prefer me to do that, I can - might just take a little longer as I recall how to add to a public branch :-)

Done.

• input:

sage: H = real_nth_root(x,2)
sage: (H.subs(x=-1))^2
---------------------------------------------------------------------------
...
ValueError: no real nth root of negative real number with even n
sage: ((H)^2).subs(x=-1)
1

I'm not sure the latter should be allowed or not.

No. real_nth_root takes real input and outputs real. For the symbolic function 1/x, we assume x takes nonzero value. Likewise, we assume real value x for real_nth_root(x, 3).

Wait, so are you saying that we should disallow that function with x=-1 to be evaluated?

Yes. What should real_nth_root(-1, 2) evaluate to?

Sage simplifies 1/x*x =1 which isn't technically correct.

Why is it not correct? It is correct but you should not evaluate 1/x at x=0.

[kcrisman]

comment:47

Done.

Thanks, I'll check that out momentarily to confirm it builds properly. Looks good so far.

• input:

sage: H = real_nth_root(x,2)
sage: (H.subs(x=-1))^2
---------------------------------------------------------------------------
...
ValueError: no real nth root of negative real number with even n
sage: ((H)^2).subs(x=-1)
1

I'm not sure the latter should be allowed or not.

No. real_nth_root takes real input and outputs real. For the symbolic function 1/x, we assume x takes nonzero value. Likewise, we assume real value x for real_nth_root(x, 3).

Wait, so are you saying that we should disallow that function with x=-1 to be evaluated?

Yes. What should real_nth_root(-1, 2) evaluate to?

I meant ((real_nth_root(x, 2))^2).subs(x=-1). Are you saying this should be allowed or not allowed?

Sage simplifies 1/x*x =1 which isn't technically correct.

Why is it not correct? It is correct but you should not evaluate 1/x at x=0.

Depends on whether such expressions are meant to include domains - but anyway that is not directly relevant.

[kcrisman]

Reviewer: Karl-Dieter Crisman, Nils Bruin, Kwanyku Lee

[kcrisman]

comment:48

Thank you for all the good work - this is really long since overdue.

Yes. What should real_nth_root(-1, 2) evaluate to?

I meant ((real_nth_root(x, 2))^2).subs(x=-1). Are you saying this should be allowed or not allowed?

Again, see comment:21. It is good that real_nth_root(-1, 2) does not evaluate, but the question is whether the symbolic simplification should do this when ((real_nth_root(x, 2))^2) doesn't really make sense for x=-1. Here is another example.

sage: ((real_nth_root(x, 2))^2)
real_nth_root(x^2, 2)
sage: f(x) = _
sage: f
x |--> real_nth_root(x^2, 2)
sage: f(-1)
1

I think if in _power_ we had something like

if exp % 2 == 1:
    return self(base**power_param, exp)
else:
 ....

that might solve it (suitable tests included). What do you think? Unfortunately the most naive things I tried for the else clause led to nasty errors; I wasted about 20 minutes on abort errors, so I have a feeling that what I did leads to Pynac problems. (Anyone else interested in this particular question other than me or Kwankyu?)

Probably we don't need to assert that n is a positive integer since the documentation clearly says so and since why would anyone use it any other way? (Famous last words.)

[kwankyu]

comment:49

Replying to @kcrisman:

Thank you for all the good work - this is really long since overdue.

Yes. What should real_nth_root(-1, 2) evaluate to?

I meant ((real_nth_root(x, 2))^2).subs(x=-1). Are you saying this should be allowed or not allowed?

Allowed. This is exactly the same problem with

sage: (1/x)*x
1
sage: _.subs(x=0)
1

which we don't bother to fix in any way. If you argue the cases are different because the difference of the domains of definition is of measure zero, then how about this:

sage: f = 1/floor(abs(x))*floor(abs(x))
sage: f.subs(x=1/2)
1

Perhaps it is important that these things do not cause real problems because the domain of definition never gets smaller from the initial domain of definition that one starts with.

[kcrisman]

comment:50

Allowed. This is exactly the same problem with

Okay, I'll roll with that. Haven't heard from any of the other people regarding that here, I don't have a big commitment here since this isn't really switchable to "regular" root functions. Thank you!

[kwankyu]

comment:51

Thank you!

[egourgoulhon]

comment:52

Thank you very much for making this happen!

[vbraun]

Changed branch from public/12074 to f8cb7a0

[mkoeppe]

Changed reviewer from Karl-Dieter Crisman, Nils Bruin, Kwanyku Lee to Karl-Dieter Crisman, Nils Bruin, Kwankyu Lee

[mkoeppe]

Changed commit from f8cb7a0 to none