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.

Author: langou

SRC/dlapy2.f

@@ -82,20 +82,32 @@ DOUBLE PRECISION FUNCTION DLAPY2( X, Y )
 *     ..
 *     .. 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
 *

SRC/slapy2.f

@@ -82,20 +82,35 @@ REAL             FUNCTION SLAPY2( X, Y )
 *     ..
 *     .. 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
 *