#1426 wrong sign of sqrt()

closed
5
2009-03-15
2008-06-02
No

Dear Developers of Maxima,

I found that sqrt() returns the incorrect expression that has the sign
opposite to the ture expression when some certain argument is given to
sqrt(). Namely, sqrt() interprets incorrectly the database that is prepared
by assume().

Here is an demonstration program:
-------------------------------------------------------------------------------
/*
* wrong_sign_of_sqrt.maxima
*
* S.Adachi 2008/06/01
*/

display2d:false;

assume(x >= 0, x <= 1);

correct_result_1:sqrt((x-1)^2);

correct_result_2:sqrt(1/(x-1)^2);

correct_result_3:sqrt(a*(x-1)^2);

incorrect_result_1:sqrt(a/(x-1)^2);

incorrect_result_2:sqrt(a^2/(x-1)^2);

incorrect_result_3:sqrt(x^2/(x-1)^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_sign_of_sqrt.maxima)

batching #p/Volumes/HFS+2/home/adachi/work/299/wrong_sign_of_sqrt.maxima
(%i2) display2d : false
(%o2) false
(%i3) assume(x >= 0,x <= 1)
(%o3) [x >= 0,x <= 1]
(%i4) correct_result_1:sqrt((x-1)^2)
(%o4) 1-x
(%i5) correct_result_2:sqrt(1/(x-1)^2)
(%o5) 1/(1-x)
(%i6) correct_result_3:sqrt(a*(x-1)^2)
(%o6) sqrt(a)*(1-x)
(%i7) incorrect_result_1:sqrt(a/(x-1)^2)
(%o7) sqrt(a)/(x-1)
(%i8) incorrect_result_2:sqrt(a^2/(x-1)^2)
(%o8) abs(a)/(x-1)
(%i9) incorrect_result_3:sqrt(x^2/(x-1)^2)
(%o9) x/(x-1)
(%o10) "wrong_sign_of_sqrt.maxima"
-------------------------------------------------------------------------------

I wonder why sqrt() returns the wrong expression if sqrt((x-1)^2) appears
in the denominator of some fraction in the argument that is more complex
than some threshold (e.g. the numerator is not a simple number or something
like that).

I think that this is a very severe bug of sqrt()
and the database that is prepared by assume().
This bug puts many user programs to the state
producing many wrong results.

Please fix this severe bug.

Sincerely yours,
Satoshi Adachi

Discussion

  • Barton Willis

    Barton Willis - 2008-06-04

    Logged In: YES
    user_id=895922
    Originator: NO

    As a workaround, you might try setting radexpand to false:

    (%i1) assume(0 < x, x <= 1)$
    (%i2) sqrt(a/(x-1)^2);
    (%o2) sqrt(a)/(x-1)
    (%i3) sqrt(a/(x-1)^2), radexpand : false;
    (%o3) sqrt(a/(x-1)^2)

     
  • Satoshi Adachi

    Satoshi Adachi - 2008-06-11

    Logged In: YES
    user_id=1953419
    Originator: YES

    Thank you for your suggestion.
    But, your suggestion just forbids Maxima to simplify certain expressions
    in order not to produce wrong results.
    Is someone trying to fix this bug now?
    If not yet, I will read lisp source code in some future (maybe, several months
    later) to try to fix the problem...

     
  • Robert Dodier

    Robert Dodier - 2008-06-23
    • labels: --> Lisp Core - Simplification
     
  • Robert Dodier

    Robert Dodier - 2008-06-23

    Logged In: YES
    user_id=501686
    Originator: NO

    Assign category = lisp core / simplification.

     
  • Raymond Toy

    Raymond Toy - 2008-07-30

    Logged In: YES
    user_id=28849
    Originator: NO

    FWIW, the code that handles this is simpexpt in simp.lisp. For the case of sqrt((x-1)^2), the e1 tag is done, and the second cond case is executed. The result is correct in this case.

    However, for sqrt(1/(x-1)^2), the e1 tag is also done, but the first cond case is executed because the (noneg (cadr gr)) returns non-nil. In the former case the noneg call returns NIL, which is correct because x-1 is not negative.

    I don't know why the two different calls to noneg return different values for x-1. NIL should be returned in each case because x <= 1.

    There is one difference, though. The z arg to simpexpt is NIL for the case that works and T for the case that is incorrect.

     
  • Dan Gildea

    Dan Gildea - 2009-02-15

    Tracing dcompare shows that the results of comparisons with x
    change in the middle of the computation.

    (%i4) sqrt(a/(x-1)^2);
    0: (DCOMPARE $X 1)
    0: DCOMPARE returned $NZ
    0: (DCOMPARE $X 1.0)
    0: DCOMPARE returned $PNZ
    0: (DCOMPARE $X 0)
    0: DCOMPARE returned $PZ
    0: (DCOMPARE $X 1)
    0: DCOMPARE returned $POS

    Similar behavior observed with bug 1419046 "sign bug".

     
  • Dan Gildea

    Dan Gildea - 2009-03-15
    • status: open --> closed
    • assigned_to: nobody --> dgildea
     
  • Dan Gildea

    Dan Gildea - 2009-03-15

    Fixed in compar.lisp rev 1.46.

    (%i2) assume(x >= 0, x <= 1);
    (%o2) [x >= 0,x <= 1]

    (%i3) correct_result_1:sqrt((x-1)^2);
    (%o3) 1-x

    (%i4) correct_result_2:sqrt(1/(x-1)^2);
    (%o4) 1/(1-x)

    (%i5) correct_result_3:sqrt(a*(x-1)^2);
    (%o5) sqrt(a)*(1-x)

    (%i6) incorrect_result_1:sqrt(a/(x-1)^2);
    (%o6) sqrt(a)/(1-x)

    (%i7) incorrect_result_2:sqrt(a^2/(x-1)^2);
    (%o7) abs(a)/(1-x)

    (%i8) incorrect_result_3:sqrt(x^2/(x-1)^2);
    (%o8) x/(1-x)

     

Log in to post a comment.