diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp
index 270f1dc3a6a8e250319f9b88bb37d23737f00add..4ccf80a8378292180c737e8011e870c3cb894c5f 100644
--- a/src/code/irrat.lisp
+++ b/src/code/irrat.lisp
@@ -187,30 +187,6 @@
     (%sqrt x))
   )
 
-;;; The standard libm routines for sin, cos, and tan on x86 (Linux,
-;;; 32-bit.  64-bit is apparently ok) and ppc are not very accurate
-;;; for large arguments when compared to sparc (and maxima).  This is
-;;; basically caused by the fact that those libraries do not do an
-;;; accurate argument reduction.  The following functions use some
-;;; routines Sun's free fdlibm library to do accurate reduction.  Then
-;;; we call the standard C functions (or vops for x86) on the reduced
-;;; argument.  This produces much more accurate values.
-;;;
-;;; You can test this by computing (cos (scale-float 1d0 120)).  The
-;;; true answer is -0.9258790228548379d0.
-
-#+(or ppc x86)
-(progn
-(declaim (inline %%ieee754-rem-pi/2))
-;; Basic argument reduction routine.  It returns two values: n and y
-;; such that (n + 8*k)*pi/2+y = x where |y|<pi/4 and n indicates in
-;; which octant the arg lies.  Y is actually computed in two parts,
-;; y[0] and y[1] such that the sum is y, for accuracy.
-
-(alien:def-alien-routine ("__ieee754_rem_pio2" %%ieee754-rem-pi/2) c-call:int
-  (x double-float)
-  (y (* double-float)))
-
 ;; Same as above, but instead of needing to pass an array in, the
 ;; output array is broken up into two output values instead.  This is
 ;; easier for the user, and we don't have to wrap calls with
@@ -221,148 +197,438 @@
   (y0 double-float :out)
   (y1 double-float :out))
 
-)
+;; Implement sin/cos/tan in Lisp.  These are based on the routines
+;; from fdlibm.
 
