[oct-cvs] Oct commit: oct qd-fun.lisp

[rtoy]

Update of /project/oct/cvsroot/oct
In directory clnet:/tmp/cvs-serv23871

Modified Files:
	qd-fun.lisp 
Log Message:
Add another implementation of REM-PI/2-INT.  This extracts just the
part of 2/pi that is needed to compute the desired result instead of
multiplying by all l584 bits of 2/pi.

Not yet used.


--- /project/oct/cvsroot/oct/qd-fun.lisp	2007/10/16 17:09:46	1.88
+++ /project/oct/cvsroot/oct/qd-fun.lisp	2007/10/17 03:44:12	1.89
@@ -682,6 +682,52 @@
 	  (values mod f)
 	  (values (mod (1+ mod) 4) (sub-qd f +qd-pi/2+))))))
 
+(defun rem-pi/2-int-b (qd)
+  (declare (type %quad-double qd))
+  (multiple-value-bind (i e s)
+      (integer-decode-qd qd)
+    ;; This computes qd = k*pi/2+s.  And because of the 2*pi
+    ;; periodicity of trig functions, we really only need mod(k, 4),
+    ;; the least 2 bits of k.
+    ;;
+    ;; Write qd = 2^e*i.  Write the value 2/pi as sum f[k]*2^(-k), k =
+    ;; 0, 1, 2,... and f[k] is 0 or 1.  Then the product, p, is
+    ;;
+    ;; p = sum i*f[k]*2^(e-k).
+    ;;
+    ;; Assume 2^m <= i < 2^(m+1).  Then it's obvious that
+    ;; i*f[k]*2^(e-k) is greater than 2 if e-k+m > 1, or k < e+m-1.
+    ;; Thus we can ignore the bits of 2/pi for k < e+m-1, because the
+    ;; contribution is greater than 2.  For k >= e+m-1, we will have
+    ;; the product that is between 0 and 4, which is the part we want.
+    ;; We just need to select how many fraction bits we want.
+    ;; 4*53=212 seems a little low, and 5*53=265 should be more than
+    ;; enough.  Let's choose 220.
+    ;;
+    ;; This code would be faster if we knew how many bits are in i.
+    ;; You might think it should be 212 (4*53), but no.  Consider #q1
+    ;; + #q1q-100.  This has way more than 212 bits.
+    (let* ((qd-bits 220)
+	   (m (integer-length i))
+	   (frac-bits (+ qd-bits (- m 2)))
+	   (bits (ldb (byte qd-bits
+			    (- (integer-length +2/pi-bits+) e frac-bits))
+		      +2/pi-bits+))
+	   ;; Multiply these together.
+	   (prod (* s i bits))
+	   ;; Extract out the integer and fractional parts
+	   (frac (ldb (byte frac-bits 0) prod))
+	   (mod (ldb (byte 2 frac-bits) prod))
+	   ;;(mod (ash prod (- frac-bits)))
+	   (f (mul-qd (scale-float-qd (rational-to-qd frac) (- frac-bits))
+		      +qd-pi/2+)))
+      ;; If we did this right, (ash prod (- frac-bits)) should be 2
+      ;; bits long at most.
+      (if (<= (abs (qd-0 f)) (/ pi 4))
+	  (values mod f)
+	  (values (mod (1+ mod) 4)
+		  (sub-qd f +qd-pi/2+))))))
+
 (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