From 5e272e0ac21298dba1225c8d6837ed15671d7872 Mon Sep 17 00:00:00 2001 From: Taylor R Campbell Date: Tue, 20 Nov 2018 06:34:38 +0000 Subject: [PATCH] Fix overflow in logistic-1/2. --- src/runtime/arith.scm | 14 +++++++++++++- tests/runtime/test-arith.scm | 6 +----- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/src/runtime/arith.scm b/src/runtime/arith.scm index 54b2322f1..ec2e183d5 100644 --- 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 f43d0aabf..31f87ee4c 100644 --- 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 -- 2.25.1