-;; If the C library is accurate, use %trig as the Lisp name.
-#-(or ppc (and sse2 (not darwin)))
-(progn
-(declaim (inline %sin %cos %tan))
-(macrolet ((frob (alien-name lisp-name)
-	     `(alien:def-alien-routine (,alien-name ,lisp-name) double-float
-		(x double-float))))
-  (frob "sin" %sin)
-  (frob "cos" %cos)
-  (frob "tan" %tan))
-)
-
-;; Make %%trig be the C library routines that don't do accurate
-;; reduction.  This is for PPC and for any SSE2 build except on
-;; Darwin. Darwin has accurate C library routines.
-#+(or ppc (and sse2 (not darwin)))
-(progn
-(declaim (inline %%sin %%cos %%tan))
-(macrolet ((frob (alien-name lisp-name)
-	     `(alien:def-alien-routine (,alien-name ,lisp-name) double-float
-		(x double-float))))
-  (frob "sin" %%sin)
-  (frob "cos" %%cos)
-  (frob "tan" %%tan))
-)
-
-;; When the C library is not accurate, define %trig to do accurate
-;; argument reduction and call the appropriate C function on the
-;; reduced arg.  For x87, we can use the x87 FPU trig instructions.
-#+(or ppc (and x86 (not darwin)))
-(macrolet
-    ((frob (sin cos tan)
-       `(progn
-	  ;; In all of the routines below, we just compute the sum of
-	  ;; y0 and y1 and use that as the (reduced) argument for the
-	  ;; trig functions.  This is slightly less accurate than what
-	  ;; fdlibm does, which calls special functions using y0 and
-	  ;; y1 separately, for greater accuracy.  This isn't
-	  ;; implemented, and some spot checks indicate that what we
-	  ;; have here is accurate.
-	  ;;
-	  ;; For x86 with an fsin/fcos/fptan instruction, the pi/4 is
-	  ;; probably too restrictive.
-	  (defun %sin (x)
-	    (declare (double-float x))
-	    (if (< (abs x) (/ pi 4))
-		(,sin x)
-		;; Argument reduction needed
-		(multiple-value-bind (n y0 y1)
-		    (%ieee754-rem-pi/2 x)
-		  (let ((reduced (+ y0 y1)))
-		    (case (logand n 3)
-		      (0 (,sin reduced))
-		      (1 (,cos reduced))
-		      (2 (- (,sin reduced)))
-		      (3 (- (,cos reduced))))))))
-	  (defun %cos (x)
-	    (declare (double-float x))
-	    (if (< (abs x) (/ pi 4))
-		(,cos x)
-		;; Argument reduction needed
-		(multiple-value-bind (n y0 y1)
-		    (%ieee754-rem-pi/2 x)
-		  (let ((reduced (+ y0 y1)))
-		    (case (logand n 3)
-		      (0 (,cos reduced))
-		      (1 (- (,sin reduced)))
-		      (2 (- (,cos reduced)))
-		      (3 (,sin reduced)))))))
-	  (defun %tan (x)
-	    (declare (double-float x))
-	    (if (< (abs x) (/ pi 4))
-		(,tan x)
-		;; Argument reduction needed
-		(multiple-value-bind (n y0 y1)
-		    (%ieee754-rem-pi/2 x)
-		  (let ((reduced (+ y0 y1)))
-		    (if (evenp n)
-			(,tan reduced)
-			(- (/ (,tan reduced)))))))))))
-  ;; Don't want %sin-quick and friends with sse2.
-  #+(and x86 (not sse2))
-  (frob %sin-quick %cos-quick %tan-quick)
-  #+(or ppc sse2)
-  (frob %%sin %%cos %%tan))
-
-;; 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)))
-)
+;; 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
+			  %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.
+;; Input y is the tail of x.
+;; Input iy indicates whether y is 0. (if iy=0, y assume to be 0). 
+;;
+;; Algorithm
+;;	1. Since sin(-x) = -sin(x), we need only to consider positive x. 
+;;	2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0.
+;;	3. sin(x) is approximated by a polynomial of degree 13 on
+;;	   [0,pi/4]
+;;		  	         3            13
+;;	   	sin(x) ~ x + S1*x + ... + S6*x
+;;	   where
+;;	
+;; 	|sin(x)         2     4     6     8     10     12  |     -58
+;; 	|----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x  +S6*x   )| <= 2
+;; 	|  x 					           | 
+;; 
+;;	4. sin(x+y) = sin(x) + sin'(x')*y
+;;		    ~ sin(x) + (1-x*x/2)*y
+;;	   For better accuracy, let 
+;;		     3      2      2      2      2
+;;		r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
+;;	   then                   3    2
+;;		sin(x) = x + (S1*x + (x *(r-y/2)+y))
+
+(declaim (ftype (function (double-float double-float fixnum)
+			  double-float)
+		kernel-sin))
+
+(defun kernel-sin (x y iy)
+  (declare (type (double-float -1d0 1d0) x y)
+	   (fixnum iy)
+	   (optimize (speed 3) (safety 0)))
+  (let ((ix (ldb (byte 31 0) (kernel:double-float-high-bits x))))
+    (when (< ix #x3e400000)
+      ;; |x| < 2^-27
+      ;; Signal inexact if x /= 0
+      (if (zerop (truncate x))
+	  (return-from kernel-sin x)
+	  (return-from kernel-sin x)))
+    (let* ((s1 -1.66666666666666324348d-01) ; #xBFC55555 #x55555549
+	   (s2  8.33333333332248946124d-03) ; #x3F811111 #x1110F8A6
+	   (s3 -1.98412698298579493134d-04) ; #xBF2A01A0 #x19C161D5
+	   (s4  2.75573137070700676789d-06) ; #x3EC71DE3 #x57B1FE7D
+	   (s5 -2.50507602534068634195d-08) ; #xBE5AE5E6 #x8A2B9CEB
+	   (s6  1.58969099521155010221d-10) ; #x3DE5D93A #x5ACFD57C
+	   (z (* x x))
+	   (v (* z x))
+	   (r (+ s2
+		 (* z
+		    (+ s3
+		       (* z
+			  (+ s4
+			     (* z
+				(+ s5
+				   (* z s6))))))))))
+      (if (zerop iy)
+	  (+ x (* v (+ s1 (* z r))))
+	  (- x (- (- (* z (- (* .5 y)
+			     (* v r)))
+		     y)
+		  (* v s1)))))))
+
+;; kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
+;; Input x is assumed to be bounded by ~pi/4 in magnitude.
+;; Input y is the tail of x. 
+;;
+;; Algorithm
+;;	1. Since cos(-x) = cos(x), we need only to consider positive x.
+;;	2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
+;;	3. cos(x) is approximated by a polynomial of degree 14 on
+;;	   [0,pi/4]
+;;		  	                 4            14
+;;	   	cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
+;;	   where the remez error is
+;;	
+;; 	|              2     4     6     8     10    12     14 |     -58
+;; 	|cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  )| <= 2
+;; 	|    					               | 
+;; 
+;; 	               4     6     8     10    12     14 
+;;	4. let r = C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  , then
+;;	       cos(x) = 1 - x*x/2 + r
+;;	   since cos(x+y) ~ cos(x) - sin(x)*y 
+;;			  ~ cos(x) - x*y,
+;;	   a correction term is necessary in cos(x) and hence
+;;		cos(x+y) = 1 - (x*x/2 - (r - x*y))
+;;	   For better accuracy when x > 0.3, let qx = |x|/4 with
+;;	   the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125.
+;;	   Then
+;;		cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)).
+;;	   Note that 1-qx and (x*x/2-qx) is EXACT here, and the
+;;	   magnitude of the latter is at least a quarter of x*x/2,
+;;	   thus, reducing the rounding error in the subtraction.
+(declaim (ftype (function (double-float double-float)
+			  double-float)
+		kernel-cos))
+
+(defun kernel-cos (x y)
+  (declare (type (double-float -1d0 1d0) x y)
+	   (optimize (speed 3) (safety 0)))
+  ;; cos(-x) = cos(x), so we just compute cos(|x|).
+  (let ((ix (ldb (byte 31 0) (kernel:double-float-high-bits x))))
+    ;; cos(x) = 1 when |x| < 2^-27
+    (when (< ix #x3e400000)
+      ;; Signal inexact if x /= 0
+      (if (zerop (truncate x))
+	  (return-from kernel-cos 1d0)
+	  (return-from kernel-cos 1d0)))
+    (let* ((c1  4.16666666666666019037d-02)
+	   (c2 -1.38888888888741095749d-03)
+	   (c3  2.48015872894767294178d-05)
+	   (c4 -2.75573143513906633035d-07)
+	   (c5  2.08757232129817482790d-09)
+	   (c6 -1.13596475577881948265d-11)
+	   (z (* x x))
+	   (r (* z
+		 (+ c1
+		    (* z
+		       (+ c2
+			  (* z
+			     (+ c3
+				(* z
+				   (+ c4
+				      (* z
+					 (+ c5
+					    (* z c6)))))))))))))
+      (cond ((< ix #x3fd33333)
+	     ;; \x| < 0.3
+	     (- 1 (- (* .5 z)
+		     (- (* z r)
+			(* x y)))))
+	    (t
+	     ;; qx = 0.28125 if |x| > 0.78125, else x/4 dropping the
+	     ;; least significant 32 bits.
+	     (let* ((qx (if (> ix #x3fe90000)
+			    0.28125d0
+			    ;; x/4, exactly, and also dropping the
+			    ;; least significant 32 bits of the
+			    ;; fraction.
+			    (make-double-float (- ix #x00200000)
+					       0)))
+		    (hz (- (* 0.5 z) qx))
+		    (a (- 1 qx)))
+	       (- a (- hz (- (* z r)
+			     (* x y))))))))))
+
+(declaim (type (simple-array double-float (*)) tan-coef))
+(defconstant tan-coef
+  (make-array 13 :element-type 'double-float
+	      :initial-contents
+	      '(3.33333333333334091986d-01
+		1.33333333333201242699d-01
+		5.39682539762260521377d-02
+		2.18694882948595424599d-02
+		8.86323982359930005737d-03
+		3.59207910759131235356d-03
+		1.45620945432529025516d-03
+		5.88041240820264096874d-04
+		2.46463134818469906812d-04
+		7.81794442939557092300d-05
+		7.14072491382608190305d-05
+		-1.85586374855275456654d-05
+		2.59073051863633712884d-05)))
+
+;; kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
+;; Input x is assumed to be bounded by ~pi/4 in magnitude.
+;; Input y is the tail of x.
+;; Input k indicates whether tan (if k = 1) or -1/tan (if k = -1) is returned.
+;;
+;; Algorithm
+;;	1. Since tan(-x) = -tan(x), we need only to consider positive x.
+;;	2. if x < 2^-28 (hx<0x3e300000 0), return x with inexact if x!=0.
+;;	3. tan(x) is approximated by a odd polynomial of degree 27 on
+;;	   [0,0.67434]
+;;		  	         3             27
+;;	   	tan(x) ~ x + T1*x + ... + T13*x
+;;	   where
+;;
+;; 	        |tan(x)         2     4            26   |     -59.2
+;; 	        |----- - (1+T1*x +T2*x +.... +T13*x    )| <= 2
+;; 	        |  x 					|
+;;
+;;	   Note: tan(x+y) = tan(x) + tan'(x)*y
+;;		          ~ tan(x) + (1+x*x)*y
+;;	   Therefore, for better accuracy in computing tan(x+y), let
+;;		     3      2      2       2       2
+;;		r = x *(T2+x *(T3+x *(...+x *(T12+x *T13))))
+;;	   then
+;;		 		    3    2
+;;		tan(x+y) = x + (T1*x + (x *(r+y)+y))
+;;
+;;      4. For x in [0.67434,pi/4],  let y = pi/4 - x, then
+;;		tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
+;;		       = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
+(declaim (ftype (function (double-float double-float fixnum)
+			  double-float)
+		kernel-tan))
+
+(defun kernel-tan (x y iy)
+  (declare (type (double-float -1d0 1d0) x y)
+	   (type (member -1 1) iy)
+	   (optimize (speed 3) (safety 0)))
+  (let* ((hx (kernel:double-float-high-bits x))
+	 (ix (logand hx #x7fffffff))
+	 (w 0d0)
+	 (z 0d0)
+	 (v 0d0)
+	 (s 0d0)
+	 (r 0d0))
+    (declare (double-float w z v s r))
+    (when (< ix #x3e300000)
+      ;; |x| < 2^-28
+      (when (zerop (truncate x))
+	(cond ((zerop (logior (logior ix (kernel:double-float-low-bits x))
+			      (+ iy 1)))
+	       ;; x = 0 (because hi and low bits are 0) and iy = -1
+	       ;; (cot)
+	       (return-from kernel-tan (/ (abs x))))
+	      ((= iy 1)
+	       (return-from kernel-tan x))
+	      (t
+	       ;; x /= 0 and iy = -1 (cot)
+	       ;; Compute -1/(x+y) carefully
+	       (let ((a 0d0)
+		     (tt 0d0))
+		 (setf w (+ x y))
+		 (setf z (make-double-float (double-float-high-bits w) 0))
+		 (setf v (- y (- z x)))
+		 (setf a (/ -1 w))
+		 (setf tt (make-double-float (double-float-high-bits a) 0))
+		 (setf s (+ 1 (* tt z)))
+		 (return-from kernel-tan (+ tt
+					    (* a (+ s (* tt v))))))))))
+    (when (>= ix #x3FE59428)
+      ;; |x| > .6744
+      (when (minusp hx)
+	(setf x (- x))
+	(setf y (- y)))
+      ;; The two constants below are such that pi/4 + pi/4_lo is pi/4
+      ;; to twice the accuracy of a double float.
+      ;;
+      ;; z = pi/4-x
+      (setf z (- (make-double-float #x3FE921FB #x54442D18) x))
+      ;; w = pi/4_lo - y.
+      (setf w (- (make-double-float #x3C81A626 #x33145C07) y))
+      (setf x (+ z w))
+      (setf y 0d0))
+    (setf z (* x x))
+    (setf w (* z z))
+    ;; Break x^5*(T[1]+x^2*T[2]+...) into
+    ;; x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
+    ;; x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
+    (setf r (+ (aref tan-coef 1)
+	       (* w
+		  (+ (aref tan-coef 3)
+		     (* w
+			(+ (aref tan-coef 5)
+			   (* w
+			      (+ (aref tan-coef 7)
+				 (* w
+				    (+ (aref tan-coef 9)
+				       (* w (aref tan-coef 11))))))))))))
+    (setf v (* z
+	       (+ (aref tan-coef 2)
+		  (* w
+		     (+ (aref tan-coef 4)
+			(* w
+			   (+ (aref tan-coef 6)
+			      (* w
+				 (+ (aref tan-coef 8)
+				    (* w
+				       (+ (aref tan-coef 10)
+					  (* w (aref tan-coef 12)))))))))))))
+    (setf s (* z x))
+    (setf r (+ y (* z (+ (* s (+ r v))
+			 y))))
+    (incf r (* s (aref tan-coef 0)))
+    (setf w (+ x r))
+    (when (>= ix #x3FE59428)
+      (let ((v (float iy 1d0)))
+	(return-from kernel-tan
+	  (* (- 1 (logand 2 (ash hx -30)))
+	     (- v
+		(* 2
+		   (- x (- (/ (* w w)
+			      (+ w v))
+			   r))))))))
+    (when (= iy 1)
+      (return-from kernel-tan w))
+    ;; Compute 1/w=1/(x+r) carefully
+    (let ((a 0d0)
+	  (tt 0d0))
+      (setf z (kernel:make-double-float (kernel:double-float-high-bits w) 0))
+      (setf v (- r (- z x)))		; z + v = r + x
+      (setf a (/ -1 w))
+      (setf tt (kernel:make-double-float (kernel:double-float-high-bits a) 0))
+      (setf s (+ 1 (* tt z)))
+      (+ tt
+	 (* a
+	    (+ s (* tt v)))))))
+
+;; Return sine function of x.
+;;
+;; kernel function:
+;;	__kernel_sin		... sine function on [-pi/4,pi/4]
+;;	__kernel_cos		... cose function on [-pi/4,pi/4]
+;;	__ieee754_rem_pio2	... argument reduction routine
+;;
+;; Method.
+;;      Let S,C and T denote the sin, cos and tan respectively on 
+;;	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 
+;;	in [-pi/4 , +pi/4], and let n = k mod 4.
+;;	We have
+;;
+;;          n        sin(x)      cos(x)        tan(x)
+;;     ----------------------------------------------------------
+;;	    0	       S	   C		 T
+;;	    1	       C	  -S		-1/T
+;;	    2	      -S	  -C		 T
+;;	    3	      -C	   S		-1/T
+;;     ----------------------------------------------------------
+;;
+;; Special cases:
+;;      Let trig be any of sin, cos, or tan.
+;;      trig(+-INF)  is NaN, with signals;
+;;      trig(NaN)    is that NaN;
+;;
+;; Accuracy:
+;;	TRIG(x) returns trig(x) nearly rounded 
+(defun %sin (x)
+  (declare (double-float x)
+	   (optimize (speed 3)))
+  (let ((ix (ldb (byte 31 0) (kernel:double-float-high-bits x))))
+    (cond
+      ((<= ix #x3fe921fb)
+       ;; |x| < pi/4, approx
+       (kernel-sin x 0d0 0))
+      ((>= ix #x7ff00000)
+       ;; sin(Inf or NaN) is NaN
+       (- x x))
+      (t
+       ;; Argument reduction needed
+       (multiple-value-bind (n y0 y1)
+	   (%ieee754-rem-pi/2 x)
+	 (case (logand n 3)
+	   (0
+	    (kernel-sin y0 y1 1))
+	   (1
+	    (kernel-cos y0 y1))
+	   (2
+	    (- (kernel-sin y0 y1 1)))
+	   (3
+	    (- (kernel-cos y0 y1)))))))))
+
+(defun %cos (x)
+  (declare (double-float x)
+	   (optimize (speed 3)))
+  (let ((ix (ldb (byte 31 0) (kernel:double-float-high-bits x))))
+    (cond
+      ((< ix #x3fe921fb)
+       ;;|x| < pi/4, approx
+       (kernel-cos x 0d0))
+      ((>= ix #x7ff00000)
+       ;; cos(Inf or NaN) is NaN
+       (- x x))
+      (t
+       ;; Argument reduction needed
+       (multiple-value-bind (n y0 y1)
+	   (%ieee754-rem-pi/2 x)
+	 (ecase (logand n 3)
+	   (0
+	    (kernel-cos y0 y1))
+	   (1
+	    (- (kernel-sin y0 y1 1)))
+	   (2
+	    (- (kernel-cos y0 y1)))
+	   (3
+	    (kernel-sin y0 y1 1))))))))
+
+(defun %tan (x)
+  (declare (double-float x)
+	   (optimize (speed 3)))
+  (let ((ix (logand #x7fffffff (kernel:double-float-high-bits x))))
+    (cond ((<= ix #x3fe921fb)
+	   ;; |x| < pi/4
+	   (kernel-tan x 0d0 1))
+	  ((>= ix #x7ff00000)
+	   ;; tan(Inf or Nan) is NaN
+	   (- x x))
+	  (t
+	   (multiple-value-bind (n y0 y1)
+	       (%ieee754-rem-pi/2 x)
+	     (let ((flag (- 1 (ash (logand n 1) 1))))
+	       ;; flag = 1 if n even, -1 if n odd
+	       (kernel-tan y0 y1 flag)))))))
+
+;; Compute sin and cos of x, simultaneously.
+(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))
 
 
 ;;;; Power functions.
@@ -1018,9 +1284,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)))
diff --git a/src/compiler/disassem.lisp b/src/compiler/disassem.lisp
index 762b77b22d072b53a849b884c8ec0a96d2b8df0d..cb04a003b2834396f70b56632ba585428e63988b 100644
--- a/src/compiler/disassem.lisp
+++ b/src/compiler/disassem.lisp
@@ -3232,12 +3232,17 @@
      dstate
      stream)))
 
-(defun disassemble-segments (segments stream dstate)
+(defun disassemble-segments (segments stream dstate &key
+						      (base 16)
+						      (case :downcase)
+						      (radix *print-radix*))
   "Disassemble the machine code instructions in each memory segment in
   SEGMENTS in turn to STREAM."
   (declare (type list segments)
 	   (type stream stream)
-	   (type disassem-state dstate))
+	   (type disassem-state dstate)
+	   (type (integer 2 36) base)
+	   (type (member :upcase :downcase :capitalize) case))
   (unless (null segments)
     (let ((first (car segments))
 	  (last (car (last segments))))
@@ -3257,15 +3262,21 @@
 	;; Initialize these to a sane value, just in case.
 	(setf vm::*note-addis-inst* nil)
 	(setf vm::*pseudo-atomic-set* nil))
-      (dolist (seg segments)
-	(disassemble-segment seg stream dstate)))))
+      (let ((*print-base* base)
+	    (*print-case* case)
+	    (*print-radix* radix))
+	(dolist (seg segments)
+	  (disassemble-segment seg stream dstate))))))
 
 ;;; ----------------------------------------------------------------
 ;;; top-level functions
 
 (defun disassemble-function (function &key (stream *standard-output*)
-				      (use-labels t)
-				      (backend c:*native-backend*))
+					(use-labels t)
+					(backend c:*native-backend*)
+					(base 16)
+					(case :downcase)
+					(radix *print-radix*))
   "Disassemble the machine code instructions for FUNCTION."
   (declare (type compiled-function function)
 	   (type stream stream)
@@ -3275,7 +3286,8 @@
 	 (segments (get-function-segments function)))
     (when use-labels
       (label-segments segments dstate))
-    (disassemble-segments segments stream dstate)))
+    (disassemble-segments segments stream dstate
+			  :base base :case case :radix radix)))
 
 (defun compile-function-lambda-expr (function)
   (declare (type function function))
@@ -3328,33 +3340,49 @@
   (declare (type (or function symbol cons) object)
 	   (type (or (member t) stream) stream)
 	   (type (member t nil) use-labels)
-	   (type c::backend backend)
-	   (type (integer 2 36) base)
-	   (type (member :upcase :downcase :capitalize) case))
-  (let ((*print-base* base)
-	(*print-case* case)
-	(*print-radix* radix)
-	(fun (compiled-function-or-lose object)))
+	   (type c::backend backend))
+  (let ((fun (compiled-function-or-lose object)))
     (if (typep fun 'kernel:byte-function)
 	(c:disassem-byte-fun fun)
 	;; we can't detect closures, so be careful
 	(disassemble-function (fun-self fun)
 			      :stream stream
 			      :use-labels use-labels
-			      :backend backend)))
+			      :backend backend
+			      :base base
+			      :case case
+			      :radix radix)))
   (values))
 
 (defun disassemble-memory (address
 			   length
 			   &key
-			   (stream *standard-output*)
-			   code-component
-			   (use-labels t)
-			   (backend c:*backend*))
-  "Disassembles the given area of memory starting at ADDRESS and LENGTH long.
-  Note that if CODE-COMPONENT is NIL and this memory could move during a GC,
-  you'd better disable it around the call to this function."
-  (declare (type (or address system:system-area-pointer) address)
+			     (stream *standard-output*)
+			     code-component
+			     (use-labels t)
+			     (backend c:*backend*)
+			     (base 16)
+			     (case :downcase)
+			     (radix *print-radix*))
+  "Disassembles the given area of memory starting at ADDRESS and
+  LENGTH (octets) long.  Note that if CODE-COMPONENT is NIL and this
+  memory could move during a GC, you'd better disable it around the
+  call to this function.  ADDRESS can be either an integer or a
+  system-area-pointer.
+
+  :Stream stream
+      The dissassembly is written to this stream.
+  :Use-labels
+      Labels are generated instead of using instruction addresses.
+  :Base
+  :Case
+  :Radix
+      The disassembler uses the specified base, case, and radix when
+      printing the disassembled code.  The default values are 16,
+      :downcase, and *print-radix*, respectively. "
+
+  (declare (type (or address
+		     system:system-area-pointer) address)
 	   (type length length)
 	   (type stream stream)
 	   (type (or null kernel:code-component) code-component)
@@ -3379,12 +3407,16 @@
 	      (list (make-memory-segment address length)))))
     (when use-labels
       (label-segments segments dstate))
-    (disassemble-segments segments stream dstate)))
+    (disassemble-segments segments stream dstate
+			  :base base :case case :radix radix)))
 
 (defun disassemble-code-component (code-component &key
-						  (stream *standard-output*)
-						  (use-labels t)
-						  (backend c:*native-backend*))
+						    (stream *standard-output*)
+						    (use-labels t)
+						    (backend c:*native-backend*)
+						    (base 16)
+						    (case :downcase)
+						    (radix *print-radix*))
   "Disassemble the machine code instructions associated with
   CODE-COMPONENT (this may include multiple entry points)."
   (declare (type (or null kernel:code-component compiled-function)
@@ -3400,7 +3432,8 @@
 	 (segments (get-code-segments code-component)))
     (when use-labels
       (label-segments segments dstate))
-    (disassemble-segments segments stream dstate)))
+    (disassemble-segments segments stream dstate
+			  :base base :case case :radix radix)))
 
 ;;; ----------------------------------------------------------------
 ;;; Code for making useful segments from arbitrary lists of code-blocks
diff --git a/src/compiler/float-tran.lisp b/src/compiler/float-tran.lisp
index a8147d9775382aaeb85d9a95c9ad5e5ba48fb15d..d123d183ea9bf53115425bdc69debb96de1490f8 100644
--- a/src/compiler/float-tran.lisp
+++ b/src/compiler/float-tran.lisp
@@ -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
diff --git a/src/i18n/locale/cmucl-x87.pot b/src/i18n/locale/cmucl-x87.pot
index 2c52230ebde2ba825adf641c78ce330bcafec20f..d0118ced7081b7d26fac2b2ebf4b15a12e8a7478 100644
--- a/src/i18n/locale/cmucl-x87.pot
+++ b/src/i18n/locale/cmucl-x87.pot
@@ -126,10 +126,6 @@ msgstr ""
 msgid "inline scalb function"
 msgstr ""
 
-#: src/compiler/x86/float.lisp
-msgid "inline log1p function"
-msgstr ""
-
 #: src/compiler/x86/float.lisp
 msgid "inline log1p with limited x range function"
 msgstr ""
@@ -146,6 +142,10 @@ msgstr ""
 msgid "inline atan2 function"
 msgstr ""
 
+#: src/compiler/x86/float.lisp
+msgid "inline log1p function"
+msgstr ""
+
 #: src/compiler/x86/float.lisp
 msgid "inline complex single-float creation"
 msgstr ""
diff --git a/src/i18n/locale/cmucl.pot b/src/i18n/locale/cmucl.pot
index 283dcaf3fc13a8899508dd9c50712df7f4ec171f..0c327c9b14909eaa345f0002bc9de83809e56fda 100644
--- a/src/i18n/locale/cmucl.pot
+++ b/src/i18n/locale/cmucl.pot
@@ -17157,10 +17157,22 @@ msgstr ""
 
 #: src/compiler/disassem.lisp
 msgid ""
-"Disassembles the given area of memory starting at ADDRESS and LENGTH long.\n"
-"  Note that if CODE-COMPONENT is NIL and this memory could move during a GC,"
+"Disassembles the given area of memory starting at ADDRESS and\n"
+"  LENGTH (octets) long.  Note that if CODE-COMPONENT is NIL and this\n"
+"  memory could move during a GC, you'd better disable it around the\n"
+"  call to this function.  ADDRESS can be either an integer or a\n"
+"  system-area-pointer.\n"
 "\n"
-"  you'd better disable it around the call to this function."
+"  :Stream stream\n"
+"      The dissassembly is written to this stream.\n"
+"  :Use-labels\n"
+"      Labels are generated instead of using instruction addresses.\n"
+"  :Base\n"
+"  :Case\n"
+"  :Radix\n"
+"      The disassembler uses the specified base, case, and radix when\n"
+"      printing the disassembled code.  The default values are 16,\n"
+"      :downcase, and *print-radix*, respectively. "
 msgstr ""
 
 #: src/compiler/disassem.lisp
diff --git a/src/lisp/Config.sparc_common b/src/lisp/Config.sparc_common
index f8f752474b7c56a754324b6a3de45cb7efac3ee7..c4a8ed6d728eae2ff7cd3d086803852e74a1d4f4 100644
--- a/src/lisp/Config.sparc_common
+++ b/src/lisp/Config.sparc_common
@@ -56,6 +56,6 @@ ASSEM_SRC = sparc-assem.S
 ARCH_SRC = sparc-arch.c
 
 DEPEND=$(CC) 
-OS_SRC = solaris-os.c os-common.c undefineds.c elf.c k_rem_pio2.c
+OS_SRC = solaris-os.c os-common.c undefineds.c elf.c e_rem_pio2.c k_rem_pio2.c
 OS_LIBS= -lsocket -lnsl -ldl
 EXEC_FINAL_OBJ = exec-final.o
diff --git a/src/lisp/Config.x86_common b/src/lisp/Config.x86_common
index 6f0ec952f8df4e665a4cc130ca38daa8d38c1adc..f263d403d62578ba20b7ac76056f812171ca7549 100644
--- a/src/lisp/Config.x86_common
+++ b/src/lisp/Config.x86_common
@@ -65,7 +65,13 @@ NM = nm -gp
 DEPEND_FLAGS = -MM
 
 # This no longer has aliasing problems, so no need to use
-# -ffloat-store and -fno-strict-aliasing anymore.
+# -fno-strict-aliasing anymore.  However, if we're building with x87,
+# we MUST use -ffloat-store to get proper double-float rounding.
 e_rem_pio2.o : e_rem_pio2.c
-	$(CC) -c $(CFLAGS) $(CPPFLAGS) $<
+	$(CC) -c $(CFLAGS) $(CPPFLAGS) $(CC_REM_PIO2) $<
+
+k_rem_pio2.o : k_rem_pio2.c
+	$(CC) -c $(CFLAGS) $(CPPFLAGS) $(CC_REM_PIO2) $<
+
+
 
diff --git a/src/lisp/Config.x86_freebsd b/src/lisp/Config.x86_freebsd
index a050e50f9dd8cb7a2fd9551ed05982b185213db5..911e2c2578f9dd6f91bc1b40652b19ba87968757 100644
--- a/src/lisp/Config.x86_freebsd
+++ b/src/lisp/Config.x86_freebsd
@@ -1,6 +1,10 @@
 # -*- Mode: makefile -*-
 include Config.x86_common
 
+# Need -ffloat-store for e_rem_pio2 and k_rem_pio2 to get properly
+# rounded double-floats while using x87 extended precision.
+CC_REM_PIO2 = -ffloat-store
+
 UNDEFSYMPATTERN = -Xlinker -u -Xlinker &
 OS_SRC += FreeBSD-os.c elf.c
 OS_LINK_FLAGS = -dynamic -export-dynamic
diff --git a/src/lisp/Config.x86_linux b/src/lisp/Config.x86_linux
index 569f4f5f9d2d5c30eb64eef6931dc17be9423928..13eb01219ac10a96bf244428c5ed89b6c46baba9 100644
--- a/src/lisp/Config.x86_linux
+++ b/src/lisp/Config.x86_linux
@@ -3,6 +3,10 @@ include Config.x86_common
 
 CPPFLAGS += -m32 -rdynamic -D__NO_CTYPE -D_GNU_SOURCE
 
+# Need -ffloat-store for e_rem_pio2 and k_rem_pio2 to get properly
+# rounded double-floats while using x87 extended precision.
+CC_REM_PIO2 = -ffloat-store
+
 UNDEFSYMPATTERN = -Xlinker -u -Xlinker &
 ASSEM_SRC +=  linux-stubs.S
 OS_SRC += Linux-os.c elf.c
diff --git a/src/lisp/os-common.c b/src/lisp/os-common.c
index 2ba4bfd8133272a8955b492d31d1c0dddbb8fdd4..90458c9509ed1e1897db74cd9d493ee108a9150f 100755
--- a/src/lisp/os-common.c
+++ b/src/lisp/os-common.c
@@ -542,3 +542,19 @@ os_guard_control_stack(int zone, int guard)
 }
 
 #endif /* not RED_ZONE_HIT */
+
+
+/* Simple interface to __ieee754_rem_pio2 */
+int ieee754_rem_pio2(double x, double *y0, double *y1)
+{
+  extern int __ieee754_rem_pio2(double x, double *y);
+
+  double y[2];
+  int n;
+
+  n = __ieee754_rem_pio2(x, y);
+  *y0 = y[0];
+  *y1 = y[1];
+
+  return n;
+}
diff --git a/src/lisp/ppc-arch.c b/src/lisp/ppc-arch.c
index 6c2ef6b663025e1164d9c8c7a62202118458d352..c46b99c8008f5c5ac21a4a507fe27c7e63d01f2a 100644
--- a/src/lisp/ppc-arch.c
+++ b/src/lisp/ppc-arch.c
@@ -723,17 +723,3 @@ arch_linkage_entry(unsigned long retaddr)
 	/ LinkageEntrySize;
 }
 #endif
-
-int ieee754_rem_pio2(double x, double *y0, double *y1)
-{
-  extern int __ieee754_rem_pio2(double x, double *y);
-
-  double y[2];
-  int n;
-
-  n = __ieee754_rem_pio2(x, y);
-  *y0 = y[0];
-  *y1 = y[1];
-
-  return n;
-}
diff --git a/src/lisp/x86-arch.c b/src/lisp/x86-arch.c
index 56703348fd1a87b3f01edafd7374b52e85224813..c694ceae943f13a3ac15570dac88c443f66b2e96 100644
--- a/src/lisp/x86-arch.c
+++ b/src/lisp/x86-arch.c
@@ -521,19 +521,3 @@ arch_linkage_entry(unsigned long retaddr)
     return ((retaddr - 5) - FOREIGN_LINKAGE_SPACE_START) / LinkageEntrySize;
 }
 #endif /* LINKAGE_TABLE */
-
-int ieee754_rem_pio2(double x, double *y0, double *y1)
-{
-  extern int __ieee754_rem_pio2(double x, double *y);
-
-  double y[2];
-  int n;
-
-  n = __ieee754_rem_pio2(x, y);
-  *y0 = y[0];
-  *y1 = y[1];
-
-  return n;
-}
-
-  
diff --git a/src/tests/trig.lisp b/src/tests/trig.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..58d10b63f29dbe7c8a6dc4170fc7d3f94946470b
--- /dev/null
+++ b/src/tests/trig.lisp
@@ -0,0 +1,217 @@
+;;; Tests for the basic trig functions, now implemented in Lisp.
+
+(defpackage :trig-tests
+  (:use :cl :lisp-unit))
+
+(in-package "TRIG-TESTS")
+
+(define-test sin.signed-zeroes
+  "Test sin for 0d0 and -0d0"
+  (:tag :sin :signed-zeroes)
+  (assert-eql 0d0 (sin 0d0))
+  (assert-eql -0d0 (sin -0d0)))
+
+
+(define-test sin.very-small
+  "Tests sin for the case of |x| < 2^-27, but not 0."
+  (:tag :sin)
+  (assert-eql (scale-float 1d0 -28)
+	      (sin (scale-float 1d0 -28))))
+
+(define-test sin.no-reduction
+  "Test sin for small args without reduction"
+  (:tag :sin)
+  (assert-eql 0.479425538604203d0
+	      (sin .5d0))
+  (assert-eql -0.479425538604203d0
+	      (sin -0.5d0)))
+
+(define-test sin.pi/2
+  "Test for arg near pi/2"
+  (:tag :sin)
+  (assert-eql 1d0 (sin (/ pi 2))))
+
+(define-test sin.arg-reduction
+  "Test for sin with arg reduction"
+  (:tag :sin)
+  ;; Test for argument reduction with n mod 4 = 0
+  (assert-eql -7.07106781186547675943154203316156531867416581156d-1
+	      (sin (* 7/4 pi)))
+    ;; Test for argument reduction with n mod 4 = 1
+  (assert-eql 7.07106781186547329560731709118834541043171055432d-1
+	      (sin (* 9/4 pi)))
+  ;; Test for argument reduction with n mod 4 = 2
+  (assert-eql 7.07106781186548390575743300374993861263439430213d-1
+	      (sin (* 11/4 pi)))
+  ;; Test for argument reduction with n mod 4 = 3
+  (assert-eql -7.07106781186547871002109559079472349116005337743d-1
+	      (sin (* 13/4 pi)))
+  ;; Test for argument reduction, big value
+  (assert-eql 0.377820109360752d0
+	      (sin (scale-float 1d0 120))))
+
+(define-test sin.exceptions
+  "Test sin for exceptional values"
+  (:tag :sin :exceptions)
+  (kernel::with-float-traps-masked ()
+    (assert-error 'floating-point-invalid-operation
+		  (sin ext:double-float-positive-infinity))
+    (assert-error 'floating-point-invalid-operation
+		  (sin ext:double-float-negative-infinity))))
+
+(define-test cos.signed-zeroes
+  "Test cos for 0d0 and -0d0"
+  (:tag :cos :signed-zeroes)
+  (assert-eql 1d0 (cos 0d0))
+  (assert-eql 1d0 (cos -0d0)))
+
+(define-test cos.very-small
+  "Test cos for |x| < 2^-27"
+  (:tag :cos)
+  (assert-eql 1d0 (cos (scale-float 1d0 -28))))
+
+(define-test cos.code-paths
+  "Tests various code paths in cos evaluation"
+  (:tag :cos)
+  ;; Test for branch |x| < .3
+  (assert-eql 0.9689124217106447d0
+	      (cos 0.25d0))
+  ;; Test for branch |x| > .3 and \x| < .78125
+  (assert-eql 8.7758256189037271611628158260382965199164519711d-1
+	      (cos 0.5d0))
+  ;; Test for branch |x| > .3 and |x| > .78125
+  (assert-eql 0.7073882691671998d0
+	      (cos 0.785d0)))
+
+(define-test cos.pi/2  
+  "Test cos(pi/2)"
+  (:tag :cos)
+  (assert-eql 6.123233995736766d-17
+	      (cos (/ pi 2))))
+
+(define-test cos.arg-reduction
+  "Test for cos with arg reduction"
+  (:tag :cos)
+    ;; Test for argument reduction with n mod 4 = 0
+  (assert-eql 7.07106781186547372858534520893509069186435867941d-1
+	      (cos (* 7/4 pi)))
+  ;; Test for argument reduction with n mod 4 = 1
+  (assert-eql 7.0710678118654771924095701509080985020443197242d-1
+	      (cos (* 9/4 pi)))
+  ;; Test for argument reduction with n mod 4 = 2
+  (assert-eql -7.07106781186546658225945423833643190916000739026d-1
+	      (cos (* 11/4 pi)))
+  ;; Test for argument reduction with n mod 4 = 3
+  (assert-eql -7.07106781186547177799579165130055836531929091466d-1
+	      (cos (* 13/4 pi)))
+  ;; Test for argument reduction
+  (assert-eql -0.9258790228548379d0
+	      (cos (scale-float 1d0 120))))
+
+(define-test cos.exceptions
+  "Test cos for exceptional values"
+  (:tag :sin :exceptions)
+  (kernel::with-float-traps-masked ()
+    (assert-error 'floating-point-invalid-operation
+		  (cos ext:double-float-positive-infinity))
+    (assert-error 'floating-point-invalid-operation
+		  (cos ext:double-float-negative-infinity))))
+
+(define-test tan.signed-zeroes
+  "Test tan for 0d0 and -0d0"
+  (:tag :tan :signed-zeroes)
+  (assert-eql 0d0 (tan 0d0))
+  (assert-eql -0d0 (tan -0d0)))
+
+(define-test tan.very-small
+  "Test for tan, |x| < 2^-28"
+  (:tag :tan)
+  (assert-eql (scale-float 1d0 -29)
+	      (tan (scale-float 1d0 -29)))
+  (assert-eql (scale-float -1d0 -29)
+	      (tan (scale-float -1d0 -29))))
+
+(define-test tan.pi/2
+  "Test for tan(pi/2)"
+  (:tag :tan)
+  (assert-eql 1.63312393531953697559677370415289165308640681049d16
+	      (tan (/ pi 2))))
+
+(define-test tan.code-paths
+  "Tests for various code paths in tan"
+  (:tag :tan)
+  ;; |x| < .6744
+  (assert-eql 5.4630248984379051325517946578028538329755172018d-1
+	      (tan 0.5d0))
+  ;; |x = 11/16 = 0.6875 > .6744
+  (assert-eql 8.21141801589894121911423965374711700875371645309d-1
+	      (tan (float 11/16 1d0)))
+  ;; This was found by maxima's testsuite.  A bug in kernel-tan when
+  ;; returning cot(x).
+  (assert-eql 2.0000000000000028604455051971538975562294147582d0
+	      (tan 1.107148717794091d0)))
+
+(define-test tan.arg-reduction
+  "Test for tan with arg reduction"
+  (:tag :tan)
+  ;; Test for argument reduction with n even
+  (assert-eql -1.00000000000000042862637970157370388940976433505d0
+	      (tan (* 7/4 pi)))
+  ;; Test for argument reduction with n odd
+  (assert-eql 9.99999999999999448908940383691222098948324989275d-1
+	      (tan (* 9/4 pi)))
+  (assert-eql -4.08066388841804238545143494525595117765084022768d-1
+	      (tan (scale-float 1d0 120))))
+
+(define-test tan.exceptions
+  "Test tan for exceptional values"
+  (:tag :sin :exceptions)
+  (kernel::with-float-traps-masked ()
+    (assert-error 'floating-point-invalid-operation
+		  (tan ext:double-float-positive-infinity))
+    (assert-error 'floating-point-invalid-operation
+		  (tan ext:double-float-negative-infinity))))
+
+(define-test sincos.signed-zeroes
+  "Test sincos at 0d0, -0d0"
+  (:tag :sincos :signed-zeroes)
+  (assert-equal '(0d0 1d0)
+		(multiple-value-list (kernel::%sincos 0d0)))
+  (assert-equal '(-0d0 1d0)
+		(multiple-value-list (kernel::%sincos -0d0))))
+
+;; Test sincos at a bunch of random points and compare the result from
+;; sin and cos.  If they differ, save the result in a list to be
+;; returned.
+(defun sincos-test (limit n)
+  (let (results)
+    (dotimes (k n)
+      (let* ((x (random limit))
+	     (s-exp (sin x))
+	     (c-exp (cos x)))
+	(multiple-value-bind (s c)
+	    (kernel::%sincos x)
+	  (unless (and (eql s s-exp)
+		       (eql c c-exp))
+	    (push (list x
+			(list s s-exp)
+			(list c c-exp))
+		  results)))))
+    results))
+
+(define-test sincos.consistent
+  "Test sincos is consistent with sin and cos"
+  (:tag :sincos)
+  ;; Small values
+  (assert-eql nil
+	      (sincos-test (/ pi 4) 1000))
+  ;; Medium
+  (assert-eql nil
+	      (sincos-test 16d0 1000))
+  ;; Large
+  (assert-eql nil
+	      (sincos-test (scale-float 1d0 120) 1000))
+  ;; Very large
+  (assert-eql nil
+	      (sincos-test (scale-float 1d0 1023) 1000)))
+