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

diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp
index fe13620..25a5f8b 100644
--- a/src/code/irrat.lisp
+++ b/src/code/irrat.lisp
@@ -25,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))
@@ -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
@@ -83,10 +95,8 @@
 
 ;;; INTEXP -- Handle the rational base, integer power case.
 
-;;; FIXME: As long as the system dies on stack overflow or memory
-;;; exhaustion, it seems reasonable to have this, but its default
-;;; should be NIL, and when it's NIL, anything should be accepted.
-(defparameter *intexp-maximum-exponent* 10000)
+(declaim (type (or integer null) *intexp-maximum-exponent*))
+(defparameter *intexp-maximum-exponent* nil)
 
 ;;; This function precisely calculates base raised to an integral
 ;;; power. It separates the cases by the sign of power, for efficiency
@@ -94,13 +104,10 @@
 ;;; a positive integer. Values of power are calculated as positive
 ;;; integers, and inverted if negative.
 (defun intexp (base power)
-  (when (> (abs power) *intexp-maximum-exponent*)
-    ;; FIXME: should be ordinary error, not CERROR. (Once we set the
-    ;; default for the variable to NIL, the un-continuable error will
-    ;; be less obnoxious.)
-    (cerror "Continue with calculation."
-	    "The absolute value of ~S exceeds ~S."
-	    power '*intexp-maximum-exponent* base power))
+  (when (and *intexp-maximum-exponent*
+	     (> (abs power) *intexp-maximum-exponent*))
+    (error "The absolute value of ~S exceeds ~S."
+	    power '*intexp-maximum-exponent*))
   (cond ((minusp power)
 	 (/ (intexp base (- power))))
 	((eql base 2)
@@ -123,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
@@ -296,7 +306,7 @@
   "Return the logarithm of NUMBER in the base BASE, which defaults to e."
   (if base-p
       (cond
-	((zerop base) base) ; ANSI spec
+	((zerop base) 0f0) ; FIXME: type
 	((and (typep number '(integer (0) *))
 	      (typep base '(integer (0) *)))
 	 (coerce (/ (log2 number) (log2 base)) 'single-float))
@@ -340,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)
@@ -436,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)
@@ -453,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)
@@ -544,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)
@@ -562,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)
@@ -642,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)
@@ -686,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
@@ -706,9 +732,11 @@
 (defun coerce-to-complex-type (x y z)
   (declare (double-float x y)
 	   (number z))
-  (if (subtypep (type-of (realpart z)) 'double-float)
+  (if (typep (realpart z) 'double-float)
       (complex x y)
-      ;; Convert anything that's not a DOUBLE-FLOAT to a SINGLE-FLOAT.
+      ;; Convert anything that's not already a DOUBLE-FLOAT (because
+      ;; the initial argument was a (COMPLEX DOUBLE-FLOAT) and we
+      ;; haven't done anything to lose precision) to a SINGLE-FLOAT.
       (complex (float x 1f0)
 	       (float y 1f0))))
 
@@ -730,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
@@ -753,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)
@@ -796,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
@@ -831,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
@@ -840,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)))
@@ -857,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))))))))
@@ -897,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
@@ -956,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)
@@ -969,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)
@@ -982,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)
@@ -995,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)))
@@ -1006,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)))
@@ -1017,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)))