Fix overflow in logistic-1/2.

diff --git a/src/runtime/arith.scm b/src/runtime/arith.scm
index 54b2322f11bc45a1e616e850f917a725f2ff413c..ec2e183d57ded9a4df4476bd09d679dfbe3de7d4 100644 (file)
--- a/src/runtime/arith.scm
+++ b/src/runtime/arith.scm
@@ -2309,7 +2309,19 @@ USA.
   ;;   2|d2| + 2|d3| + 2|d2 d3| + 2|d1| + 2|d0| + 2|d0 d1|
   ;;   <= 4 eps + 2 eps^2.
   ;;
-  (- (/ (expm1 (- x)) (* 2 (+ 1 (exp (- x)))))))
+  ;; However, if x is a large negative, e^{-x} might overflow.  When
+  ;;
+  ;;    x < log(eps) - 2 < log(eps) - log(4) = log(eps/4)
+  ;;      < log(eps/4) - log(1 - eps/4)
+  ;;      = logit(eps/4),
+  ;;
+  ;; we have logistic(x) < eps/4.  Hence the relative error of
+  ;; logistic(x) - 1/2 from -1/2 is bounded by eps/2, and so the
+  ;; relative error of -1/2 from logistic(x) - 1/2 is bounded by eps.
+  ;;
+  (if (< x (flo:- flo:log-error-bound 2.))
+      -0.5
+      (- (/ (expm1 (- x)) (* 2 (+ 1 (exp (- x))))))))
 
 (define-integrable logit-boundary-lo   ;logistic(-1)
   (flo:/ (flo:exp -1.) (flo:+ 1. (flo:exp -1.))))

diff --git a/tests/runtime/test-arith.scm b/tests/runtime/test-arith.scm
index f43d0aabf2dd9a25750366556ee65ffcfbc52729..31f87ee4c4c0371df24148c52c209611af21a7ba 100644 (file)
--- a/tests/runtime/test-arith.scm
+++ b/tests/runtime/test-arith.scm
@@ -401,11 +401,7 @@ USA.
             (assert-<= (relerr x (logit1/2+ p)) 1e-15)
             (assert-= (- (logit1/2+ p)) (logit1/2+ (- p)))))
       (assert-<= (relerr p (logistic-1/2 x)) 1e-15)
-      ;; Currently logistic-1/2 takes no measures to avoid overflow in
-      ;; (exp (- x)).
-      (with-expected-failure (if (>= x 710) 'xfail #f)
-        (lambda ()
-          (assert-= (- (logistic-1/2 x)) (logistic-1/2 (- x))))))))
+      (assert-= (- (logistic-1/2 x)) (logistic-1/2 (- x))))))
 
 (define-enumerated-test 'expt-exact
   (vector