Fix #147: xlapy2 not propagating nans

xLAPY2 now returns a NaN whenever input variables X or Y are NaNs.

The previous xLAPY2 was relying on FORTRAN INTRINSIC MAX and MIN to behave in a
certain way with NaNs (i.e. return a NaN whenever X or Y are NaN) and this
behavior is not observed on some (most?) compilers.

We handle the NaN behavior of xLAPY2 by checking for NaNs at the start of the
function.

Thanks to Andreas Noack for providing report and sample code to demonstrate the
problem.

diff --git a/SRC/dlapy2.f b/SRC/dlapy2.f
index 3861b1d..b40c46a 100644 (file)
--- a/SRC/dlapy2.f
+++ b/SRC/dlapy2.f
@@ -82,20 +82,32 @@
 *     ..
 *     .. Local Scalars ..
       DOUBLE PRECISION   W, XABS, YABS, Z
+      LOGICAL            X_IS_NAN, Y_IS_NAN
+*     ..
+*     .. External Functions ..
+      LOGICAL            DISNAN
+      EXTERNAL           DISNAN
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          ABS, MAX, MIN, SQRT
 *     ..
 *     .. Executable Statements ..
 *
-      XABS = ABS( X )
-      YABS = ABS( Y )
-      W = MAX( XABS, YABS )
-      Z = MIN( XABS, YABS )
-      IF( Z.EQ.ZERO ) THEN
-         DLAPY2 = W
-      ELSE
-         DLAPY2 = W*SQRT( ONE+( Z / W )**2 )
+      X_IS_NAN = DISNAN( X )
+      Y_IS_NAN = DISNAN( Y )
+      IF ( X_IS_NAN ) DLAPY2 = X
+      IF ( Y_IS_NAN ) DLAPY2 = Y
+*
+      IF ( .NOT.( X_IS_NAN.OR.Y_IS_NAN ) ) THEN
+         XABS = ABS( X )
+         YABS = ABS( Y )
+         W = MAX( XABS, YABS )
+         Z = MIN( XABS, YABS )
+         IF( Z.EQ.ZERO ) THEN
+            DLAPY2 = W
+         ELSE
+            DLAPY2 = W*SQRT( ONE+( Z / W )**2 )
+         END IF
       END IF
       RETURN
 *

diff --git a/SRC/slapy2.f b/SRC/slapy2.f
index 13e2198..a8ee8b5 100644 (file)
--- a/SRC/slapy2.f
+++ b/SRC/slapy2.f
@@ -82,20 +82,35 @@
 *     ..
 *     .. Local Scalars ..
       REAL               W, XABS, YABS, Z
+      LOGICAL            X_IS_NAN, Y_IS_NAN
+*     ..
+*     .. External Functions ..
+      LOGICAL            SISNAN
+      EXTERNAL           SISNAN
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          ABS, MAX, MIN, SQRT
 *     ..
 *     .. Executable Statements ..
 *
-      XABS = ABS( X )
-      YABS = ABS( Y )
-      W = MAX( XABS, YABS )
-      Z = MIN( XABS, YABS )
-      IF( Z.EQ.ZERO ) THEN
-         SLAPY2 = W
-      ELSE
-         SLAPY2 = W*SQRT( ONE+( Z / W )**2 )
+*     ..
+*     .. Executable Statements ..
+*
+      X_IS_NAN = SISNAN( X )
+      Y_IS_NAN = SISNAN( Y )
+      IF ( X_IS_NAN ) SLAPY2 = X
+      IF ( Y_IS_NAN ) SLAPY2 = Y
+*
+      IF ( .NOT.( X_IS_NAN.OR.Y_IS_NAN ) ) THEN
+         XABS = ABS( X )
+         YABS = ABS( Y )
+         W = MAX( XABS, YABS )
+         Z = MIN( XABS, YABS )
+         IF( Z.EQ.ZERO ) THEN
+            SLAPY2 = W
+         ELSE
+            SLAPY2 = W*SQRT( ONE+( Z / W )**2 )
+         END IF
       END IF
       RETURN
 *