X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcode%2Firrat.lisp;h=7b18f5a7856b53bc65b7aba06ed5ad31dc15b835;hb=aa61c7571b33b86981301f34d3acdb66666f53a3;hp=555c952173267ee1d4157cf78207738540085d1b;hpb=ce62508ec1a0f39008c18a2a5a06461eabe662c0;p=sbcl.git

diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp
index 555c952..7b18f5a 100644
--- a/src/code/irrat.lisp
+++ b/src/code/irrat.lisp
@@ -67,81 +67,6 @@
 #!-x86 (def-math-rtn "sqrt" 1)
 (def-math-rtn "hypot" 2)
 #!-(or hpux x86) (def-math-rtn "log1p" 1)
-
-#!+x86 ;; These are needed for use by byte-compiled files.
-(progn
-  (defun %sin (x)
-    (declare (double-float x)
-	     (values double-float))
-    (%sin x))
-  (defun %sin-quick (x)
-    (declare (double-float x)
-	     (values double-float))
-    (%sin-quick x))
-  (defun %cos (x)
-    (declare (double-float x)
-	     (values double-float))
-    (%cos x))
-  (defun %cos-quick (x)
-    (declare (double-float x)
-	     (values double-float))
-    (%cos-quick x))
-  (defun %tan (x)
-    (declare (double-float x)
-	     (values double-float))
-    (%tan x))
-  (defun %tan-quick (x)
-    (declare (double-float x)
-	     (values double-float))
-    (%tan-quick x))
-  (defun %atan (x)
-    (declare (double-float x)
-	     (values double-float))
-    (%atan x))
-  (defun %atan2 (x y)
-    (declare (double-float x y)
-	     (values double-float))
-    (%atan2 x y))
-  (defun %exp (x)
-    (declare (double-float x)
-	     (values double-float))
-    (%exp x))
-  (defun %log (x)
-    (declare (double-float x)
-	     (values double-float))
-    (%log x))
-  (defun %log10 (x)
-    (declare (double-float x)
-	     (values double-float))
-    (%log10 x))
-  #+nil ;; notyet
-  (defun %pow (x y)
-    (declare (type (double-float 0d0) x)
-	     (double-float y)
-	     (values (double-float 0d0)))
-    (%pow x y))
-  (defun %sqrt (x)
-    (declare (double-float x)
-	     (values double-float))
-    (%sqrt x))
-  (defun %scalbn (f ex)
-    (declare (double-float f)
-	     (type (signed-byte 32) ex)
-	     (values double-float))
-    (%scalbn f ex))
-  (defun %scalb (f ex)
-    (declare (double-float f ex)
-	     (values double-float))
-    (%scalb f ex))
-  (defun %logb (x)
-    (declare (double-float x)
-	     (values double-float))
-    (%logb x))
-  (defun %log1p (x)
-    (declare (double-float x)
-	     (values double-float))
-    (%log1p x))
-  ) ; progn
 
 ;;;; power functions
 
@@ -194,7 +119,7 @@
 ;;; from the general complex case.
 (defun expt (base power)
   #!+sb-doc
-  "Returns BASE raised to the POWER."
+  "Return BASE raised to the POWER."
   (if (zerop power)
       (1+ (* base power))
     (labels (;; determine if the double float is an integer.
@@ -387,7 +312,7 @@
 
 (defun abs (number)
   #!+sb-doc
-  "Returns the absolute value of the number."
+  "Return the absolute value of the number."
   (number-dispatch ((number number))
     (((foreach single-float double-float fixnum rational))
      (abs number))
@@ -458,7 +383,7 @@
 
 (defun cis (theta)
   #!+sb-doc
-  "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
+  "Return cos(Theta) + i sin(Theta), i.e. exp(i Theta)."
   (declare (type real theta))
   (complex (cos theta) (sin theta)))
 
@@ -532,12 +457,6 @@
 ;; error. These bad definitions also mean that sin and cos for
 ;; complex numbers can also lose big.
 
-#+nil
-(defun sinh (number)
-  #!+sb-doc
-  "Return the hyperbolic sine of NUMBER."
-  (/ (- (exp number) (exp (- number))) 2))
-
 (defun sinh (number)
   #!+sb-doc
   "Return the hyperbolic sine of NUMBER."
@@ -549,12 +468,6 @@
        (complex (* (sinh x) (cos y))
 		(* (cosh x) (sin y)))))))
 
-#+nil
-(defun cosh (number)
-  #!+sb-doc
-  "Return the hyperbolic cosine of NUMBER."
-  (/ (+ (exp number) (exp (- number))) 2))
-
 (defun cosh (number)
   #!+sb-doc
   "Return the hyperbolic cosine of NUMBER."
@@ -617,8 +530,11 @@
     ((complex)
      (complex-atanh number))))
 
-;;; HP-UX does not supply a C version of log1p, so
-;;; use the definition.
+;;; HP-UX does not supply a C version of log1p, so use the definition.
+;;; 
+;;; FIXME: This is really not a good definition. As per Raymond Toy
+;;; working on CMU CL, "The definition really loses big-time in
+;;; roundoff as x gets small."
 #!+hpux
 #!-sb-fluid (declaim (inline %log1p))
 #!+hpux
@@ -627,7 +543,7 @@
 	   (optimize (speed 3) (safety 0)))
   (the double-float (log (the (double-float 0d0) (+ number 1d0)))))
 
