X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcode%2Firrat.lisp;h=5e2c8d0bc5ec37c63b435cc3130be3cf79b64342;hb=dc9fb9111cb1b645aaede0d3ec019c0f78200be0;hp=5f901c910233db1e52550551e29bae529234b101;hpb=1b87bfbe482c555879d2a902e88f9d147ead394e;p=sbcl.git diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp index 5f901c9..5e2c8d0 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))) @@ -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 ;;;;