Omit unnecessary case for logistic. Add a test for this case.

The relative error is small in naive evaluation of 1/(1 + e^{-x})
even if x is very negative.

diff --git a/src/runtime/arith.scm b/src/runtime/arith.scm
index 658d98c1b36a8f940d5fb8b95c3f2d6381ba7d9a..134e9990851be813d8e62bea1a929744ae0b34f2 100644 (file)
--- a/src/runtime/arith.scm
+++ b/src/runtime/arith.scm
@@ -2078,44 +2078,12 @@ USA.
          ;;      < eps.
          ;;
          (exp x))
-        ((< x 0)
-         ;; e^x < 1 < 1 + e^x, so e^x/(1 + e^x) < 1, so, if exp has
-         ;; relative error d0, + has relative error d1, and / has
-         ;; relative error d2, then we get
-         ;;
-         ;;     (1 + d2) (1 + d0) e^x/[(1 + (1 + d0) e^x)(1 + d1)]
-         ;;     = (1 + d0 + d2 + d0 d2) e^x/[(1 + (1 + d0) e^x)(1 + d1)]
-         ;;     = (1 + d') e^x/[(1 + (1 + d0) e^x)(1 + d1)]
-         ;;     = (1 + d') e^x/[(1 + e^x + d0 e^x)(1 + d1)]
-         ;;     = (1 + d') e^x/[1 + e^x + d0 e^x + d1 + d1 e^x + d1 d0 e^x]
-         ;;     = (1 + d') e^x/[1 + e^x + d0 e^x + d1 (1 + e^x) + d1 d0 e^x]
-         ;;     = (1 + d') e^x/[(1 + e^x)(1 + d0 e^x/(1 + e^x)
-         ;;                                 + d1
-         ;;                                 + d1 d0 e^x/(1 + e^x))]
-         ;;     = (1 + d') e^x/[(1 + e^x)(1 + d'')]
-         ;;     = [e^x/(1 + e^x)] (1 + d')/(1 + d'')
-         ;;
-         ;; where
-         ;;
-         ;;     d' = d0 + d2 + d0 d2,
-         ;;     d'' = d0 e^x/(1 + e^x) + d1 + d1 d0 e^x/(1 + e^x).
-         ;;
-         ;; Note that since e^x/(1 + e^x) < 1,
-         ;;
-         ;;     |d''| <= |d0| + |d1| + |d1 d0| <= 2 eps + eps^2,
-         ;;
-         ;; so by Lemma 2 the relative error is bounded by
-         ;;
-         ;;     2|d' - d''|
-         ;;     <= 2|d0| + 2|d2| + 2|d0 d2| + 2|d0| + 2|d1| + 2|d1 d0|
-         ;;     <= 8 eps + 4 eps^2.
-         ;;
-         (let ((e^x (exp x)))
-           (/ e^x (+ 1 e^x))))
        ((<= x 37)
-         ;; e^{-x} < 1 < 1 + e^{-x}, so 1/(1 + e^{-x}) < 1, so, if exp
-         ;; has relative error d0, + has relative error d1, and / has
-         ;; relative error d2, then we get
+         ;; e^{-x} > 0, so 1 + e^{-x} > 1, and 0 < 1/(1 + e^{-x}) < 1;
+         ;; further, since e^{-x} < 1 + e^{-x}, we also have 0 <
+         ;; e^{-x}/(1 + e^{-x}) < 1.  Thus, if exp has relative error
+         ;; d0, + has relative error d1, and / has relative error d2,
+         ;; then we get
          ;;
          ;;     (1 + d2)/[(1 + (1 + d0) e^{-x})(1 + d1)]
          ;;     = (1 + d2)/[1 + e^{-x} + d0 e^{-x}

diff --git a/tests/runtime/test-arith.scm b/tests/runtime/test-arith.scm
index df7eddcab4bb73e189de4b933f348991cd771f18..7017c7e0e76e8ea302925666eb6a072fe721c2c4 100644 (file)
--- a/tests/runtime/test-arith.scm
+++ b/tests/runtime/test-arith.scm
@@ -192,6 +192,9 @@ USA.
 
 (define-enumerated-test 'logit-logistic
   (vector
+   (vector -36.7368005696771
+           1.1102230246251565e-16
+           (log 1.1102230246251565e-16))
    (vector -1.0000001 0.2689414017088022 (log .2689414017088022))
    (vector -0.9999999 0.26894144103118883 (log .26894144103118883))
    (vector -1 0.2689414213699951 (log 0.2689414213699951))