[oct-cvs] Oct commit: oct qd-fun.lisp
rtoy
rtoy at common-lisp.net
Wed Oct 17 03:44:12 UTC 2007
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
More information about the oct-cvs
mailing list