-;;;; OLD-SPECFUN stuff
+;;;; not-OLD-SPECFUN stuff
 ;;;;
 ;;;; (This was conditional on #-OLD-SPECFUN in the CMU CL sources,
 ;;;; but OLD-SPECFUN was mentioned nowhere else, so it seems to be
@@ -691,10 +607,8 @@
 ;;; should be used instead?
 
 (declaim (inline square))
-(declaim (ftype (function (double-float) (double-float 0d0)) square))
 (defun square (x)
-  (declare (double-float x)
-	   (values (double-float 0d0)))
+  (declare (double-float x))
   (* x x))
 
 ;;; original CMU CL comment, apparently re. SCALB and LOGB and
@@ -711,6 +625,18 @@
 	   (type double-float-exponent n))
   (scale-float x n))
 
+;;; This is like LOGB, but X is not infinity and non-zero and not a
+;;; NaN, so we can always return an integer.
+(declaim (inline logb-finite))
+(defun logb-finite (x)
+  (declare (type double-float x))
+  (multiple-value-bind (signif exponent sign)
+      (decode-float x)
+    (declare (ignore signif sign))
+    ;; DECODE-FLOAT is almost right, except that the exponent is off
+    ;; by one.
+    (1- exponent)))
+
 ;;; Compute an integer N such that 1 <= |2^N * x| < 2.
 ;;; For the special cases, the following values are used:
 ;;;    x             logb
@@ -725,17 +651,11 @@
 	 sb!ext:double-float-positive-infinity)
 	((zerop x)
 	 ;; The answer is negative infinity, but we are supposed to
-	 ;; signal divide-by-zero.
-	 ;; (error 'division-by-zero :operation 'logb :operands (list x))
+          ;; signal divide-by-zero, so do the actual division
 	 (/ -1.0d0 x)
 	 )
 	(t
-	 (multiple-value-bind (signif expon sign)
-	     (decode-float x)
-	   (declare (ignore signif sign))
-	   ;; DECODE-FLOAT is almost right, except that the exponent
-	   ;; is off by one.
-	   (1- expon)))))
+          (logb-finite x))))
 
 ;;; This function is used to create a complex number of the
 ;;; appropriate type:
