* Bug #13111 fixed - Elementary_functions: fix 0 and -0 for sqrt. 00/13300/1
Paul Bignier [Mon, 2 Dec 2013 10:07:25 +0000 (11:07 +0100)]
sqrt returned different results when imaginary part was -0 versus 0.

Change-Id: Ic82dd9b0f2e8a155fac815df89f782820db934da

scilab/CHANGES_5.5.X
scilab/modules/elementary_functions/src/fortran/wsqrt.f
scilab/modules/elementary_functions/tests/nonreg_tests/bug_13111.dia.ref [new file with mode: 0644]
scilab/modules/elementary_functions/tests/nonreg_tests/bug_13111.tst [new file with mode: 0644]

index ace2884..493e319 100644 (file)
@@ -275,6 +275,8 @@ Scilab Bug Fixes
 
 * Bug #13109 fixed - pll2str now supports polynomials with complex coefficients and hypermatrices.
 
+* Bug #13111 fixed - sqrt returned different results when imaginary part was -0 versus 0.
+
 * Bug #13114 fixed - clear_pixmap should have been removed in Scilab 5.4.1.
 
 * Bug #13116 fixed - qpsolve now respects upper-bounds constraints.
index e9d797a..32898d3 100644 (file)
@@ -1,10 +1,10 @@
 c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 c Copyright (C) Bruno Pincon
-c 
+c
 c This file must be used under the terms of the CeCILL.
 c This source file is licensed as described in the file COPYING, which
 c you should have received as part of this distribution.  The terms
-c are also available at    
+c are also available at
 c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
 
       subroutine wsqrt(xr,xi,yr,yi)
@@ -23,7 +23,7 @@ c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
 *
 *     ALGORITHM
 *        essentially the classic one which consists in
-*        choosing the good formula such as avoid cancellation ; 
+*        choosing the good formula such as avoid cancellation ;
 *        Also rare spurious overflow are treated with a
 *        "manual" method. For some more "automated" methods
 *        (but currently difficult to implement in a portable
@@ -33,20 +33,20 @@ c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
 *          "Implementing Complex Elementary Functions Using
 *          Exception Handling", ACM TOMS, Vol. 20 (1994), pp 215-244
 *
-*        for xr > 0 :                            
+*        for xr > 0 :
 *          yr = sqrt(2( xr + sqrt( xr^2 + xi^2)) )/ 2
 *          yi = xi / sqrt(2(xr + sqrt(xr^2 + xi^2)))
 *
-*        and for xr < 0 :                           
+*        and for xr < 0 :
 *          yr = |xi| / sqrt( 2( -xr + sqrt( xr^2 + xi^2 )) )
 *          yi = sign(xi) sqrt(2(-xr + sqrt( xr^2 + xi^2))) / 2
-*     
-*        for xr = 0 use          
+*
+*        for xr = 0 use
 *          yr = sqrt(0.5)*sqrt(|xi|)  when |xi| is such that 0.5*|xi| may underflow
 *             = sqrt(0.5*|xi|)        else
 *          yi = sign(xi) yr
 *
-*        Noting t = sqrt( 2( |xr| + sqrt( xr^2 + yr^2)) ) 
+*        Noting t = sqrt( 2( |xr| + sqrt( xr^2 + yr^2)) )
 *                 = sqrt( 2( |xr| + pythag(xr,xi) ) )
 *        it comes :
 *
@@ -54,9 +54,9 @@ c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
 *         --------------+---------------------
 *           yr = 0.5*t  |  yr =  |xi| / t
 *           yi = xi / t |  yi = sign(xi)*0.5* t
-*    
-*        as the function pythag must not underflow (and overflow only 
-*        if sqrt(x^2+y^2) > rmax) only spurious (rare) case of overflow 
+*
+*        as the function pythag must not underflow (and overflow only
+*        if sqrt(x^2+y^2) > rmax) only spurious (rare) case of overflow
 *        occurs in which case a scaling is done.
 *
 *     AUTHOR
@@ -71,7 +71,7 @@ c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
       double precision a, b, t
 *     EXTERNAL
       double precision pythag, dlamch
-      external         pythag, dlamch 
+      external         pythag, dlamch
       integer          isanan
       external         isanan
 
@@ -81,7 +81,7 @@ c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
       save    first
       data    first /.true./
       save             RMAX, BRMIN
-      
+
 
       if (first) then
          RMAX = dlamch('O')
@@ -112,10 +112,14 @@ c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
             yi = b/t
          else
             yr = abs(b)/t
-            yi = sign(0.5d0,b)*t
+            if ( b .ge. 0 ) then
+               yi = 0.5d0*t
+            else
+               yi = -0.5d0*t
+            endif
          endif
       else
-*        Here we treat the special cases where a and b are +- 00 or NaN. 
+*        Here we treat the special cases where a and b are +- 00 or NaN.
 *        The following is the treatment recommended by the C99 standard
 *        with the simplification of returning NaN + i NaN if the
 *        the real part or the imaginary part is NaN (C99 recommends
@@ -131,7 +135,11 @@ c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
          elseif ( a .lt. -RMAX ) then
 *           here a is -Inf and b is finite
             yr = 0.d0
-            yi = sign(1.d0,b)*abs(a)
+            if ( b .ge. 0 ) then
+               yi = 1.d0*abs(a)
+            else
+               yi = -1.d0*abs(a)
+            endif
          else
 *           here a is +Inf and b is finite
             yr = a
@@ -151,7 +159,11 @@ c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
          yi = 4.d0*b/t
       else
          yr = 4.d0*abs(b)/t
-         yi = sign(2.d0,b)*t
+         if ( b .ge. 0 ) then
+            yi = 2.d0*t
+         else
+            yi = -2.d0*t
+         endif
       endif
 
       end
diff --git a/scilab/modules/elementary_functions/tests/nonreg_tests/bug_13111.dia.ref b/scilab/modules/elementary_functions/tests/nonreg_tests/bug_13111.dia.ref
new file mode 100644 (file)
index 0000000..f36838c
--- /dev/null
@@ -0,0 +1,20 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 13111 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=13111
+//
+// <-- Short Description -->
+// sqrt returned different results when imaginary part was -0 versus 0.
+assert_checkequal(sqrt(-4-0*%i), 2*%i);
+assert_checkequal(sqrt(-4-0*%i), sqrt(-4+0*%i));
+assert_checkequal(sqrt(-%inf-0*%i), complex(0, %inf));
+assert_checkequal(sqrt(-%inf-0*%i), sqrt(-%inf+0*%i));
diff --git a/scilab/modules/elementary_functions/tests/nonreg_tests/bug_13111.tst b/scilab/modules/elementary_functions/tests/nonreg_tests/bug_13111.tst
new file mode 100644 (file)
index 0000000..89ded30
--- /dev/null
@@ -0,0 +1,22 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 13111 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=13111
+//
+// <-- Short Description -->
+// sqrt returned different results when imaginary part was -0 versus 0.
+
+assert_checkequal(sqrt(-4-0*%i), 2*%i);
+assert_checkequal(sqrt(-4-0*%i), sqrt(-4+0*%i));
+
+assert_checkequal(sqrt(-%inf-0*%i), complex(0, %inf));
+assert_checkequal(sqrt(-%inf-0*%i), sqrt(-%inf+0*%i));