Sqrt simplification

[Eric Gourgoulhon]

Hi,

Here is some behavior of Sage which can be quite disturbing for a new user and probably would be considered as a bug:

sage: assume(x<0)
sage: sqrt(x^2).simplify_full()
x

It appears that the instruction "assume(x<0)" is not sufficient to enforce the desired behavior: x is still considered as a complex number and sqrt(x^2) returns *a* square root of x^2, not necessarily the positive one. To enforce x real, one should set:

sage: maxima_calculus.eval('domain:real')
'real'

Then the output of simplify_full() is correct:

sage: sqrt(x^2).simplify_full()         
-x

But the following problem remains:

sage: sqrt((x-1)^2).simplify_full()  # correct result:
-x + 1
sage: sqrt(x^2-2*x+1).simplify_full() # wrong result (and not consistent with the above):
x - 1

I know that the problem is due to Maxima function randcan and has been much discussed, in particular here:
http://www.math.utexas.edu/pipermail/maxima/2011/026097.html
What is the latest status of this ? In particular, is there any possibility to enforce the correct behavior by passing some options to Maxima ? (I've noticed that maxima_calculus.eval('radexpand:true') does not help).
Thank you for your help.

Eric.

[Michael Orlitzky]

On 06/15/2013 11:31 AM, Eric Gourgoulhon wrote:
> Hi,
>
> Here is some behavior of Sage which can be quite disturbing for a new
> user and probably would be considered as a bug:
>
> sage: assume(x<0)
> sage: sqrt(x^2).simplify_full()
> x
>
>

> I know that the problem is due to Maxima function randcan and has been
> much discussed, in particular here:
> http://www.math.utexas.edu/pipermail/maxima/2011/026097.html
> What is the latest status of this ? In particular, is there any
> possibility to enforce the correct behavior by passing some options to
> Maxima ? (I've noticed that maxima_calculus.eval('radexpand:true') does
> not help).

We must remove simplify_radical() from simplify_full() and rename it. It
doesn't do simplification, and as long as it has "simplify" in the name,
users are going to (correctly) conclude that this is a bug.

There are plenty of tickets about this. Most of the info can be found here:

http://trac.sagemath.org/sage_trac/ticket/12737

[kcrisman]

What is the latest status of this ? In particular, is there any possibility to enforce the correct behavior by passing some options to Maxima ? (I've noticed that maxima_calculus.eval('radexpand:true') does not help).

Yes, I think that is broken.  Or at any rate the documentation is wrong.  Maybe this would be a way to preserve some of radcan for simplify_full, and then simplify_radical could use radexpand. 

[rjf]

Since your belief about what the correct behavior is, is at best
arguable, and at worst, wrong, I would not expect it to
be "corrected"

Sqrt(anything)   has TWO values.
Choosing one should depend on the choice of square root branch, e.g.
assume(sqrt(E))

assuming something about E doesn't say much.

e.g. even sqrt(9)  is +-3.
who cares about the "sign" of 9.

In some engineering and physical circumstances, a branch of the sqrt
can be intuited.  But the math issue is not resolved by "assume(x>0)"  etc

[Eric Gourgoulhon]

I do agree that a complex number has two square roots. Let me illustrate my point of view on a simple example, on which actually I faced the problem.
I am using Sage to do differential geometry and in particular to manipulate various charts on a manifold.
Let us consider the hyperbolic space H^2. The polar chart (R,phi) associated with the Poincaré disk model of H^2 is related to the chart (r,phi') associated with the hyperboloidal model by the following transformations: phi'=phi and
r = 2*R/(1-R^2)  <==> R = r / (1 + sqrt(1+r^2))
R in [0,1), r in [0, +infinity)

Now in Sage:

sage: R = var('R')
sage: assume(R>=0) ; assume(R<1)
sage: r = 2*R/(1-R^2)
sage: R == ( r / (1+sqrt(1+r^2)) ).simplify_full()
R == -1/R

which is false and is very annoying when manipulating fields defined on H^2, when passing from one chart to the other.

A solution could be to define two functions:
- sqrt(x): the current one, which returns one of the two square roots of x, with no control of which one
- sqrt_real(x), which assumes x is real and positive and always returns the positive root of x, so that sqrt_real(x^2) = abs(x) would be guaranteed.

Eric.

[rjf]

If you don't like the semantics for square root, or for radcan  you can certainly
define another function and try to make its simplification and evaluation
conform to your needs.   Your "simple example" explanation is not
sufficiently simple for me, but if you do this in Maxima:

 E :  r/ (1+sqrt(1+r^2)) , r=2*R/(1-R^2);

F: radcan(E)

you get a value for F of  -1/R,  which is true, as R-> infinity.
Since      RADCAN IGNORES ASSUMPTIONS and treats all variables
in the radicals as they increase beyond bounds to choose a sign,   this is the right answer.

If you want to evaluate the expression E,
(2*R)/((1-R^2)*(sqrt((4*R^2)/(1-R^2)^2+1)+1))

without such simplification at a point p between 0 and 1,
you will  find that E evaluates to p.  (actually between -1 and 1).

Also, you might find it potentially useful to use radcan as intended, and note that

radcan( substitute (1/w, R, E))  is  1/w   [[ meaning that near zero, E(p)=p ]]


RJF

[Eric Gourgoulhon]

Thanks for your explanation regarding how radcan works.

If radcan ignores assumptions, how can one explain the following behavior:

sage: assume(x<0)
sage: sqrt(x^2).simplify_radical()
x

sage: maxima_calculus.eval('domain:real')
'real'

sage: sqrt(x^2).simplify_radical()
-x

It seems that assuming the domain is 'real' had an effect on radcan.

[Michael Orlitzky]

The code for simplify() just sends an expression to Maxima and back:

return self._parent(self._maxima_())

If you set the domain to real, and send that expression to Maxima and
back, you get abs(x):

sage: _ = maxima_lib.eval('domain:real')
sage: sqrt(x^2).simplify()
abs(x)

Now if you assume x<0, you get -x.

In simplify_radical(), we necessarily have to send the expression
through Maxima in order to call radcan() on it. As soon as the
expression hits Maxima, the above simplification happens (before radcan
can do anything).

[Eric Gourgoulhon]

Thanks for your explanation. So it is clear that radcan ignore assumptions.
As already mentioned by kcrisman, radcan does not take into account the option variable radexpand, which, if set to 'true' (and domain to 'real'),  would do the job, as suggested in the documentation
http://maxima.sourceforge.net/docs/manual/en/maxima_9.html

Prof. Richard Fateman, since it seems that you are the primary author of radcan (am I correct ?), would it require a lot of work to make a new version of radcan (with a new name, in order not to alter the current one), so that it takes radexpand into account ?

Eric.

[rjf]

On Thursday, June 20, 2013 1:04:51 AM UTC-7, Eric Gourgoulhon wrote:

Thanks for your explanation. So it is clear that radcan ignore assumptions.
As already mentioned by kcrisman, radcan does not take into account the option variable radexpand, which, if set to 'true' (and domain to 'real'),  would do the job, as suggested in the documentation
http://maxima.sourceforge.net/docs/manual/en/maxima_9.html

First of all,  what radexpand:true  does, seems to me to be nonsense.
I didn't write that and don't endorse its use.  As I've explained, sqrt(9) has 2 values, and
it doesn't depend on whether that 9 was created by squaring 3  or -3.
If you want to define a function that returns the positive square root of a positive real number, and call
that sqrt also, as the author of the radexpand related simplification seems to prefer,
that is possible, but not necessarily mathematically sensible. In some physical
circumstances it may be an appropriate model to use.  But to make the choice
of which branch to take for sqrt(E)  by some computation that says E is a perfect
square of something which in some putative form is a square of something which
is positive or negative, is hard to justify, as you have just seen.
sqrt(x*2)    vs  sqrt((-x)^2)   which requires the simplifier to say, in effect,
no no no  don't change (-x)^2 to x^2 until I get the whiff of its sign, because
that taint will affect the meaning of the sqrt....

Prof. Richard Fateman, since it seems that you are the primary author of radcan (am I correct ?), would it require a lot of work to make a new version of radcan (with a new name, in order not to alter the current one), so that it takes radexpand into account ?

Yes, I wrote radcan.  The full source of it is available. I think it is in file rat3e.lisp. Look for
(defun $radcan
or perhaps (defmfun $radcan

I think that if you wish to modify it to make a function with another name, you will first
have to specify what that function should do.  It can then be read in to the maxima
image.   (you can presumably prod maxima to do so by having it execute
"load("yourfilename.lisp") "     from Sage).

Radcan was written sometime before 1971.  Other people have thought about
what might be "better" for 40+  years.  If you come up with a truly 100%
satisfactory replacement, more power to you. ")
RJF


 

[Eric Gourgoulhon]

Le vendredi 21 juin 2013 07:00:54 UTC+2, rjf a écrit :

Yes, I wrote radcan.  The full source of it is available. I think it is in file rat3e.lisp. Look for
(defun $radcan
or perhaps (defmfun $radcan

I think that if you wish to modify it to make a function with another name, you will first
have to specify what that function should do.  It can then be read in to the maxima
image.   (you can presumably prod maxima to do so by having it execute
"load("yourfilename.lisp") "     from Sage).

Radcan was written sometime before 1971.  Other people have thought about
what might be "better" for 40+  years.  If you come up with a truly 100%
satisfactory replacement, more power to you. ")
RJF

Thanks for your answer.
I will have to remember my lisp courses then...
Meanwhile, I've written the following workaround in python: basically, it scans all the sqrt's in a given symbolic expression, send them to radcan, takes the absolute value of the output and then calls simplify(). It is not optimal but it works:

        sage: assume(x<0)     
        sage: simplify_sqrt_real( sqrt(x^2) )
        -x
        sage: simplify_sqrt_real( sqrt(x^2-2*x+1) )
        -x + 1
        sage: simplify_sqrt_real( sqrt(x^2) + sqrt(x^2-2*x+1) )
        -2*x + 1

Here is the source code:

def simplify_sqrt_real(expr):
    r"""
    Simplify sqrt in symbolic expressions in the real domain.
   
    EXAMPLES:
   
    Simplifications of basic expressions::
   
        sage: assume(x<0)     
        sage: simplify_sqrt_real( sqrt(x^2) )
        -x
        sage: simplify_sqrt_real( sqrt(x^2-2*x+1) )
        -x + 1
        sage: simplify_sqrt_real( sqrt(x^2) + sqrt(x^2-2*x+1) )
        -2*x + 1

    """
    from sage.symbolic.ring import SR
    from sage.calculus.calculus import maxima
    # 1/ Search for the sqrt's in expr
    sexpr = str(expr)
    if 'sqrt(' not in sexpr:  # no sqrt to simplify
        return expr
    pos_sqrts = []   # positions of the sqrt's in sexpr
    the_sqrts = []   # the sqrt sub-expressions in sexpr
    for pos in range(len(sexpr)):
        if sexpr[pos:pos+5] == 'sqrt(':
            pos_sqrts.append(pos)
            parenth = 1
            scan = pos+5
            while parenth != 0:
                if sexpr[scan] == '(': parenth += 1
                if sexpr[scan] == ')': parenth -= 1
                scan += 1
            the_sqrts.append( sexpr[pos:scan] )
    # 2/ Simplifications of the sqrt's
    new_expr = ""    # will contain the result
    pos0 = 0
    for i, pos in enumerate(pos_sqrts):
        # radcan is called on each sqrt:
        x = SR(the_sqrts[i])
        simpl = SR(x._maxima_().radcan())
        # the absolute value of radcan's output is taken, the call to simplify()
        # taking into account possible assumptions regarding the sign of simpl:
        new_expr += sexpr[pos0:pos] + '(' + str(abs(simpl).simplify()) + ')'
        pos0 = pos + len(the_sqrts[i])
    new_expr += sexpr[pos0:]
    return SR(new_expr)

[Burcin Erocal]

Instead of going through strings, you might want to try pattern
matching:

sage: t = sin(x)*sqrt(x^2)*e^(sqrt(x+1))
sage: w = SR.wild()
sage: t.find(w^(1/2))
[sqrt(x + 1), sqrt(x^2)]

Note that this won't catch sqrt() in the denominator. You will need to
specify a negative exponent for that:

sage: u = 1/sqrt(x+1)
sage: u.find(w^(-1/2))
[1/sqrt(x + 1)]

> # 2/ Simplifications of the sqrt's
> new_expr = "" # will contain the result
> pos0 = 0
> for i, pos in enumerate(pos_sqrts):
> # radcan is called on each sqrt:
> x = SR(the_sqrts[i])
> simpl = SR(x._maxima_().radcan())

This is x.maxima_methods().radcan().

> # the absolute value of radcan's output is taken, the call to
> simplify()
> # taking into account possible assumptions regarding the sign
> of simpl:
> new_expr += sexpr[pos0:pos] + '(' +
> str(abs(simpl).simplify()) + ')' pos0 = pos + len(the_sqrts[i])
> new_expr += sexpr[pos0:]
> return SR(new_expr)

Putting the new values back can be done with subs().

sage: t
sqrt(x^2)*e^(sqrt(x + 1))*sin(x)
sage: t.subs({sqrt(x^2): x})
x*e^(sqrt(x + 1))*sin(x)


Cheers,
Burcin

[rjf]

On Sunday, June 30, 2013 1:23:50 PM UTC-7, Eric Gourgoulhon wrote:

...
 

. It is not optimal but it works:

It works on this example in your opinion.  If you want to deal with square roots in this
case using some bogus relationship between the choice of branch and some alleged
internal "sign", I suppose that is your choice.

If it were mathematically meaningful to do so, it would generalize to say, cube roots,
don't you think?

What you've written is just a hack.  Not that maxima's general simplifier inspires confidence.

Maybe consult someone who has taken college courses in complex variables, conformal mapping,
and modern algebra.

      

[Eric Gourgoulhon]

Le mardi 2 juillet 2013 02:38:44 UTC+2, rjf a écrit :

What you've written is just a hack. 

 
Of course it's a hack; this is why I did not submit it as a patch for Sage.
As far as one restricts oneself to the REAL DOMAIN, I think it works. Please show me a counter-example.

Again, let me insist:
it is elementary mathematics that
sqrt: R+ --> R+, x |--> sqrt(x)
is a ONE-TO-ONE map, the reverse of which is
R+ --> R+, x |--> x^2
I think everybody agrees that when working only with real variables, this is the standard meaning of the sqrt function.
Then it must obey
sqrt(x^2) = abs(x)
and the above hack simply ensures this. This is why I called it 'simplify_sqrt_real' to insist that it is valid only for REAL-valued expressions of REAL variables.

[rjf]

If you called the function in question  "real positive branch of the square root of a positive number" rather than sqrt,
maybe you might get agreement.  Let's call that RPBSRPN.

  Your statement  then translates to RPBSRPN(x^2) = abs(x) .
But then if it ir R+-->R+,  the abs() is unnecessary,  and RPBSRPN(x^2) = x.

Surely you don't believe  that sqrt of positive numbers are always  positive.

(for the movie clip from Airplane -- Surely you can't be serious -- see
http://www.ign.com/top/movie-moments/88  )

 For example you can then derive the quadratic formula and show that a quadratic
equation has only one root.   (Or only has a single root if b^2-4ac > 0, your constraint.

As I suggested, you are free to come up with some other function, perhaps absqrt()
with your specified behavior,  but I would expect it to be of limited use.

[Eric Gourgoulhon]

Le mercredi 3 juillet 2013 01:07:35 UTC+2, rjf a écrit :

  Your statement  then translates to RPBSRPN(x^2) = abs(x) .
But then if it ir R+-->R+,  the abs() is unnecessary,  and RPBSRPN(x^2) = x.

No, the abs is necessary: consider the following function:
f : R --> R+, x |-->  RPBSRPN(x^2)
then, for x<0, f(x) = abs(x) = -x.

Surely you don't believe  that sqrt of positive numbers are always  positive.

Yes I do: although a positive number has two square roots, a positive one and a negative one, the standard expectation regarding the sqrt function is that sqrt(x) is THE positive square root of x if x is a positive number.
Hence my initial post in this thread.
 

[rjf]

On Wednesday, July 3, 2013 2:05:05 AM UTC-7, Eric Gourgoulhon wrote:

Le mercredi 3 juillet 2013 01:07:35 UTC+2, rjf a écrit :

  Your statement  then translates to RPBSRPN(x^2) = abs(x) .
But then if it ir R+-->R+,  the abs() is unnecessary,  and RPBSRPN(x^2) = x.

No, the abs is necessary: consider the following function:
f : R --> R+, x |-->  RPBSRPN(x^2)
then, for x<0, f(x) = abs(x) = -x.

If you change the function, sure.   If I change the function, say to
f:C --> R+
then sqrt( (-i)^2)   comes out as abs(i)  whatever that might be.  Probably 1.

 

Surely you don't believe  that sqrt of positive numbers are always  positive.

Yes I do:

Then perhaps I should not depend on your personal belief system and point out that
using such a simplistic approach would be a disaster for a computer system that
does algebra.

For fun, I looked at the wikipedia article on the topic.  It refers to
the principal square root of a positive number.  It is a relatively simple treatment;
using it unaltered as a basis for a computer system in which symbols (and rational powers other
than 1/2) exist would be problematical.

 

although a positive number has two square roots, a positive one and a negative one, the standard expectation regarding the sqrt function is that sqrt(x) is THE positive square root of x if x is a positive number.
Hence my initial post in this thread.

If you wish to name your function something other than sqrt, go for it.