@@ -751,46 +671,47 @@
   (if (subtypep (type-of (realpart z)) 'double-float)
       (complex x y)
       ;; Convert anything that's not a DOUBLE-FLOAT to a SINGLE-FLOAT.
-      (complex (float x 1.0)
-	       (float y 1.0))))
+      (complex (float x 1f0)
+	       (float y 1f0))))
 
 ;;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
 ;;; result is r + i*k, where k is an integer.
 #!+long-float (eval-when (:compile-toplevel :load-toplevel :execute)
 		(error "needs work for long float support"))
 (defun cssqs (z)
-  ;; Save all FP flags
   (let ((x (float (realpart z) 1d0))
-	(y (float (imagpart z) 1d0))
-	(k 0)
-	(rho 0d0))
-    (declare (double-float x y)
-	     (type (double-float 0d0) rho)
-	     (fixnum k))
+	(y (float (imagpart z) 1d0)))
     ;; Would this be better handled using an exception handler to
     ;; catch the overflow or underflow signal?  For now, we turn all
     ;; traps off and look at the accrued exceptions to see if any
     ;; signal would have been raised.
     (with-float-traps-masked (:underflow :overflow)
-      (setf rho (+ (square x) (square y)))
+      (let ((rho (+ (square x) (square y))))
+       (declare (optimize (speed 3) (space 0)))
       (cond ((and (or (float-nan-p rho)
 		      (float-infinity-p rho))
 		  (or (float-infinity-p (abs x))
 		      (float-infinity-p (abs y))))
-	     (setf rho sb!ext:double-float-positive-infinity))
+              (values sb!ext:double-float-positive-infinity 0))
 	    ((let ((threshold #.(/ least-positive-double-float
 				   double-float-epsilon))
 		   (traps (ldb sb!vm::float-sticky-bits
 			       (sb!vm:floating-point-modes))))
-	       ;; overflow raised or (underflow raised and rho < lambda/eps)
+                ;; Overflow raised or (underflow raised and rho <
+                ;; lambda/eps)
 	       (or (not (zerop (logand sb!vm:float-overflow-trap-bit traps)))
 		   (and (not (zerop (logand sb!vm:float-underflow-trap-bit
 					    traps)))
 			(< rho threshold))))
-	     (setf k (logb (max (abs x) (abs y))))
-	     (setf rho (+ (square (scalb x (- k)))
-			  (square (scalb y (- k))))))))
-    (values rho k)))
+              ;; If we're here, neither x nor y are infinity and at
+              ;; least one is non-zero.. Thus logb returns a nice
+              ;; integer.
+              (let ((k (- (logb-finite (max (abs x) (abs y))))))
+                (values (+ (square (scalb x k))
+                           (square (scalb y k)))
+                        (- k))))
+             (t
+              (values rho 0)))))))
 
 ;;; principal square root of Z
 ;;;
@@ -799,14 +720,18 @@
   (declare (number z))
   (multiple-value-bind (rho k)
       (cssqs z)
-    (declare (type (double-float 0d0) rho)
-	     (fixnum k))
+    (declare (type (or (member 0d0) (double-float 0d0)) rho)
+             (type fixnum k))
     (let ((x (float (realpart z) 1.0d0))
 	  (y (float (imagpart z) 1.0d0))
 	  (eta 0d0)
 	  (nu 0d0))
       (declare (double-float x y eta nu))
 
+      (locally
+         ;; space 0 to get maybe-inline functions inlined.
+         (declare (optimize (speed 3) (space 0)))
+
       (if (not (float-nan-p x))
 	  (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
 
@@ -827,7 +752,7 @@
 	    (when (< x 0d0)
 		  (setf eta (abs nu))
 		  (setf nu (float-sign y rho))))
-      (coerce-to-complex-type eta nu z))))
+       (coerce-to-complex-type eta nu z)))))
     
 ;;; Compute log(2^j*z).
 ;;;
@@ -848,23 +773,21 @@
 	(y (float (imagpart z) 1.0d0)))
     (multiple-value-bind (rho k)
 	(cssqs z)
-      (declare (type (double-float 0d0) rho)
-	       (fixnum k))
+      (declare (optimize (speed 3)))
       (let ((beta (max (abs x) (abs y)))
 	    (theta (min (abs x) (abs y))))
-	(declare (type (double-float 0d0) beta theta))
-	(if (and (zerop k)
+        (coerce-to-complex-type (if (and (zerop k)
 		 (< t0 beta)
 		 (or (<= beta t1)
 		     (< rho t2)))
-	    (setf rho (/ (%log1p (+ (* (- beta 1.0d0)
+                                  (/ (%log1p (+ (* (- beta 1.0d0)
 				       (+ beta 1.0d0))
 				    (* theta theta)))
-			 2d0))
-	    (setf rho (+ (/ (log rho) 2d0)
-			 (* (+ k j) ln2))))
-	(setf theta (atan y x))
-	(coerce-to-complex-type rho theta z)))))
+                                     2d0)
+                                  (+ (/ (log rho) 2d0)
+                                     (* (+ k j) ln2)))
+                                (atan y x)
+                                z)))))
 
 ;;; log of Z = log |Z| + i * arg Z
 ;;;
