X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcode%2Firrat.lisp;h=bab76615a6dd170b40994be6539439d47b54ec36;hb=bb7c5beef3a2c45f0ff99f8038409dc4787aa295;hp=1da4fb7d3e5e18d8c460bedc75adc09010c5f91f;hpb=667ec9d494530079bef28e8589dd0d3274b935ec;p=sbcl.git

diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp
index 1da4fb7..bab7661 100644
--- a/src/code/irrat.lisp
+++ b/src/code/irrat.lisp
@@ -25,7 +25,7 @@
   (let ((function (symbolicate "%" (string-upcase name))))
     `(progn
        (proclaim '(inline ,function))
-       (sb!alien:def-alien-routine (,name ,function) double-float
+       (sb!alien:define-alien-routine (,name ,function) double-float
          ,@(let ((results nil))
              (dotimes (i num-args (nreverse results))
                (push (list (intern (format nil "ARG-~D" i))
@@ -41,8 +41,9 @@
 ) ; EVAL-WHEN
 
 ;;;; stubs for the Unix math library
-
-;;; Please refer to the Unix man pages for details about these routines.
+;;;;
+;;;; Many of these are unnecessary on the X86 because they're built
+;;;; into the FPU.
 
 ;;; trigonometric
 #!-x86 (def-math-rtn "sin" 1)
@@ -268,18 +269,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) base) ; ANSI spec
+	((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
@@ -435,7 +472,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)))
@@ -451,11 +488,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
@@ -595,7 +632,7 @@
 ;;;;
 ;;;; The original CMU CL code requested:
 ;;;;   Please send any bug reports, comments, or improvements to
-;;;;   Raymond Toy at toy@rtp.ericsson.se.
+;;;;   Raymond Toy at <email address deleted during 2002 spam avalanche>.
 
 ;;; FIXME: In SBCL, the floating point infinity constants like
 ;;; SB!EXT:DOUBLE-FLOAT-POSITIVE-INFINITY aren't available as
@@ -869,7 +906,8 @@
 	      #-(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))
+					(log most-positive-double-float))
+				     4d0))
               (coerce-to-complex-type (float-sign x)
                                       (float-sign y) z))
 	  (t