diff --git a/qd-fun.lisp b/qd-fun.lisp
index d131a7a642e3fc017d6c8499568943a4f7077350..91c0b329e5c6b8dc40b6f4846fd3a1695f8fff61 100644
--- a/qd-fun.lisp
+++ b/qd-fun.lisp
@@ -696,38 +696,46 @@ is the cosine of A"
 	  (values mod f)
 	  (values (mod (1+ mod) 4) (sub-qd f +qd-pi/2+))))))
 
-(defun rem-pi/2-int (qd)
-  ;; Compute qd rem pi/2 = k*pi/2+y.  So we compute k + y*2/pi =
-  ;; qd*2/pi.
-  ;;
-  ;; First convert qd to 2^e*I.  We already have 2/pi in the form
-  ;; 2^-1584*J.  Then qd*2/pi = 2^(e-1584)*I*J.  Extract out the
-  ;; integer and fractional parts of this.  For the integer part we
-  ;; only need it modulo 4, because of the periodicity.  For the
-  ;; fractional part, we only need 212 (or so bits of fraction).
-  ;;
-  ;; FIXME: But we don't really need to compute all the bits of I*J.
-  ;; In the product, we really only need the 2 bits to the left of the
-  ;; binary point, and then 212 bits to the right.  This doesn't
-  ;; require doing the full multiplication.
-  (multiple-value-bind (i e s)
+(defun split-prod-2/pi (qd)
+  ;; Compute k and f such that k + f = qd*2/pi, where k is an intege
+  ;; and f < 1. Since we do not need more than the the least
+  ;; two bits from k, only the least two bits of k are returned.
+  (multiple-value-bind (i e)
       (integer-decode-qd qd)
+    ;; First convert qd to 2^e*I.  We already have 2/pi in the form
+    ;; 2^-1584*J.  Then qd*2/pi = 2^(e-1584)*I*J.  Extract out the
+    ;; integer and fractional parts of this.  For the integer part we
+    ;; only need it modulo 4, because of the periodicity. For the
+    ;; fractional bits, we return the actual fraction as a ratio.
+    ;;
+    ;; FIXME: We compute the entire product I*J, but we only need the
+    ;; 2 bits to the left of the binary point.  But we need many of
+    ;; the bits to the right of the binary point becaue we have no
+    ;; lower bound on the difference between a quad-double and pi.
     (let* ((exp (- e (integer-length +2/pi-bits+)))
-	   (prod (* (* s i) +2/pi-bits+))
-	   (mod (ldb (byte 2 (- exp)) prod))
-	   ;; A quad-double has about 212 bits, but we add another 53
-	   ;; (5 doubles) for some extra accuracty.
-	   (qd-bits 265)
-	   (frac (ldb (byte qd-bits (- (- exp) qd-bits)) prod))
-	   (f (mul-qd (scale-float-qd (rational-to-qd frac) (- qd-bits))
-		      +qd-pi/2+)))
-      ;; We want the remainder part to be <= pi/4 because the trig
-      ;; functions want that.  So if the fraction is too big, adjust
-      ;; it, and mod value.
-      (if (qd-<= (abs-qd f) +qd-pi/4+)
-	  (values mod f)
-	  (values (mod (1+ mod) 4) (sub-qd f +qd-pi/2+))))))
+	   (binary-point (- exp))
+	   (prod (* i +2/pi-bits+))
+	   (int (ldb (byte 2 binary-point) prod))
+	   (frac (/ (ldb (byte binary-point 0) prod)
+		    (ash 1 binary-point))))
+      (values int frac))))
 
+(defun rem-pi/2-int (qd)
+  "Compute qd rem pi/2 = k + f, where k is an integer and |f| <
+  pi/4. Two values are returned: k mod 4 and f."
+  (if (minusp-qd qd)
+      (multiple-value-bind (k f)
+	  (rem-pi/2-int (neg-qd qd))
+	(values (mod (- k) 4) (neg-qd f)))
+      (multiple-value-bind (k f)
+	  (split-prod-2/pi qd)
+	;; The trig functions want the arg to be less than pi/4.  That
+	;; is f < 1/2.  Thus adjust k and f accordingly.
+	(when (> f 1/2)
+	  (setf k (mod (1+ k) 4))
+	  (setf f (- f 1)))
+	(values k (mul-qd +qd-pi/2+ (rational-to-qd f))))))
+	    
 (defun rem-pi/2-int-b (qd)
   (declare (type %quad-double qd))
   (multiple-value-bind (i e s)
@@ -785,6 +793,13 @@ is the cosine of A"
 		   (sub-qd a (mul-qd-d +qd-pi/2+ (float quot 1d0))))))
 	(t
 	 (rem-pi/2-int a))))
+
+(defun rem-pi/2 (a)
+  ;; If the number is small enough, we don't need to use the full
+  ;; precision algorithm to compute the remainder.  The value of 1024
+  ;; here is rather arbitrary.  We should do an analysis to figure
+  ;; where the breakpoint should be.
+  (rem-pi/2-int a))
       
 
 (defun sin-qd (a)
diff --git a/rt-tests.lisp b/rt-tests.lisp
index caf568fafe2422c82bc6fb8ac94aea347a7f80ef..69f6ee7c7bb54b34dae794a14ab117353c65c859 100644
--- a/rt-tests.lisp
+++ b/rt-tests.lisp
@@ -187,6 +187,18 @@
       (check-accuracy 212 val true))
   nil)
 
+;;; sin
+(rt:deftest oct.sin.pi
+    (not (zerop (sin +pi+)))
+  t)
+
+;;; cos
+(rt:deftest oct.cos.big
+  (let* ((val (cos (scale-float #q1 120)))
+	 (err (abs (- val -0.9258790228548379d0))))
+    (<= err 5.2d-17))
+  t)
+
 ;;; Tests of atan where we know the analytical result
 (rt:deftest oct.atan.1
     (let* ((arg (/ (sqrt #q3)))
@@ -1731,6 +1743,6 @@
   nil)
 
 (rt:deftest erfc
-  (check-accuracy 210 (erfc #q-4)
+  (check-accuracy 198 (erfc #q-4)
 		  #q1.9999999845827420997199811478403265131159514278547464108088316570950057869589732)
   nil)