Mask underflow exception in logsumexp because it doesn't matter.
authorTaylor R Campbell <campbell@mumble.net>
Tue, 30 Oct 2018 15:59:04 +0000 (15:59 +0000)
committerTaylor R Campbell <campbell@mumble.net>
Tue, 30 Oct 2018 15:59:04 +0000 (15:59 +0000)
src/runtime/arith.scm
tests/runtime/test-arith.scm

index a129257fd0cb8c08df4fadbe267c24ca2a4a8f0e..14dd4a00ee381a15c022be1507b22c7aee8ab45f 100644 (file)
@@ -2064,7 +2064,18 @@ USA.
                     (not (= (- m) (reduce min #f l)))
                     (not (any nan? l)))
                m
-               (+ m (log (reduce + 0 (map (lambda (x) (exp (- x m))) l))))))
+               ;; Overflow is not possible because everything is
+               ;; normalized to be below zero.  Underflow can be
+               ;; safely ignored because it can't change the outcome:
+               ;; even if you had 2^64 copies of the largest subnormal
+               ;; in the sum, 2^64 * largest subnormal < 2^900 <<<
+               ;; epsilon = 2^-53, and at least one addend in the sum
+               ;; is 1.
+               (flo:with-exceptions-untrapped (flo:exception:underflow)
+                 (lambda ()
+                   (+ m
+                      (log
+                       (reduce + 0 (map (lambda (x) (exp (- x m))) l))))))))
          (car l))
       (flo:-inf.0)))
 
index ebce10c033909ab0f19cb47aa9203f5b949eddb5..5237363576fc6d93f494b52ad7dce547165980ce 100644 (file)
@@ -192,6 +192,7 @@ USA.
 
 (define-enumerated-test 'logsumexp-values
   (vector
+   (vector (iota 1000) 999.45867514538713)
    (vector '(999 1000) 1000.3132616875182)
    (vector '(-1000 -1000) (+ -1000 (log 2)))
    (vector '(0 0) (log 2)))