Define log1p and expm1 with no range limitations.

diff --git a/src/runtime/arith.scm b/src/runtime/arith.scm
index 97304bbf75ed330d3d87c4a56007df4f78469a53..1beb8fff63e7990269025431a934c650f4f937e2 100644 (file)
--- 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 7f9c8abb430d976f97b0bf3b286955b2b266209e..693ff665d9a590f640b5924287ce31ea4b1cb0f6 100644 (file)
--- 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 f781f3c97b988934ccf8015f622deffaac4aced5..290966bc7006acfdba5d1bb745a762bac534d7b1 100644 (file)
--- 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