From eb9cd33fe73c202a449d1154c816465a98cecea3 Mon Sep 17 00:00:00 2001
From: Taylor R Campbell <campbell@mumble.net>
Date: Sun, 21 Oct 2018 05:43:17 +0000
Subject: [PATCH] Define log1p and expm1 with no range limitations.

---
 src/runtime/arith.scm        | 27 +++++++++++++++++++++
 src/runtime/runtime.pkg      |  2 ++
 tests/runtime/test-arith.scm | 47 +++++++++++++++++++++++++++++++++++-
 3 files changed, 75 insertions(+), 1 deletion(-)

diff --git a/src/runtime/arith.scm b/src/runtime/arith.scm
index 97304bbf7..1beb8fff6 100644
--- a/src/runtime/arith.scm
+++ b/src/runtime/arith.scm
@@ -1175,6 +1175,8 @@ USA.
 
 (define-transcendental-unary real:exp real:exact0= 1 flo:exp)
 (define-transcendental-unary real:log real:exact1= 0 flo:log)
+(define-transcendental-unary real:expm1 real:exact0= 0 flo:expm1-guarded)
+(define-transcendental-unary real:log1p real:exact0= 0 flo:log1p-guarded)
 (define-transcendental-unary real:sin real:exact0= 0 flo:sin)
 (define-transcendental-unary real:cos real:exact0= 1 flo:cos)
 (define-transcendental-unary real:tan real:exact0= 0 flo:tan)
@@ -1182,6 +1184,21 @@ USA.
 (define-transcendental-unary real:acos real:exact1= 0 flo:acos)
 (define-transcendental-unary real:atan real:exact0= 0 flo:atan)
 
+(define-integrable flo:log2 (flo:log 2.))
+(define-integrable flo:1-sqrt1/2 (flo:- 1. (flo:sqrt 0.5)))
+
+(declare (integrate flo:expm1-guarded))
+(define (flo:expm1-guarded x)
+  (if (flo:< (flo:abs x) flo:log2)
+      (flo:expm1 x)
+      (flo:- (flo:exp x) 1.)))
+
+(declare (integrate flo:log1p-guarded))
+(define (flo:log1p-guarded x)
+  (if (flo:< (flo:abs x) flo:1-sqrt1/2)
+      (flo:log1p x)
+      (flo:log (flo:+ 1. x))))
+
 (define (real:atan2 y x)
   (if (and (real:exact0= y)
 	   (real:exact? x))
@@ -1643,6 +1660,16 @@ USA.
 	(else
 	 ((copy real:log) z))))
 
+(define (complex:expm1 z)
+  (if (recnum? z)
+      (complex:- (complex:exp z) 1)	;XXX
+      ((copy real:expm1) z)))
+
+(define (complex:log1p z)
+  (if (recnum? z)
+      (complex:log (complex:+ z 1))	;XXX
+      ((copy real:log1p) z)))
+
 (define (complex:sin z)
   (if (recnum? z)
       (complex:/ (let ((iz (complex:+i* z)))
diff --git a/src/runtime/runtime.pkg b/src/runtime/runtime.pkg
index 7f9c8abb4..693ff665d 100644
--- a/src/runtime/runtime.pkg
+++ b/src/runtime/runtime.pkg
@@ -3299,6 +3299,7 @@ USA.
 	  (exact-rational? rat:rational?)
 	  (exact? complex:exact?)
 	  (exp complex:exp)
+	  (expm1 complex:expm1)
 	  (expt complex:expt)
 	  (finite? complex:finite?)
 	  (floor complex:floor)
@@ -3314,6 +3315,7 @@ USA.
 	  (integer-truncate complex:quotient)
 	  (integer? complex:integer?)
 	  (log complex:log)
+	  (log1p complex:log1p)
 	  (magnitude complex:magnitude)
 	  (make-polar complex:make-polar)
 	  (make-rectangular complex:make-rectangular)
diff --git a/tests/runtime/test-arith.scm b/tests/runtime/test-arith.scm
index f781f3c97..290966bc7 100644
--- a/tests/runtime/test-arith.scm
+++ b/tests/runtime/test-arith.scm
@@ -115,4 +115,49 @@ USA.
 
 (define-enumerated-test 'X/NAN
   (vector (flo:-inf.0) -2. -1 -0. 0 +0. +1 +2. (flo:+inf.0))
-  (lambda (v) (assert-nan (/ v (flo:nan.0)))))
\ No newline at end of file
+  (lambda (v) (assert-nan (/ v (flo:nan.0)))))
+
+(define-enumerated-test 'log1p-exact
+  (vector
+   (cons 0 0)
+   (cons 1 (log 2))
+   (cons (- (exp 1) 1) 1.)
+   (cons (expt 2 -53) (expt 2. -53))
+   (cons (- (expt 2 -53)) (- (expt 2. -53))))
+  (lambda (v)
+    (assert-eqv (cdr v) (log1p (car v)))))
+
+(define-enumerated-test 'expm1-exact
+  (vector
+   (cons 0 0)
+   (cons (log 2) 1.)
+   (cons 1 (- (exp 1) 1))
+   (cons (expt 2 -53) (expt 2. -53))
+   (cons (- (expt 2 -53)) (- (expt 2. -53))))
+  (lambda (v)
+    (assert-eqv (cdr v) (expm1 (car v)))))
+
+(define (relerr a e)
+  (abs (/ (- a e) a)))
+
+(define-enumerated-test 'expm1-approx
+  (vector
+   (cons -0.7 -0.5034146962085905)
+   (cons (- (log 2)) -0.5)
+   (cons -0.6 -0.45118836390597356)
+   (cons 0.6 .8221188003905089)
+   (cons (log 2) 1.)
+   (cons 0.7 1.0137527074704766))
+  (lambda (v)
+    (assert->= 1e-15 (relerr (cdr v) (expm1 (car v))))))
+
+(define-enumerated-test 'log1p-approx
+  (vector
+   (cons -0.3 -.35667494393873245)
+   (cons (- (sqrt 1/2) 1) -0.34657359027997264)
+   (cons -0.25 -.2876820724517809)
+   (cons 0.25 .22314355131420976)
+   (cons (- 1 (sqrt 1/2)) 0.25688251232181475)
+   (cons 0.3 .26236426446749106))
+  (lambda (v)
+    (assert->= 1e-15 (relerr (cdr v) (log1p (car v))))))
\ No newline at end of file
-- 
2.25.1