diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp index 5944182ca66af76571a96269e72e874d22077458..d1e9a26add398f9b4fa44719dc3e30c1ee17bb49 100644 --- a/src/code/irrat.lisp +++ b/src/code/irrat.lisp @@ -309,7 +309,49 @@ #+(or ppc sse2) (frob %%sin %%cos %%tan)) - +;; Linux has a sincos function in the C library. Use it. But we need +;; to do pi reduction ourselves because the C library doesn't do +;; accurate reduction. +#+(or (and linux x86)) +(progn +(declaim (inline %%sincos)) +(export '%%sincos) +(alien:def-alien-routine ("sincos" %%sincos) c-call:void + (x double-float) + (sin double-float :out) + (cos double-float :out)) + +(defun %sincos (theta) + (declare (double-float theta)) + ;; Accurately reduce theta. + (multiple-value-bind (n y0 y1) + (%ieee754-rem-pi/2 theta) + (multiple-value-bind (ignore s c) + (%%sincos y0) + (declare (ignore ignore)) + ;; Figure out which quadrant to use, and finish out the + ;; computation using y1. This is done by using a 1st-order + ;; Taylor expansion about y0. + (flet ((sin2 (s c y) + ;; sin(x+y) = sin(x) + cos(x)*y + (+ s (* c y))) + (cos2 (s c y) + ;; cos(x+y) = cos(x) - sin(x)*y + (- c (* s y)))) + (case (logand n 3) + (0 + (values (cos2 s c y1) + (sin2 s c y1))) + (1 + (values (- (sin2 s c y1)) + (cos2 s c y1))) + (2 + (values (- (cos2 s c y1)) + (- (sin2 s c y1)))) + (3 + (values (sin2 s c y1) + (- (cos2 s c y1))))))))) +) ;;;; Power functions. @@ -964,7 +1006,24 @@ "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)." (if (complexp theta) (error (intl:gettext "Argument to CIS is complex: ~S") theta) - (complex (cos theta) (sin theta)))) + #-(or (and linux x86)) + (complex (cos theta) (sin theta)) + #+(or (and linux x86)) + (number-dispatch ((theta real)) + ((rational) + (let ((arg (coerce theta 'double-float))) + (multiple-value-bind (s c) + (%sincos arg) + (complex (coerce s 'single-float) + (coerce c 'single-float))))) + (((foreach single-float double-float)) + (multiple-value-bind (s c) + (%sincos (coerce theta 'double-float)) + (complex (coerce s '(dispatch-type theta)) + (coerce c '(dispatch-type theta))))) + #+double-double + ((double-double-float) + (complex (cos theta) (sin theta)))))) (defun asin (number) "Return the arc sine of NUMBER." diff --git a/src/compiler/float-tran.lisp b/src/compiler/float-tran.lisp index 5fbc2e03b79247dde76990d2fb9e5736d3d7b524..ecca5f126867a5fb62e464ad5856b360dd27c106 100644 --- a/src/compiler/float-tran.lisp +++ b/src/compiler/float-tran.lisp @@ -731,6 +731,30 @@ (deftransform name ((x) '(double-float) rtype :eval-name t :when :both) `(,prim x)))) +#+(or (and linux x86)) +(progn +(defknown (kernel::%sincos) + (double-float) (values double-float double-float) + (movable foldable flushable)) + +(deftransform cis ((x) (single-float) * :when :both) + `(multiple-value-bind (s c) + (kernel::%sincos (coerce x 'double-float)) + (complex (coerce s 'single-float) + (coerce c 'single-float)))) + +(deftransform cis ((x) (double-float) * :when :both) + `(multiple-value-bind (ignore s c) + (kernel::%%sincos x) + (declare (ignore ignore)) + (complex s c))) + +#+double-double +(deftransform cis ((z) (double-double-float) *) + ;; Cis. + '(complex (cos z) (sin z))) +) + ;;; The argument range is limited on the x86 FP trig. functions. A ;;; post-test can detect a failure (and load a suitable result), but ;;; this test is avoided if possible. @@ -1777,6 +1801,7 @@ (deftransform * ((z w) (,real-type (complex ,type)) *) ;; Real * complex '(complex (* z (realpart w)) (* z (imagpart w)))) + #-(or (and linux x86)) (deftransform cis ((z) ((,type)) *) ;; Cis. '(complex (cos z) (sin z)))