[Maxima-bugs] [ maxima-Bugs-1981623 ] wrong sign of sqrt()

SourceForge Headquarters 225 Broadway Suite 1600 San Diego, CA 92101 +1 (858) 454-5900

Bugs item #1981623, was opened at 2008-06-01 23:39
Message generated for change (Comment added) made by dgildea
You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=104933&aid=1981623&group_id=4933

Please note that this message will contain a full copy of the comment thread,
including the initial issue submission, for this request,
not just the latest update.
Category: Lisp Core - Simplification
Group: None
>Status: Closed
>Resolution: Fixed
Priority: 5
Private: No
Submitted By: Satoshi Adachi (satoshi_adachi)
>Assigned to: Dan Gildea (dgildea)
Summary: wrong sign of sqrt()

Initial Comment:
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

----------------------------------------------------------------------

>Comment By: Dan Gildea (dgildea)
Date: 2009-03-15 17:18

Message:
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)

----------------------------------------------------------------------

Comment By: Dan Gildea (dgildea)
Date: 2009-02-15 16:06

Message:
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".

----------------------------------------------------------------------

Comment By: Raymond Toy (rtoy)
Date: 2008-07-30 14:15

Message:
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.

----------------------------------------------------------------------

Comment By: Robert Dodier (robert_dodier)
Date: 2008-06-23 10:26

Message:
Logged In: YES 
user_id=501686
Originator: NO

Assign category = lisp core / simplification.

----------------------------------------------------------------------

Comment By: Satoshi Adachi (satoshi_adachi)
Date: 2008-06-11 02:50

Message:
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...

----------------------------------------------------------------------

Comment By: Barton Willis (willisbl)
Date: 2008-06-04 07:20

Message:
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)

----------------------------------------------------------------------

You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=104933&aid=1981623&group_id=4933

[Maxima-bugs] [ maxima-Bugs-1981623 ] wrong sign of sqrt()

Computer Algebra System written in Common Lisp

[Maxima-bugs] [ maxima-Bugs-1981623 ] wrong sign of sqrt()