X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcode%2Firrat.lisp;h=05410a6888464e252f458fda8cd9509d643b6c64;hb=43a526583b7015e6b9945d16e31da72fda1325f5;hp=c4e9ad46387daaf3b3f2b13a16cee2efb1e2493e;hpb=fcd65db754f3a5062fccf136bc633e658e4967b3;p=sbcl.git diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp index c4e9ad4..05410a6 100644 --- a/src/code/irrat.lisp +++ b/src/code/irrat.lisp @@ -23,15 +23,18 @@ (eval-when (:compile-toplevel :execute) (sb!xc:defmacro def-math-rtn (name num-args) - (let ((function (symbolicate "%" (string-upcase name)))) + (let ((function (symbolicate "%" (string-upcase name))) + (args (loop for i below num-args + collect (intern (format nil "ARG~D" i))))) `(progn (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)) - 'double-float) - results))))))) + (defun ,function ,args + (alien-funcall + (extern-alien ,name + (function double-float + ,@(loop repeat num-args + collect 'double-float))) + ,@args))))) (defun handle-reals (function var) `((((foreach fixnum single-float bignum ratio)) @@ -70,20 +73,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))) @@ -107,8 +122,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 @@ -158,10 +181,15 @@ #!+sb-doc "Return BASE raised to the POWER." (if (zerop power) - (let ((result (1+ (* base power)))) - (if (and (floatp result) (float-nan-p result)) - (float 1 result) - result)) + (if (and (zerop base) (floatp power)) + (error 'arguments-out-of-domain-error + :operands (list base power) + :operation 'expt + :references (list '(:ansi-cl :function expt))) + (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 @@ -616,18 +644,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 ;;;; @@ -800,8 +816,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 < @@ -878,10 +899,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) @@ -979,16 +1012,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