0.8.1.20:

[sbcl.git] / src / code / irrat.lisp

diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp
index ac8b805..771795a 100644 (file)
--- a/src/code/irrat.lisp
+++ b/src/code/irrat.lisp
@@ -14,8 +14,9 @@
 \f
 ;;;; miscellaneous constants, utility functions, and macros
 
-(defconstant pi 3.14159265358979323846264338327950288419716939937511L0)
-;(defconstant e 2.71828182845904523536028747135266249775724709369996L0)
+(defconstant pi
+  #!+long-float 3.14159265358979323846264338327950288419716939937511l0
+  #!-long-float 3.14159265358979323846264338327950288419716939937511d0)
 
 ;;; Make these INLINE, since the call to C is at least as compact as a
 ;;; Lisp call, and saves number consing to boot.
@@ -24,7 +25,7 @@
 (sb!xc:defmacro def-math-rtn (name num-args)
   (let ((function (symbolicate "%" (string-upcase name))))
     `(progn
-       (proclaim '(inline ,function))
+       (declaim (inline ,function))
        (sb!alien:define-alien-routine (,name ,function) double-float
          ,@(let ((results nil))
              (dotimes (i num-args (nreverse results))
@@ -269,18 +270,54 @@
             (* base power)
             (exp (* power (log base)))))))))
 
+;;; FIXME: Maybe rename this so that it's clearer that it only works
+;;; on integers?
+(defun log2 (x)
+  (declare (type integer x))
+  ;; CMUCL comment:
+  ;;
+  ;;   Write x = 2^n*f where 1/2 < f <= 1.  Then log2(x) = n +
+  ;;   log2(f).  So we grab the top few bits of x and scale that
+  ;;   appropriately, take the log of it and add it to n.
+  ;;
+  ;; Motivated by an attempt to get LOG to work better on bignums.
+  (let ((n (integer-length x)))
+    (if (< n sb!vm:double-float-digits)
+       (log (coerce x 'double-float) 2.0d0)
+       (let ((f (ldb (byte sb!vm:double-float-digits
+                           (- n sb!vm:double-float-digits))
+                     x)))
+         (+ n (log (scale-float (coerce f 'double-float)
+                                (- sb!vm:double-float-digits))
+                   2.0d0))))))
+
 (defun log (number &optional (base nil base-p))
   #!+sb-doc
   "Return the logarithm of NUMBER in the base BASE, which defaults to e."
   (if base-p
-      (if (zerop base)
-         base                          ; ANSI spec
-         (/ (log number) (log base)))
+      (cond
+       ((zerop base) 0f0) ; FIXME: type
+       ((and (typep number '(integer (0) *))
+             (typep base '(integer (0) *)))
+        (coerce (/ (log2 number) (log2 base)) 'single-float))
+       (t (/ (log number) (log base))))
       (number-dispatch ((number number))
-       (((foreach fixnum bignum ratio))
+       (((foreach fixnum bignum))
+        (if (minusp number)
+            (complex (log (- number)) (coerce pi 'single-float))
+            (coerce (/ (log2 number) (log (exp 1.0d0) 2.0d0)) 'single-float)))
+       ((ratio)
         (if (minusp number)
             (complex (log (- number)) (coerce pi 'single-float))
-            (coerce (%log (coerce number 'double-float)) 'single-float)))
+            (let ((numerator (numerator number))
+                  (denominator (denominator number)))
+              (if (= (integer-length numerator)
+                     (integer-length denominator))
+                  (coerce (%log1p (coerce (- number 1) 'double-float))
+                          'single-float)
+                  (coerce (/ (- (log2 numerator) (log2 denominator))
+                             (log (exp 1.0d0) 2.0d0))
+                          'single-float)))))
        (((foreach single-float double-float))
         ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
         ;; Since this doesn't seem to be an implementation issue
@@ -436,7 +473,7 @@
                           (float-sign y pi))
                       (float-sign y (/ pi 2)))
                   (%atan2 y x))))
-       (number-dispatch ((y number) (x number))
+       (number-dispatch ((y real) (x real))
          ((double-float
            (foreach double-float single-float fixnum bignum ratio))
           (atan2 y (coerce x 'double-float)))
@@ -452,11 +489,11 @@
        ((complex)
         (complex-atan y)))))
 
-;; It seems that everyone has a C version of sinh, cosh, and
-;; tanh. Let's use these for reals because the original
-;; implementations based on the definitions lose big in round-off
-;; error. These bad definitions also mean that sin and cos for
-;; complex numbers can also lose big.
+;;; It seems that every target system has a C version of sinh, cosh,
+;;; and tanh. Let's use these for reals because the original
+;;; implementations based on the definitions lose big in round-off
+;;; error. These bad definitions also mean that sin and cos for
+;;; complex numbers can also lose big.
 
 (defun sinh (number)
   #!+sb-doc
@@ -867,13 +904,18 @@
       ;; space 0 to get maybe-inline functions inlined
       (declare (optimize (speed 3) (space 0)))
     (cond ((> (abs x)
-             #-(or linux hpux) #.(/ (asinh most-positive-double-float) 4d0)
-             ;; This is more accurate under linux.
-             #+(or linux hpux) #.(/ (+ (log 2.0d0)
-                                       (log most-positive-double-float))
-                                    4d0))
-              (coerce-to-complex-type (float-sign x)
-                                      (float-sign y) z))
+             ;; FIXME: this form is hideously broken wrt
+             ;; cross-compilation portability.  Much else in this
+             ;; file is too, of course, sometimes hidden by
+             ;; constant-folding, but this one in particular clearly
+             ;; depends on host and target
+             ;; MOST-POSITIVE-DOUBLE-FLOATs being equal.  -- CSR,
+             ;; 2003-04-20
+             #.(/ (+ (log 2.0d0)
+                     (log most-positive-double-float))
+                  4d0))
+          (coerce-to-complex-type (float-sign x)
+                                  (float-sign y) z))
          (t
           (let* ((tv (%tan y))
                  (beta (+ 1.0d0 (* tv tv)))