@@ -881,20 +804,22 @@
 (defun complex-atanh (z)
   (declare (number z))
   (let* (;; constants
-	 (theta #.(/ (sqrt most-positive-double-float) 4.0d0))
-	 (rho #.(/ 4.0d0 (sqrt most-positive-double-float)))
-	 (half-pi #.(/ pi 2.0d0))
+         (theta (/ (sqrt most-positive-double-float) 4.0d0))
+         (rho (/ 4.0d0 (sqrt most-positive-double-float)))
+         (half-pi (/ pi 2.0d0))
 	 (rp (float (realpart z) 1.0d0))
 	 (beta (float-sign rp 1.0d0))
 	 (x (* beta rp))
 	 (y (* beta (- (float (imagpart z) 1.0d0))))
 	 (eta 0.0d0)
 	 (nu 0.0d0))
-    (declare (double-float theta rho half-pi rp beta y eta nu)
-	     (type (double-float 0d0) x))
+    ;; Shouldn't need this declare.
+    (declare (double-float x y))
+    (locally
+       (declare (optimize (speed 3)))
     (cond ((or (> x theta)
 	       (> (abs y) theta))
-	   ;; to avoid overflow...
+            ;; To avoid overflow...
 	   (setf eta (float-sign y half-pi))
 	   ;; nu is real part of 1/(x + iy).  This is x/(x^2+y^2),
 	   ;; which can cause overflow.  Arrange this computation so
@@ -902,7 +827,6 @@
 	   (setf nu (let* ((x-bigger (> x (abs y)))
 			   (r (if x-bigger (/ y x) (/ x y)))
 			   (d (+ 1.0d0 (* r r))))
-		      (declare (double-float r d))
 		      (if x-bigger
 			  (/ (/ x) d)
 			  (/ (/ r y) d)))))
@@ -912,19 +836,14 @@
 	   ;; tanh(176) is 1.0d0 within working precision.
 	   (let ((t1 (+ 4d0 (square y)))
 		 (t2 (+ (abs y) rho)))
-	     (declare (type (double-float 0d0) t1 t2))
-	     #+nil
 	     (setf eta (log (/ (sqrt (sqrt t1)))
 			    (sqrt t2)))
-	     (setf eta (* 0.5d0 (log (the (double-float 0.0d0)
-					  (/ (sqrt t1) t2)))))
 	     (setf nu (* 0.5d0
 			 (float-sign y
 				     (+ half-pi (atan (* 0.5d0 t2))))))))
 	  (t
 	   (let ((t1 (+ (abs y) rho)))
-	     (declare (double-float t1))
-	     ;; normal case using log1p(x) = log(1 + x)
+              ;; Normal case using log1p(x) = log(1 + x)
 	     (setf eta (* 0.25d0
 			  (%log1p (/ (* 4.0d0 x)
 				     (+ (square (- 1.0d0 x))
@@ -936,29 +855,29 @@
 				  (square t1))))))))
     (coerce-to-complex-type (* beta eta)
 			    (- (* beta nu))
-			    z)))
+                             z))))
 
 ;;; Compute tanh z = sinh z / cosh z.
 (defun complex-tanh (z)
   (declare (number z))
   (let ((x (float (realpart z) 1.0d0))
 	(y (float (imagpart z) 1.0d0)))
-    (declare (double-float x y))
+    (locally
+      ;; 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))
-	   (complex (float-sign x)
-		    (float-sign y 0.0d0)))
+              (coerce-to-complex-type (float-sign x)
+                                      (float-sign y) z))
 	  (t
 	   (let* ((tv (%tan y))
 		  (beta (+ 1.0d0 (* tv tv)))
 		  (s (sinh x))
 		  (rho (sqrt (+ 1.0d0 (* s s)))))
-	     (declare (double-float tv s)
-		      (type (double-float 0.0d0) beta rho))
 	     (if (float-infinity-p (abs tv))
 		 (coerce-to-complex-type (/ rho s)
 					 (/ tv)
@@ -967,7 +886,7 @@
 		   (coerce-to-complex-type (/ (* beta rho s)
 					      den)
 					   (/ tv den)
-					   z))))))))
+                                            z)))))))))
 
 ;;; Compute acos z = pi/2 - asin z.
 ;;;