X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcode%2Firrat.lisp;h=55dbb831a1d2bf5ea1ed8eace0b6125c4f3da7d9;hb=9c71b0ad73bd4597d3130553bfbe70c172fe0501;hp=5477d6d666eec37347512c2040cfff1ea5a8ca6e;hpb=ddaf14f0438b5252d4d5d149f65921795c9a771d;p=sbcl.git diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp index 5477d6d..55dbb83 100644 --- a/src/code/irrat.lisp +++ b/src/code/irrat.lisp @@ -70,20 +70,32 @@ #!-x86 (def-math-rtn "sin" 1) #!-x86 (def-math-rtn "cos" 1) #!-x86 (def-math-rtn "tan" 1) -(def-math-rtn "asin" 1) -(def-math-rtn "acos" 1) #!-x86 (def-math-rtn "atan" 1) #!-x86 (def-math-rtn "atan2" 2) -(def-math-rtn "sinh" 1) -(def-math-rtn "cosh" 1) #!-win32 (progn + (def-math-rtn "acos" 1) + (def-math-rtn "asin" 1) + (def-math-rtn "cosh" 1) + (def-math-rtn "sinh" 1) (def-math-rtn "tanh" 1) (def-math-rtn "asinh" 1) (def-math-rtn "acosh" 1) (def-math-rtn "atanh" 1)) #!+win32 (progn + (declaim (inline %asin)) + (defun %asin (number) + (%atan (/ number (sqrt (- 1 (* number number)))))) + (declaim (inline %acos)) + (defun %acos (number) + (- (/ pi 2) (%asin number))) + (declaim (inline %cosh)) + (defun %cosh (number) + (/ (+ (exp number) (exp (- number))) 2)) + (declaim (inline %sinh)) + (defun %sinh (number) + (/ (- (exp number) (exp (- number))) 2)) (declaim (inline %tanh)) (defun %tanh (number) (/ (%sinh number) (%cosh number))) @@ -95,7 +107,7 @@ (log (+ number (sqrt (- (* number number) 1.0d0))) #.(exp 1.0d0))) (declaim (inline %atanh)) (defun %atanh (number) - (let ((ratio (/ (1+ number) (1- number)))) + (let ((ratio (/ (+ 1 number) (- 1 number)))) ;; Were we effectively zero? (if (= ratio -1.0d0) 0.0d0 @@ -107,8 +119,16 @@ #!-x86 (def-math-rtn "log10" 1) #!-win32(def-math-rtn "pow" 2) #!-(or x86 x86-64) (def-math-rtn "sqrt" 1) -(def-math-rtn "hypot" 2) -#!-(or hpux x86) (def-math-rtn "log1p" 1) +#!-win32 (def-math-rtn "hypot" 2) +#!-x86 (def-math-rtn "log1p" 1) + +#!+win32 +(progn + ;; FIXME: libc hypot "computes the sqrt(x*x+y*y) without undue overflow or underflow" + ;; ...we just do the stupid simple thing. + (declaim (inline %hypot)) + (defun %hypot (x y) + (sqrt (+ (* x x) (* y y))))) ;;;; power functions @@ -211,6 +231,7 @@ (when (zerop (logior y-ihi y-lo)) (return-from real-expt (coerce 1d0 rtype))) ;; +-NaN return x+y + ;; FIXME: Hardcoded qNaN/sNaN values are not portable. (when (or (> x-ihi #x7ff00000) (and (= x-ihi #x7ff00000) (/= x-lo 0)) (> y-ihi #x7ff00000) @@ -334,11 +355,20 @@ "Return the logarithm of NUMBER in the base BASE, which defaults to e." (if base-p (cond - ((zerop base) 0f0) ; FIXME: type + ((zerop base) + (if (or (typep number 'double-float) (typep base 'double-float)) + 0.0d0 + 0.0f0)) ((and (typep number '(integer (0) *)) (typep base '(integer (0) *))) (coerce (/ (log2 number) (log2 base)) 'single-float)) - (t (/ (log number) (log base)))) + ((and (typep number 'integer) (typep base 'double-float)) + ;; No single float intermediate result + (/ (log2 number) (log base 2.0d0))) + ((and (typep number 'double-float) (typep base 'integer)) + (/ (log number 2.0d0) (log2 base))) + (t + (/ (log number) (log base)))) (number-dispatch ((number number)) (((foreach fixnum bignum)) (if (minusp number) @@ -606,18 +636,6 @@ ((complex) (complex-atanh number)))) -;;; 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 -(defun %log1p (number) - (declare (double-float number) - (optimize (speed 3) (safety 0))) - (the double-float (log (the (double-float 0d0) (+ number 1d0))))) ;;;; not-OLD-SPECFUN stuff ;;;; @@ -790,8 +808,13 @@ (values (double-from-bits 0 (1+ sb!vm:double-float-normal-exponent-max) 0) 0)) - ((let ((threshold #.(/ least-positive-double-float - double-float-epsilon)) + ((let ((threshold + ;; (/ least-positive-double-float double-float-epsilon) + (load-time-value + #!-long-float + (sb!kernel:make-double-float #x1fffff #xfffffffe) + #!+long-float + (error "(/ least-positive-long-float long-float-epsilon)"))) (traps (ldb sb!vm::float-sticky-bits (sb!vm:floating-point-modes)))) ;; Overflow raised or (underflow raised and rho < @@ -868,10 +891,22 @@ ;; influences the choices of these constants but doesn't say how to ;; choose them. We'll just assume his choices matches our ;; implementation of log1p. - (let ((t0 #.(/ 1 (sqrt 2.0d0))) + (let ((t0 (load-time-value + #!-long-float + (sb!kernel:make-double-float #x3fe6a09e #x667f3bcd) + #!+long-float + (error "(/ (sqrt 2l0))"))) + ;; KLUDGE: if repeatable fasls start failing under some weird + ;; xc host, this 1.2d0 might be a good place to examine: while + ;; it _should_ be the same in all vaguely-IEEE754 hosts, 1.2 + ;; is not exactly representable, so something could go wrong. (t1 1.2d0) (t2 3d0) - (ln2 #.(log 2d0)) + (ln2 (load-time-value + #!-long-float + (sb!kernel:make-double-float #x3fe62e42 #xfefa39ef) + #!+long-float + (error "(log 2l0)"))) (x (float (realpart z) 1.0d0)) (y (float (imagpart z) 1.0d0))) (multiple-value-bind (rho k) @@ -969,16 +1004,11 @@ ;; space 0 to get maybe-inline functions inlined (declare (optimize (speed 3) (space 0))) (cond ((> (abs x) - ;; 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)) + (load-time-value + #!-long-float + (sb!kernel:make-double-float #x406633ce #x8fb9f87e) + #!+long-float + (error "(/ (+ (log 2l0) (log most-positive-long-float)) 4l0)"))) (coerce-to-complex-type (float-sign x) (float-sign y) z)) (t