X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcode%2Firrat.lisp;h=f2947a7e0be502925d77addc8ce319a6c7fbe75e;hb=01044af1b8d69fc3899dc0417064c1512223223d;hp=555c952173267ee1d4157cf78207738540085d1b;hpb=ce62508ec1a0f39008c18a2a5a06461eabe662c0;p=sbcl.git diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp index 555c952..f2947a7 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)) @@ -41,8 +42,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) @@ -67,81 +69,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 @@ -156,10 +83,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 @@ -167,13 +92,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) @@ -194,7 +116,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. @@ -343,18 +265,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 (%log (coerce number 'double-float)) 'single-float))) + (coerce (/ (log2 number) (log (exp 1.0d0) 2.0d0)) 'single-float))) + ((ratio) + (if (minusp number) + (complex (log (- number)) (coerce pi '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 @@ -387,7 +345,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 +416,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))) @@ -510,7 +468,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))) @@ -526,17 +484,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. - -#+nil -(defun sinh (number) - #!+sb-doc - "Return the hyperbolic sine of NUMBER." - (/ (- (exp number) (exp (- number))) 2)) +;;; 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 @@ -549,12 +501,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 +563,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 +576,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 @@ -679,7 +628,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 @@ -688,13 +637,13 @@ ;;; 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) (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 +660,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 @@ -722,20 +683,15 @@ (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. - ;; (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: @@ -748,49 +704,55 @@ (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. - (complex (float x 1.0) - (float y 1.0)))) + ;; 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)))) ;;; 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)) + ;; 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 (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 +761,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 +793,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 +814,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 +845,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 +868,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 +877,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 +896,34 @@ (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))) + ;; 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))) (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 +932,7 @@ (coerce-to-complex-type (/ (* beta rho s) den) (/ tv den) - z)))))))) + z))))))))) ;;; Compute acos z = pi/2 - asin z. ;;;