X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcode%2Firrat.lisp;h=25a5f8b5adf720dcbabe6e16b57e676d415bec83;hb=55e6ffb0b21df99a1908f1c0bc00f0baf4322f92;hp=c22fd9311fee9f5244896123f33ead151bacab85;hpb=223e5d1767a058283731879358279385d871ec50;p=sbcl.git

diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp
index c22fd93..25a5f8b 100644
--- a/src/code/irrat.lisp
+++ b/src/code/irrat.lisp
@@ -41,6 +41,18 @@
 
 ) ; EVAL-WHEN
 
+#!+x86 ;; for constant folding
+(macrolet ((def (name ll)
+	     `(defun ,name ,ll (,name ,@ll))))
+  (def %atan2 (x y))
+  (def %atan (x))
+  (def %tan-quick (x))
+  (def %cos-quick (x))
+  (def %sin-quick (x))
+  (def %sqrt (x))
+  (def %log (x))
+  (def %exp (x)))
+
 ;;;; stubs for the Unix math library
 ;;;;
 ;;;; Many of these are unnecessary on the X86 because they're built
@@ -118,7 +130,10 @@
   #!+sb-doc
   "Return BASE raised to the POWER."
   (if (zerop power)
-      (1+ (* base power))
+      (let ((result (1+ (* base power))))
+	(if (and (floatp result) (float-nan-p result))
+	    (float 1 result)
+	    result))
     (labels (;; determine if the double float is an integer.
 	     ;;  0 - not an integer
 	     ;;  1 - an odd int
@@ -335,7 +350,7 @@
 	 (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
     (((foreach single-float double-float))
      (if (minusp number)
-	 (complex-sqrt number)
+	 (complex-sqrt (complex number))
 	 (coerce (%sqrt (coerce number 'double-float))
 		 '(dispatch-type number))))
      ((complex)
@@ -431,7 +446,7 @@
     (((foreach single-float double-float))
      (if (or (> number (coerce 1 '(dispatch-type number)))
 	     (< number (coerce -1 '(dispatch-type number))))
-	 (complex-asin number)
+	 (complex-asin (complex number))
 	 (coerce (%asin (coerce number 'double-float))
 		 '(dispatch-type number))))
     ((complex)
@@ -448,7 +463,7 @@
     (((foreach single-float double-float))
      (if (or (> number (coerce 1 '(dispatch-type number)))
 	     (< number (coerce -1 '(dispatch-type number))))
-	 (complex-acos number)
+	 (complex-acos (complex number))
 	 (coerce (%acos (coerce number 'double-float))
 		 '(dispatch-type number))))
     ((complex)
@@ -539,7 +554,7 @@
 	 (coerce (%acosh (coerce number 'double-float)) 'single-float)))
     (((foreach single-float double-float))
      (if (< number (coerce 1 '(dispatch-type number)))
-	 (complex-acosh number)
+	 (complex-acosh (complex number))
 	 (coerce (%acosh (coerce number 'double-float))
 		 '(dispatch-type number))))
     ((complex)
@@ -557,7 +572,7 @@
     (((foreach single-float double-float))
      (if (or (> number (coerce 1 '(dispatch-type number)))
 	     (< number (coerce -1 '(dispatch-type number))))
-	 (complex-atanh number)
+	 (complex-atanh (complex number))
 	 (coerce (%atanh (coerce number 'double-float))
 		 '(dispatch-type number))))
     ((complex)
@@ -637,7 +652,22 @@
 ;;; they're effectively implemented as special variable references,
 ;;; and the code below which uses them might be unnecessarily
 ;;; inefficient. Perhaps some sort of MAKE-LOAD-TIME-VALUE hackery
-;;; should be used instead?
+;;; should be used instead?  (KLUDGED 2004-03-08 CSR, by replacing the
+;;; special variable references with (probably equally slow)
+;;; constructors)
+;;;
+;;; FIXME: As of 2004-05, when PFD noted that IMAGPART and COMPLEX
+;;; differ in their interpretations of the real line, IMAGPART was
+;;; patch, which without a certain amount of effort would have altered
+;;; all the branch cut treatment.  Clients of these COMPLEX- routines
+;;; were patched to use explicit COMPLEX, rather than implicitly
+;;; passing in real numbers for treatment with IMAGPART, and these
+;;; COMPLEX- functions altered to require arguments of type COMPLEX;
+;;; however, someone needs to go back to Kahan for the definitive
+;;; answer for treatment of negative real floating point numbers and
+;;; branch cuts.  If adjustment is needed, it is probably the removal
+;;; of explicit calls to COMPLEX in the clients of irrational
+;;; functions.  -- a slightly bitter CSR, 2004-05-16
 
 (declaim (inline square))
 (defun square (x)
@@ -681,7 +711,8 @@
   (cond ((float-nan-p x)
 	 x)
 	((float-infinity-p x)
-	 sb!ext:double-float-positive-infinity)
+	 ;; DOUBLE-FLOAT-POSITIVE-INFINITY
+	 (double-from-bits 0 (1+ sb!vm:double-float-normal-exponent-max) 0))
 	((zerop x)
 	 ;; The answer is negative infinity, but we are supposed to
           ;; signal divide-by-zero, so do the actual division
@@ -727,7 +758,10 @@
 		      (float-infinity-p rho))
 		  (or (float-infinity-p (abs x))
 		      (float-infinity-p (abs y))))
-              (values sb!ext:double-float-positive-infinity 0))
+	     ;; DOUBLE-FLOAT-POSITIVE-INFINITY
+	     (values
+	      (double-from-bits 0 (1+ sb!vm:double-float-normal-exponent-max) 0)
+	      0))
 	    ((let ((threshold #.(/ least-positive-double-float
 				   double-float-epsilon))
 		   (traps (ldb sb!vm::float-sticky-bits
@@ -750,9 +784,15 @@
 
 ;;; principal square root of Z
 ;;;
-;;; Z may be any NUMBER, but the result is always a COMPLEX.
+;;; Z may be RATIONAL or COMPLEX; the result is always a COMPLEX.
 (defun complex-sqrt (z)
-  (declare (number z))
+  ;; KLUDGE: Here and below, we can't just declare Z to be of type
+  ;; COMPLEX, because one-arg COMPLEX on rationals returns a rational.
+  ;; Since there isn't a rational negative zero, this is OK from the
+  ;; point of view of getting the right answer in the face of branch
+  ;; cuts, but declarations of the form (OR RATIONAL COMPLEX) are
+  ;; still ugly.  -- CSR, 2004-05-16
+  (declare (type (or complex rational) z))
   (multiple-value-bind (rho k)
       (cssqs z)
     (declare (type (or (member 0d0) (double-float 0d0)) rho)
@@ -793,7 +833,7 @@
 ;;;
 ;;; This is for use with J /= 0 only when |z| is huge.
 (defun complex-log-scaled (z j)
-  (declare (number z)
+  (declare (type (or rational complex) z)
 	   (fixnum j))
   ;; The constants t0, t1, t2 should be evaluated to machine
   ;; precision.  In addition, Kahan says the accuracy of log1p
@@ -828,7 +868,7 @@
 ;;;
 ;;; Z may be any number, but the result is always a complex.
 (defun complex-log (z)
-  (declare (number z))
+  (declare (type (or rational complex) z))
   (complex-log-scaled z 0))
 	       
 ;;; KLUDGE: Let us note the following "strange" behavior. atanh 1.0d0
@@ -837,7 +877,7 @@
 ;;; i*y is never 0 since we have positive and negative zeroes. -- rtoy
 ;;; Compute atanh z = (log(1+z) - log(1-z))/2.
 (defun complex-atanh (z)
-  (declare (number z))
+  (declare (type (or rational complex) z))
   (let* (;; constants
          (theta (/ (sqrt most-positive-double-float) 4.0d0))
          (rho (/ 4.0d0 (sqrt most-positive-double-float)))
@@ -854,25 +894,25 @@
        (declare (optimize (speed 3)))
     (cond ((or (> x theta)
 	       (> (abs y) theta))
-            ;; 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),
+	   ;; To avoid overflow...
+	   (setf nu (float-sign y half-pi))
+	   ;; ETA is real part of 1/(x + iy).  This is x/(x^2+y^2),
 	   ;; which can cause overflow.  Arrange this computation so
 	   ;; that it won't overflow.
-	   (setf nu (let* ((x-bigger (> x (abs y)))
-			   (r (if x-bigger (/ y x) (/ x y)))
-			   (d (+ 1.0d0 (* r r))))
-		      (if x-bigger
-			  (/ (/ x) d)
-			  (/ (/ r y) d)))))
+	   (setf eta (let* ((x-bigger (> x (abs y)))
+			    (r (if x-bigger (/ y x) (/ x y)))
+			    (d (+ 1.0d0 (* r r))))
+		       (if x-bigger
+			   (/ (/ x) d)
+			   (/ (/ r y) d)))))
 	  ((= x 1.0d0)
 	   ;; Should this be changed so that if y is zero, eta is set
 	   ;; to +infinity instead of approx 176?  In any case
 	   ;; tanh(176) is 1.0d0 within working precision.
 	   (let ((t1 (+ 4d0 (square y)))
 		 (t2 (+ (abs y) rho)))
-	     (setf eta (log (/ (sqrt (sqrt t1)))
-			    (sqrt t2)))
+	     (setf eta (log (/ (sqrt (sqrt t1))
+			       (sqrt t2))))
 	     (setf nu (* 0.5d0
 			 (float-sign y
 				     (+ half-pi (atan (* 0.5d0 t2))))))))
@@ -894,7 +934,7 @@
 
 ;;; Compute tanh z = sinh z / cosh z.
 (defun complex-tanh (z)
-  (declare (number z))
+  (declare (type (or rational complex) z))
   (let ((x (float (realpart z) 1.0d0))
 	(y (float (imagpart z) 1.0d0)))
     (locally
@@ -953,7 +993,7 @@
   ;;
   ;; and these two expressions are equal if and only if arg conj z =
   ;; -arg z, which is clearly true for all z.
-  (declare (number z))
+  (declare (type (or rational complex) z))
   (let ((sqrt-1+z (complex-sqrt (+ 1 z)))
 	(sqrt-1-z (complex-sqrt (- 1 z))))
     (with-float-traps-masked (:divide-by-zero)
@@ -966,7 +1006,7 @@
 ;;;
 ;;; Z may be any NUMBER, but the result is always a COMPLEX.
 (defun complex-acosh (z)
-  (declare (number z))
+  (declare (type (or rational complex) z))
   (let ((sqrt-z-1 (complex-sqrt (- z 1)))
 	(sqrt-z+1 (complex-sqrt (+ z 1))))
     (with-float-traps-masked (:divide-by-zero)
@@ -979,7 +1019,7 @@
 ;;;
 ;;; Z may be any NUMBER, but the result is always a COMPLEX.
 (defun complex-asin (z)
-  (declare (number z))
+  (declare (type (or rational complex) z))
   (let ((sqrt-1-z (complex-sqrt (- 1 z)))
 	(sqrt-1+z (complex-sqrt (+ 1 z))))
     (with-float-traps-masked (:divide-by-zero)
@@ -992,7 +1032,7 @@
 ;;;
 ;;; Z may be any number, but the result is always a complex.
 (defun complex-asinh (z)
-  (declare (number z))
+  (declare (type (or rational complex) z))
   ;; asinh z = -i * asin (i*z)
   (let* ((iz (complex (- (imagpart z)) (realpart z)))
 	 (result (complex-asin iz)))
@@ -1003,7 +1043,7 @@
 ;;;
 ;;; Z may be any number, but the result is always a complex.
 (defun complex-atan (z)
-  (declare (number z))
+  (declare (type (or rational complex) z))
   ;; atan z = -i * atanh (i*z)
   (let* ((iz (complex (- (imagpart z)) (realpart z)))
 	 (result (complex-atanh iz)))
@@ -1014,7 +1054,7 @@
 ;;;
 ;;; Z may be any number, but the result is always a complex.
 (defun complex-tan (z)
-  (declare (number z))
+  (declare (type (or rational complex) z))
   ;; tan z = -i * tanh(i*z)
   (let* ((iz (complex (- (imagpart z)) (realpart z)))
 	 (result (complex-tanh iz)))