X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcode%2Firrat.lisp;h=f2947a7e0be502925d77addc8ce319a6c7fbe75e;hb=01044af1b8d69fc3899dc0417064c1512223223d;hp=7b18f5a7856b53bc65b7aba06ed5ad31dc15b835;hpb=3c65762b927af861c9c8bc416e4cbac9a14ec0c3;p=sbcl.git diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp index 7b18f5a..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) @@ -81,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 @@ -92,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) @@ -268,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 (/ (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 @@ -435,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))) @@ -451,11 +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. +;;; 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 @@ -595,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 @@ -604,7 +637,9 @@ ;;; 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)) (defun square (x) @@ -648,7 +683,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 +704,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 +730,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 @@ -866,13 +907,18 @@ ;; 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)))