Commit b79c2872 authored by Raymond Toy's avatar Raymond Toy

Implement sincos using the new Lisp trig routines. This can now be

used for all platforms.

 code/irrat.lisp::
 * Implement %SINCOS

 compiler/float-tran.lisp::
 * Update deftransforms for CIS.  %SINCOS can be used on any platform.

 tests/trig.lisp:
 * Add tests to verify %sincos returns exactly the same values as for
   sin and cos.
parent e6a9577f
......@@ -202,7 +202,9 @@
;; Block compile so the trig routines don't cons their args when
;; calling the kernel trig routines.
(declaim (ext:start-block kernel-sin kernel-cos kernel-tan %sin %cos %tan))
(declaim (ext:start-block kernel-sin kernel-cos kernel-tan
%sin %cos %tan
%sincos))
;; kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
;; Input x is assumed to be bounded by ~pi/4 in magnitude.
......@@ -592,63 +594,32 @@
;; flag = 1 if n even, -1 if n odd
(kernel-tan y0 y1 flag)))))))
(defun %sincos (x)
(declare (double-float x)
(optimize (speed 3)))
(cond ((<= (abs x) (/ pi 4))
(values (kernel-sin x 0d0 0)
(kernel-cos x 0d0)))
(t
;; Argument reduction needed
(multiple-value-bind (n y0 y1)
(%ieee754-rem-pi/2 x)
(case (logand n 3)
(0
(values (kernel-sin y0 y1 1)
(kernel-cos y0 y1)))
(1
(values (kernel-cos y0 y1)
(- (kernel-sin y0 y1 1))))
(2
(values (- (kernel-sin y0 y1 1))
(- (kernel-cos y0 y1))))
(3
(values (- (kernel-cos y0 y1))
(kernel-sin y0 y1 1))))))))
(declaim (ext:end-block))
;; Linux and sparc have a sincos function in the C library. Use it.
;; But on linux we need to do pi reduction ourselves because the C
;; library doesn't do accurate reduction. Sparc does accurate pi
;; reduction, so we don't need to do it ourselves.
#+(or (and linux x86) sparc)
(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))
#+(and linux x86)
(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 (sin2 s c y1)
(cos2 s c y1)))
(1
(values (cos2 s c y1)
(- (sin2 s c y1))))
(2
(values (- (sin2 s c y1))
(- (cos2 s c y1))))
(3
(values (- (cos2 s c y1))
(sin2 s c y1))))))))
#+sparc
(declaim (inline %sinccos))
#+sparc
(defun %sincos (theta)
(multiple-value-bind (ignore s c)
(%%sincos theta)
(declare (ignore ignore))
(values s c)))
)
;;;; Power functions.
......@@ -1303,9 +1274,6 @@
"Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
(if (complexp theta)
(error (intl:gettext "Argument to CIS is complex: ~S") theta)
#-(or (and linux x86) sparc)
(complex (cos theta) (sin theta))
#+(or (and linux x86) sparc)
(number-dispatch ((theta real))
((rational)
(let ((arg (coerce theta 'double-float)))
......
......@@ -731,8 +731,6 @@
(deftransform name ((x) '(double-float) rtype :eval-name t :when :both)
`(,prim x))))
#+(or (and linux x86) sparc)
(progn
(defknown (kernel::%sincos)
(double-float) (values double-float double-float)
(movable foldable flushable))
......@@ -752,7 +750,7 @@
(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
......
......@@ -140,4 +140,89 @@
(rt:deftest tan.misc.1
(tan (scale-float 1d0 120))
-4.08066388841804238545143494525595117765084022768d-1)
(rt:deftest sincos.1
(let (results)
(dotimes (k 1000)
(let* ((x (random (/ pi 4)))
(s-exp (sin x))
(c-exp (cos x)))
(multiple-value-bind (s c)
(kernel::%sincos x)
(unless (and (= s s-exp)
(= c c-exp))
(push (list x
(list s s-exp)
(list c c-exp))
results)))))
results)
nil)
(rt:deftest sincos.2
(let (results)
(dotimes (k 1000)
(let* ((x (random 16d0))
(s-exp (sin x))
(c-exp (cos x)))
(multiple-value-bind (s c)
(kernel::%sincos x)
(unless (and (= s s-exp)
(= c c-exp))
(push (list x
(list s s-exp)
(list c c-exp))
results)))))
results)
nil)
(rt:deftest sincos.3
(let (results)
(dotimes (k 1000)
(let* ((x (random (scale-float 1d0 120)))
(s-exp (sin x))
(c-exp (cos x)))
(multiple-value-bind (s c)
(kernel::%sincos x)
(unless (and (= s s-exp)
(= c c-exp))
(push (list x
(list s s-exp)
(list c c-exp))
results)))))
results)
nil)
(rt:deftest sincos.3a
(let (results)
(dotimes (k 1000)
(let* ((x (- (random (scale-float 1d0 120))))
(s-exp (sin x))
(c-exp (cos x)))
(multiple-value-bind (s c)
(kernel::%sincos x)
(unless (and (= s s-exp)
(= c c-exp))
(push (list x
(list s s-exp)
(list c c-exp))
results)))))
results)
nil)
(rt:deftest sincos.4
(let (results)
(dotimes (k 1000)
(let* ((x (random (scale-float 1d0 1023)))
(s-exp (sin x))
(c-exp (cos x)))
(multiple-value-bind (s c)
(kernel::%sincos x)
(unless (and (= s s-exp)
(= c c-exp))
(push (list x
(list s s-exp)
(list c c-exp))
results)))))
results)
nil)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment