#1515 Wrong factorization of sqrt()


Dear Developers of Maxima,

I found a wrong behavihor of sqrt().
It is a principle that sqrt(X*Y) should not be factorized to sqrt(X)*sqrt(Y)
unless X and Y can be judged to be non-negative.
I found a case in which this principle is broken.

Here is a program for demonstration:
* wrong_factorization_of_sqrt.maxima:
* S.Adachi 2008/10/01


/* real for 1 <= x <= 2 */

/* real for 2-sqrt(2) <= x <= 2+sqrt(2) */



/* END */

The result of execution is as follows:
Maxima 5.14.0cvs http://maxima.sourceforge.net
Using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (aka GCL)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) batch(wrong_factorization_of_sqrt.maxima)

batching #p/Volumes/HFS+2/home/adachi/work/307/wrong_factorization_of_sqrt.maxima
(%i2) display2d : false
(%o2) false
(%i3) correct:sqrt((1-1/x)*(2/x-1))
(%o3) sqrt((1-1/x)*(2/x-1))
(%i4) INCORRECT:sqrt((1-(2-sqrt(2))/x)*((sqrt(2)+2)/x-1))
(%o4) sqrt((2-sqrt(2))/x-1)*sqrt(1-(sqrt(2)+2)/x)
(%i5) INCORRECT_AT_2:at(INCORRECT,[x = 2])
(%o5) sqrt((2-sqrt(2))/2-1)*sqrt(1-(sqrt(2)+2)/2)
(%o6) -0.70710678118655
(%o7) "wrong_factorization_of_sqrt.maxima"

At first, as is seen from (%i2) and (%o2),
sqrt(1-1/x)*(2/x-1)) is not factorized to sqrt(...)*sqrt(...)
in any automatic way. This is the correct behavihor.

Next, as is seen from (%i3) and (%o3),
sqrt(1-(2-sqrt(2))/x)*((2+sqrt(2))/x-1)) is factorized to
sqrt((2-sqrt(2))/x-1)*sqrt(1-(sqrt(2)+2)/x) in an automatic way.

This expression is non-negative for 2-sqrt(2) <= x <= 2+sqrt(2).
However, as is demonstrated by (%i4)-(%o5),
the calculated expression takes a negative value at x=2.

Please stop the incorrect factorization of sqrt(X*Y) to sqrt(X)*sqrt(Y).

Sincerely yours,

Satoshi Adachi


  • boud

    boud - 2009-02-19

    There's another bug related to this one IMHO:
    [ 2202764 ] Taylor series of sqrt(1+x-y)
    [This is why i came here - i had a Taylor series bug that i think is

    The bug is solved for me (maxima 5.10.0, debian-etch) by setting

    radexpand : false;

    However, i'm not convinced that this is a sufficient response to the
    bug. i don't really know maxima well enough to know what would be
    most reasonable.

    (1) Set radexpand to false by default? And add more warnings in the
    documentation (or do we expect users to learn mathematics elsewhere
    than in software documentation?)

    (2) Tell maxima not to rewrite sqrt(X) as %i sqrt(-X) when domain:real ?

    (3) Extend the "is X positive, negative or zero?" question to the
    cases where it still needs to be asked but so far does not get asked,
    i.e. when deciding whether

    sqrt(X*Y) = sqrt(-X)*sqrt(-Y)


    sqrt(X*Y) = - sqrt(-X)*sqrt(-Y) ?

    (4) Decide not to invert the order of an expression inside sqrt( ... )
    if domain:real, except when it's sure that the expression is positive?

  • Dieter Kaiser

    Dieter Kaiser - 2009-10-07

    Maxima can simplify the sqrt function correctly.

    No simplification, when the sign of a is not known:

    (%i1) sqrt(a*z);
    (%o1) sqrt(a*z)

    Simplification for a a positive or b a negative value:

    (%i2) assume(a>0)$
    (%i3) sqrt(a*z);
    (%o3) sqrt(a)*sqrt(z)

    (%i4) assume(b<0)$
    (%i5) sqrt(b*z);
    (%o5) sqrt(-b)*sqrt(-z)

    The problem of this bug report is one of the factors of the example. The sign is wrongly determined to be negative:

    (%i8) sign(1-(2-sqrt(2))/x);
    (%o8) neg

    Therefore Maxima does the above simplification, which is wrong.

    Dieter Kaiser

  • Dieter Kaiser

    Dieter Kaiser - 2009-11-11

    For the record. Again two example for a wrong sign:

    (%i3) sign(1-(1+sqrt(2))/x);
    (%o3) neg
    (%i4) sign(1-(1+sqrt(2))*x);
    (%o4) neg

    For these expressions the functions splitprod and splitsum are called. I have not figured out up to now the reason for the bug, but the variable x is lost somewhere and not taken into account to determine the sign.

    Dieter Kaiser

  • Dieter Kaiser

    Dieter Kaiser - 2009-11-11

    I think I have found the bug in the routine compsplt.

    (defun compsplt (x)
    (cond ((atom x) (setq lhs x rhs 0))
    ((atom (car x)) (setq lhs x rhs 0))
    ---> ((not (null (cdr (symbols x)))) (compsplt2 x))
    (t (compsplt1 x))))

    At the marked line Maxima tests the CDR of the result of the function SYMBOLS. SYMBOLS returns a list of all variables of an expression. The list is of the form '($X $Y ...), where $X, $Y are the variables. Because COMPSPLT takes the CDR of SYMBOLS, one symbol in an expression is not recognized. Later in the algorithm the sign of this variable is not taken into account.

    For almost all expression this behavior seems to be no problem, but not for the examples of this bug report.

    When we change the test of the marked line:

    ((not (null (symbols x))) ...

    the sign of the given examples of this bug report are correct:

    (%i3) sign(1-(2+sqrt(2))*x);
    (%o3) pnz
    (%i4) sign(1-(2+sqrt(2))/x);
    (%o4) pnz

    and the reported expression does not longer factorize wrongly:

    (%i5) sqrt((1-(2-sqrt(2))/x)*((2+sqrt(2))/x-1));
    (%o5) sqrt((1-(2-sqrt(2))/x)*((sqrt(2)+2)/x-1))

    One example of the testsuite no longer works. It is one of the new examples for the sign of the abs function. This is the example:

    (%i2) assume(abs(x)<%e/2+sin(1));
    (%o2) [sin(1)+%e/2 > abs(x)]
    (%i3) is(x>-2*%e);
    (%o3) true

    This will be the new result:

    (%i6) assume(abs(x)<%e/2+sin(1));
    (%o6) [-abs(x)+sin(1)+%e/2 > 0]

    (%i7) is(x>-2*%e);
    (%o7) unknown

    This test no longer works, because the ordering of the terms will change. The reason is that COMPSPLT calls splitsum for one variable x in the expression too.

    So, with one exception the testsuite and the share_testsuite have no problems with the suggested change in COMPSPLT.

    Dieter Kaiser

  • Dieter Kaiser

    Dieter Kaiser - 2009-11-14
    • status: open --> closed
  • Dieter Kaiser

    Dieter Kaiser - 2009-11-14

    Fixed in compar.lisp revision 1.62. Closing this bug report as fixed.
    Dieter Kaiser


Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

No, thanks