X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcode%2Firrat.lisp;h=25a5f8b5adf720dcbabe6e16b57e676d415bec83;hb=55e6ffb0b21df99a1908f1c0bc00f0baf4322f92;hp=1da4fb7d3e5e18d8c460bedc75adc09010c5f91f;hpb=667ec9d494530079bef28e8589dd0d3274b935ec;p=sbcl.git diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp index 1da4fb7..25a5f8b 100644 --- a/src/code/irrat.lisp +++ b/src/code/irrat.lisp @@ -14,8 +14,9 @@ ;;;; 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,8 +25,8 @@ (sb!xc:defmacro def-math-rtn (name num-args) (let ((function (symbolicate "%" (string-upcase name)))) `(progn - (proclaim '(inline ,function)) - (sb!alien:def-alien-routine (,name ,function) double-float + (declaim (inline ,function)) + (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)) @@ -40,9 +41,22 @@ ) ; EVAL-WHEN -;;;; stubs for the Unix math library +#!+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))) -;;; Please refer to the Unix man pages for details about these routines. +;;;; stubs for the Unix math library +;;;; +;;;; Many of these are unnecessary on the X86 because they're built +;;;; into the FPU. ;;; trigonometric #!-x86 (def-math-rtn "sin" 1) @@ -81,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 @@ -92,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) @@ -121,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 @@ -268,18 +280,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 @@ -302,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) @@ -398,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) @@ -415,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) @@ -435,7 +483,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 +499,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 @@ -506,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) @@ -524,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) @@ -595,7 +643,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 . ;;; FIXME: In SBCL, the floating point infinity constants like ;;; SB!EXT:DOUBLE-FLOAT-POSITIVE-INFINITY aren't available as @@ -604,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) @@ -648,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 @@ -668,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)))) @@ -692,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 @@ -715,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) @@ -758,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 @@ -793,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 @@ -802,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))) @@ -819,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)))))))) @@ -859,19 +934,25 @@ ;;; 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 ;; 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))) @@ -912,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) @@ -925,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) @@ -938,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) @@ -951,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))) @@ -962,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))) @@ -973,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)))