Mercurial > octave-nkf

changeset 14392:d17237256856

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make SLATEC-FN atanh/erfc functions more tolerant about edge cases like numerical underflow or NaN values. * slatec-fn/atanh,f (ATANH): Returns infinity for +-1 and NaN for >1 * slatec-fn/datanh.f (DATANH): Likewise. * slatec-fn/erfc.f (ERFC): Returns NaN for Nan input. * slatec-fn/derfc.f (DERFC): Likewise.
author Michael Goffioul <michael.goffioul@gmail.com>
date Thu, 23 Feb 2012 09:12:47 +0000
parents c9ec21bef97a
children e41e538e9d03
files libcruft/slatec-fn/atanh.f libcruft/slatec-fn/datanh.f libcruft/slatec-fn/derfc.f libcruft/slatec-fn/erfc.f
diffstat 4 files changed, 28 insertions(+), 8 deletions(-) [+]
line wrap: on
line diff
--- a/libcruft/slatec-fn/atanh.f	Wed Feb 22 20:06:05 2012 -0500
+++ b/libcruft/slatec-fn/atanh.f	Thu Feb 23 09:12:47 2012 +0000
@@ -57,8 +57,14 @@
       FIRST = .FALSE.
 C
       Y = ABS(X)
-      IF (Y .GE. 1.0) CALL XERMSG ('SLATEC', 'ATANH', 'ABS(X) GE 1', 2,
-     +   2)
+      IF (Y .GE. 1.0) THEN
+         IF (Y .GT. 1.0) THEN 
+            ATANH = (X - X) / (X - X)
+         ELSE
+            ATANH = X / 0.0
+         ENDIF
+         RETURN
+      ENDIF
 C
       IF (1.0-Y .LT. DXREL) CALL XERMSG ('SLATEC', 'ATANH',
      +   'ANSWER LT HALF PRECISION BECAUSE ABS(X) TOO NEAR 1', 1, 1)
--- a/libcruft/slatec-fn/datanh.f	Wed Feb 22 20:06:05 2012 -0500
+++ b/libcruft/slatec-fn/datanh.f	Thu Feb 23 09:12:47 2012 +0000
@@ -68,8 +68,14 @@
       FIRST = .FALSE.
 C
       Y = ABS(X)
-      IF (Y .GE. 1.D0) CALL XERMSG ('SLATEC', 'DATANH', 'ABS(X) GE 1',
-     +   2, 2)
+      IF (Y .GE. 1.D0) THEN
+         IF (Y .GT. 1.D0) THEN 
+            DATANH = (X - X) / (X - X)
+         ELSE
+            DATANH = X / 0.D0
+         ENDIF
+         RETURN
+      ENDIF
 C
       IF (1.D0-Y .LT. DXREL) CALL XERMSG ('SLATEC', 'DATANH',
      +   'ANSWER LT HALF PRECISION BECAUSE ABS(X) TOO NEAR 1', 1, 1)
--- a/libcruft/slatec-fn/derfc.f	Wed Feb 22 20:06:05 2012 -0500
+++ b/libcruft/slatec-fn/derfc.f	Thu Feb 23 09:12:47 2012 +0000
@@ -191,6 +191,11 @@
       ENDIF
       FIRST = .FALSE.
 C
+      IF (ISNAN(X)) THEN
+         DERFC = X
+         RETURN
+      ENDIF
+C
       IF (X.GT.XSML) GO TO 20
 C
 C ERFC(X) = 1.0 - ERF(X)  FOR  X .LT. XSML
@@ -219,8 +224,7 @@
       IF (X.LT.0.D0) DERFC = 2.0D0 - DERFC
       RETURN
 C
- 40   CALL XERMSG ('SLATEC', 'DERFC', 'X SO BIG ERFC UNDERFLOWS', 1, 1)
-      DERFC = 0.D0
+ 40   DERFC = 0.D0
       RETURN
 C
       END
--- a/libcruft/slatec-fn/erfc.f	Wed Feb 22 20:06:05 2012 -0500
+++ b/libcruft/slatec-fn/erfc.f	Thu Feb 23 09:12:47 2012 +0000
@@ -121,6 +121,11 @@
       ENDIF
       FIRST = .FALSE.
 C
+      IF (ISNAN(X)) THEN
+         ERFC = X
+         RETURN
+      ENDIF
+C
       IF (X.GT.XSML) GO TO 20
 C
 C ERFC(X) = 1.0 - ERF(X) FOR X .LT. XSML
@@ -149,8 +154,7 @@
       IF (X.LT.0.) ERFC = 2.0 - ERFC
       RETURN
 C
- 40   CALL XERMSG ('SLATEC', 'ERFC', 'X SO BIG ERFC UNDERFLOWS', 1, 1)
-      ERFC = 0.
+ 40   ERFC = 0.
       RETURN
 C
       END