(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))
#!-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
+#!-(and win32 x86)
(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
+ (def-math-rtn "asinh" 1)
+ (def-math-rtn "acosh" 1)
+ (def-math-rtn "atanh" 1)))
#!+win32
(progn
- (declaim (inline %tanh))
- (defun %tanh (number)
- (/ (%sinh number) (%cosh number)))
+ #!-x86-64
+ (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))))
(declaim (inline %asinh))
(defun %asinh (number)
(log (+ number (sqrt (+ (* number number) 1.0d0))) #.(exp 1.0d0)))
#!-x86 (def-math-rtn "exp" 1)
#!-x86 (def-math-rtn "log" 1)
#!-x86 (def-math-rtn "log10" 1)
-#!-win32(def-math-rtn "pow" 2)
+#!-(and win32 x86) (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
+ ;; This is written in a peculiar way to avoid overflow. Note that in
+ ;; sqrt(x^2 + y^2), either square or the sum can overflow.
+ ;;
+ ;; Factoring x^2 out of sqrt(x^2 + y^2) gives us the expression
+ ;; |x|sqrt(1 + (y/x)^2), which, assuming |x| >= |y|, can only overflow
+ ;; if |x| is sufficiently large.
+ ;;
+ ;; The ZEROP test suffices (y is non-negative) to guard against
+ ;; divisions by zero: x >= y > 0.
+ (declaim (inline %hypot))
+ (defun %hypot (x y)
+ (declare (type double-float x y))
+ (let ((x (abs x))
+ (y (abs y)))
+ (when (> y x)
+ (rotatef x y))
+ (if (zerop y)
+ x
+ (let ((y/x (/ y x)))
+ (* x (sqrt (1+ (* y/x y/x)))))))))
\f
;;;; power functions
#!+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
(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)
(coerce (* pow (%cos y*pi))
rtype)
(coerce (* pow (%sin y*pi))
- rtype)))))))))))))
- (declare (inline real-expt))
+ rtype))))))))))))
+ (complex-expt (base power)
+ (if (and (zerop base) (plusp (realpart power)))
+ (* base power)
+ (exp (* power (log base))))))
+ (declare (inline real-expt complex-expt))
(number-dispatch ((base number) (power number))
(((foreach fixnum (or bignum ratio) (complex rational)) integer)
(intexp base power))
(real-expt base power 'double-float))
((double-float single-float)
(real-expt base power 'double-float))
- (((foreach (complex rational) (complex float)) rational)
+ ;; Handle (expt <complex> <rational>), except the case dealt with
+ ;; in the first clause above, (expt <(complex rational)> <integer>).
+ (((foreach (complex rational) (complex single-float)
+ (complex double-float))
+ rational)
(* (expt (abs base) power)
(cis (* power (phase base)))))
- (((foreach fixnum (or bignum ratio) single-float double-float)
- complex)
- (if (and (zerop base) (plusp (realpart power)))
- (* base power)
- (exp (* power (log base)))))
- (((foreach (complex float) (complex rational))
+ ;; The next three clauses handle (expt <real> <complex>).
+ (((foreach fixnum (or bignum ratio) single-float)
+ (foreach (complex single-float) (complex rational)))
+ (complex-expt base power))
+ (((foreach fixnum (or bignum ratio) single-float)
+ (complex double-float))
+ (complex-expt (coerce base 'double-float) power))
+ ((double-float complex)
+ (complex-expt base power))
+ ;; The next three clauses handle (expt <complex> <float>) and
+ ;; (expt <complex> <complex>).
+ (((foreach (complex single-float) (complex rational))
+ (foreach (complex single-float) (complex rational) single-float))
+ (complex-expt base power))
+ (((foreach (complex single-float) (complex rational))
+ (foreach (complex double-float) double-float))
+ (complex-expt (coerce base '(complex double-float)) power))
+ (((complex double-float)
(foreach complex double-float single-float))
- (if (and (zerop base) (plusp (realpart power)))
- (* base power)
- (exp (* power (log base)))))))))
+ (complex-expt base power))))))
;;; FIXME: Maybe rename this so that it's clearer that it only works
;;; on integers?
"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)
((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)))))
\f
;;;; not-OLD-SPECFUN stuff
;;;;
(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 <
;; 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)
;; 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
projects / sbcl.git / blobdiff