oct-cvs
Threads by month
- ----- 2025 -----
- July
- June
- May
- April
- March
- February
- January
- ----- 2024 -----
- December
- November
- October
- September
- August
- July
- June
- May
- April
- March
- February
- January
- ----- 2023 -----
- December
- November
- October
- September
- August
- July
- June
- May
- April
- March
- February
- January
- ----- 2022 -----
- December
- November
- October
- September
- August
- July
- June
- May
- April
- March
- February
- January
- ----- 2021 -----
- December
- November
- October
- September
- August
- July
- June
- May
- April
- March
- February
- January
- ----- 2020 -----
- December
- November
- October
- September
- August
- July
- June
- May
- April
- March
- February
- January
- ----- 2019 -----
- December
- November
- October
- September
- August
- July
- June
- May
- April
- March
- February
- January
- ----- 2018 -----
- December
- November
- October
- September
- August
- July
- June
- May
- April
- March
- February
- January
- ----- 2017 -----
- December
- November
- October
- September
- August
- July
- June
- May
- April
- March
- February
- January
- ----- 2016 -----
- December
- November
- October
- September
- August
- July
- June
- May
- April
- March
- February
- January
- ----- 2015 -----
- December
- November
- October
- September
- August
- July
- June
- May
- April
- March
- February
- January
- ----- 2014 -----
- December
- November
- October
- September
- August
- July
- June
- May
- April
- March
- February
- January
- ----- 2013 -----
- December
- November
- October
- September
- August
- July
- June
- May
- April
- March
- February
- January
- ----- 2012 -----
- December
- November
- October
- September
- August
- July
- June
- May
- April
- March
- February
- January
- ----- 2011 -----
- December
- November
- October
- September
- August
- July
- June
- May
- April
- March
- February
- January
- ----- 2010 -----
- December
- November
- October
- September
- August
- July
- June
- May
- April
- March
- February
- January
- ----- 2009 -----
- December
- November
- October
- September
- August
- July
- June
- May
- April
- March
- February
- January
- ----- 2008 -----
- December
- November
- October
- September
- August
- July
- June
- May
- April
- March
- February
- January
- ----- 2007 -----
- December
- November
- October
- September
- August
- 202 discussions
[oct-cvs] Oct commit: oct qd-const.lisp qd-dd.lisp qd-fun.lisp qd-io.lisp qd-rep.lisp qd-test.lisp qd.lisp
by rtoy 15 Oct '07
by rtoy 15 Oct '07
15 Oct '07
Update of /project/oct/cvsroot/oct
In directory clnet:/tmp/cvs-serv23636
Modified Files:
qd-const.lisp qd-dd.lisp qd-fun.lisp qd-io.lisp qd-rep.lisp
qd-test.lisp qd.lisp
Log Message:
o Oops. Fix up a few IN-PACKAGE's for the new package names.
qd-fun.lisp:
o Comment out the old sin/cos routines
o Fix a few mistakes in accurate-sincos-qd
o Rename accurate-sincos-qd to sincos-qd.
--- /project/oct/cvsroot/oct/qd-const.lisp 2007/10/14 18:38:14 1.18
+++ /project/oct/cvsroot/oct/qd-const.lisp 2007/10/15 18:53:43 1.19
@@ -23,7 +23,7 @@
;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
;;;; OTHER DEALINGS IN THE SOFTWARE.
-(in-package #:qdi)
+(in-package #:octi)
(defconstant +qd-zero+
(make-qd-d 0d0))
--- /project/oct/cvsroot/oct/qd-dd.lisp 2007/09/16 05:04:04 1.10
+++ /project/oct/cvsroot/oct/qd-dd.lisp 2007/10/15 18:53:44 1.11
@@ -23,7 +23,7 @@
;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
;;;; OTHER DEALINGS IN THE SOFTWARE.
-(in-package #:qdi)
+(in-package #:octi)
;;; double-double float routines needed for quad-double.
;;;
--- /project/oct/cvsroot/oct/qd-fun.lisp 2007/10/15 18:21:46 1.84
+++ /project/oct/cvsroot/oct/qd-fun.lisp 2007/10/15 18:53:44 1.85
@@ -32,7 +32,7 @@
;;; argument is real and the result is real. Behavior is undefined if
;;; this doesn't hold.
-(in-package #:qdi)
+(in-package #:octi)
(defun logb-finite (x)
"Same as logb but X is not infinity and non-zero and not a NaN, so
@@ -370,7 +370,11 @@
(declare (type %quad-double a b))
(let ((n (nint-qd (div-qd a b))))
(values n (sub-qd a (mul-qd n b)))))
-
+
+;; Old, original routines. These are correct, but they don't handle
+;; large args because drem-qd isn't accurate enough.
+#+(or)
+(progn
(defun sin-qd (a)
"Sin(a)"
(declare (type %quad-double a))
@@ -611,6 +615,7 @@
;; cos(j*pi/2) = -1, sin(j*pi/2) = 0
(values (neg-qd s)
(neg-qd c))))))))))))
+)
;; A more accurate implementation of sin and cos. We use a more
;; accurate argument reduction using 1584 bits of 2/pi. This makes
@@ -848,21 +853,20 @@
;; cos(j*pi/2) = 0, sin(j*pi/2) = -1
s))))))))
-(defun accurate-sincos-qd (a)
+(defun sincos-qd (a)
(declare (type %quad-double a))
(when (zerop-qd a)
- (return-from accurate-sincos-qd
+ (return-from sincos-qd
(values +qd-zero+
+qd-one+)))
(multiple-value-bind (j r)
(rem-pi/2 a)
(multiple-value-bind (k tmp)
- (divrem-qd tmp +qd-pi/1024+)
+ (divrem-qd r +qd-pi/1024+)
(let* ((k (truncate (qd-0 k)))
- (abs-j (abs j))
(abs-k (abs k)))
- (assert (<= abs-j 2))
+ (assert (<= 0 j 3))
(assert (<= abs-k 256))
;; Compute sin(s) and cos(s)
(multiple-value-bind (sin-t cos-t)
@@ -900,7 +904,7 @@
(format t "c = ~/qd::qd-format/~%" c))
;; sin(x) = sin(s+k*pi/1024) * cos(j*pi/2)
;; + cos(s+k*pi/1024) * sin(j*pi/2)
- (cond ((zerop abs-j)
+ (cond ((zerop j)
;; cos(j*pi/2) = 1, sin(j*pi/2) = 0
(values s c))
((= j 1)
@@ -912,7 +916,7 @@
(neg-qd c)))
((= j 3)
;; cos(j*pi/2) = 0, sin(j*pi/2) = -1
- (values (neg-qd c) s))))))))))
+ (values (neg-qd c) s)))))))))
(defun atan2-qd/newton (y x)
(declare (type %quad-double y x)
--- /project/oct/cvsroot/oct/qd-io.lisp 2007/10/15 18:21:46 1.19
+++ /project/oct/cvsroot/oct/qd-io.lisp 2007/10/15 18:53:44 1.20
@@ -23,7 +23,7 @@
;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
;;;; OTHER DEALINGS IN THE SOFTWARE.
-(in-package #:qdi)
+(in-package #:octi)
;; Smallest exponent for a double-float.
(eval-when (:compile-toplevel :load-toplevel :execute)
--- /project/oct/cvsroot/oct/qd-rep.lisp 2007/09/18 11:20:16 1.7
+++ /project/oct/cvsroot/oct/qd-rep.lisp 2007/10/15 18:53:44 1.8
@@ -23,7 +23,7 @@
;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
;;;; OTHER DEALINGS IN THE SOFTWARE.
-(in-package #:qdi)
+(in-package #:octi)
;;; This file contains the actual representation of a %quad-double
;;; number. The only real requirement for a %quad-double number is an
--- /project/oct/cvsroot/oct/qd-test.lisp 2007/10/15 18:21:47 1.20
+++ /project/oct/cvsroot/oct/qd-test.lisp 2007/10/15 18:53:44 1.21
@@ -24,7 +24,7 @@
;;;; OTHER DEALINGS IN THE SOFTWARE.
-(in-package #:qdi)
+(in-package #:octi)
;; Compute to how many bits EST and TRUE are equal. If they are
;; identical, return T.
--- /project/oct/cvsroot/oct/qd.lisp 2007/10/15 15:45:33 1.56
+++ /project/oct/cvsroot/oct/qd.lisp 2007/10/15 18:53:44 1.57
@@ -36,7 +36,7 @@
;;; to support quad-doubles.
-(in-package #:qdi)
+(in-package #:octi)
#+cmu
(eval-when (:compile-toplevel)
1
0
[oct-cvs] Oct commit: oct qd-class.lisp qd-format.lisp qd-fun.lisp qd-io.lisp qd-methods.lisp qd-package.lisp qd-test.lisp rt-tests.lisp tests.lisp timing.lisp
by rtoy 15 Oct '07
by rtoy 15 Oct '07
15 Oct '07
Update of /project/oct/cvsroot/oct
In directory clnet:/tmp/cvs-serv20299
Modified Files:
qd-class.lisp qd-format.lisp qd-fun.lisp qd-io.lisp
qd-methods.lisp qd-package.lisp qd-test.lisp rt-tests.lisp
tests.lisp timing.lisp
Log Message:
o Rename QUAD-DOUBLE-INTERNAL package to OCT-INTERNAL, with nickname
OCTI instead of QDI.
o Rename OCT package to NET.COMMON-LISP.OCT, with a nickname of OCT
o Remove nickname of QD. (Conflicts with other packages dealing with
quad-doubles.)
o Update all uses of QDI: to OCTI:
qd-fun.lisp:
o Add REM-PI/2 to do a simpler computation if the arg is small
enough. Otherwise, use the accurate but expensive rem operation.
o Renamed ACCURATE-SIN-QD to SIN-QD, etc.
o Update SIN-QD etc to use REM-PI/2.
--- /project/oct/cvsroot/oct/qd-class.lisp 2007/10/13 02:14:43 1.27
+++ /project/oct/cvsroot/oct/qd-class.lisp 2007/10/15 18:21:46 1.28
@@ -47,15 +47,15 @@
#-cmu
(defmethod print-object ((qd qd-real) stream)
- (format stream "~/qdi::qd-format/" (qd-value qd)))
+ (format stream "~/octi::qd-format/" (qd-value qd)))
#+cmu
(defun print-qd (q stream)
(declare (type %quad-double q))
(if (or (ext:float-infinity-p (qd-0 q))
(ext:float-nan-p (qd-0 q)))
- (format stream "~/qdi::qd-format/" q)
- (format stream "#q~/qdi::qd-format/" q)))
+ (format stream "~/octi::qd-format/" q)
+ (format stream "#q~/octi::qd-format/" q)))
#+cmu
(defmethod print-object ((qd qd-real) stream)
(print-qd (qd-value qd) stream))
@@ -67,7 +67,7 @@
(make-instance 'qd-real :value (qd-value x)))
(defmethod print-object ((qd qd-complex) stream)
- (format stream "#q(~<#q~/qdi::qd-format/ #q~/qdi::qd-format/~:@>)"
+ (format stream "#q(~<#q~/octi::qd-format/ #q~/octi::qd-format/~:@>)"
(list (qd-real qd)
(qd-imag qd))))
--- /project/oct/cvsroot/oct/qd-format.lisp 2007/10/13 02:14:43 1.7
+++ /project/oct/cvsroot/oct/qd-format.lisp 2007/10/15 18:21:46 1.8
@@ -70,7 +70,7 @@
1 0))
nil)))
(multiple-value-bind (fstr flen lpoint tpoint)
- (qdi::qd-to-string (qd-value num) spaceleft fdig k fmin)
+ (octi::qd-to-string (qd-value num) spaceleft fdig k fmin)
(when (and d (zerop d)) (setq tpoint nil))
(when w
(cl:decf spaceleft flen)
--- /project/oct/cvsroot/oct/qd-fun.lisp 2007/10/15 03:26:38 1.83
+++ /project/oct/cvsroot/oct/qd-fun.lisp 2007/10/15 18:21:46 1.84
@@ -646,7 +646,52 @@
(values mod f)
(values (mod (1+ mod) 4) (sub-qd f +qd-pi/2+))))))
-(defun accurate-sin-qd (a)
+(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)
+ (integer-decode-qd qd)
+ (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 (<= (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
+ ;; here is rather arbitrary. We should do an analysis to figure
+ ;; where the breakpoint should be.
+ (cond ((<= (abs (qd-0 a)) 256)
+ (let ((quot (truncate (qd-0 (nint-qd (div-qd a +qd-pi/2+))))))
+ (values (mod quot 4)
+ (sub-qd a (mul-qd-d +qd-pi/2+ (float quot 1d0))))))
+ (t
+ (rem-pi/2-int a))))
+
+
+(defun sin-qd (a)
"Sin(a)"
(declare (type %quad-double a))
;; To compute sin(x), choose integers a, b so that
@@ -666,10 +711,10 @@
;; cos(s+k*pi/1024) = cos(s)*cos(k*pi/1024)
;; - sin(s)*sin(k*pi/1024)
(when (zerop-qd a)
- (return-from accurate-sin-qd a))
+ (return-from sin-qd a))
(multiple-value-bind (j r)
- (rem-pi/2-int a)
+ (rem-pi/2 a)
(multiple-value-bind (k tmp)
(divrem-qd r +qd-pi/1024+)
(let* ((k (truncate (qd-0 k)))
@@ -723,7 +768,7 @@
;; cos(j*pi/2) = 0, sin(j*pi/2) = -1
(neg-qd c)))))))))
-(defun accurate-cos-qd (a)
+(defun cos-qd (a)
"Cos(a)"
;; Just like sin-qd, but for cos.
(declare (type %quad-double a))
@@ -744,10 +789,10 @@
;; cos(s+k*pi/1024) = cos(s)*cos(k*pi/1024)
;; - sin(s)*sin(k*pi/1024)
(when (zerop-qd a)
- (return-from accurate-cos-qd +qd-one+))
+ (return-from cos-qd +qd-one+))
(multiple-value-bind (j r)
- (rem-pi/2-int a)
+ (rem-pi/2 a)
(multiple-value-bind (k tmp)
(divrem-qd r +qd-pi/1024+)
(let* ((k (truncate (qd-0 k)))
@@ -811,7 +856,7 @@
+qd-one+)))
(multiple-value-bind (j r)
- (rem-pi/2-int a)
+ (rem-pi/2 a)
(multiple-value-bind (k tmp)
(divrem-qd tmp +qd-pi/1024+)
(let* ((k (truncate (qd-0 k)))
@@ -930,9 +975,9 @@
(yy (div-qd y r)))
#+nil
(progn
- (format t "r = ~/qdi::qd-format/~%" r)
- (format t "xx = ~/qdi::qd-format/~%" xx)
- (format t "yy = ~/qdi::qd-format/~%" yy))
+ (format t "r = ~/octi::qd-format/~%" r)
+ (format t "xx = ~/octi::qd-format/~%" xx)
+ (format t "yy = ~/octi::qd-format/~%" yy))
;; Compute double-precision approximation to atan
(let ((z (make-qd-d (atan (qd-0 y) (qd-0 x))))
--- /project/oct/cvsroot/oct/qd-io.lisp 2007/10/10 15:24:07 1.18
+++ /project/oct/cvsroot/oct/qd-io.lisp 2007/10/15 18:21:46 1.19
@@ -376,7 +376,7 @@
(scale-float (float q1 1d0)
(cl:- len (* 2 53))))
#+(or)
- (format t "~/qdi::qd-format/~%" (mul-qd xx yy)))
+ (format t "~/octi::qd-format/~%" (mul-qd xx yy)))
(if (minusp sign)
(neg-qd (mul-qd xx yy))
(mul-qd xx yy))))))))
--- /project/oct/cvsroot/oct/qd-methods.lisp 2007/10/13 02:14:43 1.63
+++ /project/oct/cvsroot/oct/qd-methods.lisp 2007/10/15 18:21:47 1.64
@@ -26,23 +26,23 @@
(in-package #:oct)
(defconstant +pi+
- (make-instance 'qd-real :value qdi:+qd-pi+)
+ (make-instance 'qd-real :value octi:+qd-pi+)
"Quad-double value of pi")
(defconstant +pi/2+
- (make-instance 'qd-real :value qdi:+qd-pi/2+)
+ (make-instance 'qd-real :value octi:+qd-pi/2+)
"Quad-double value of pi/2")
(defconstant +pi/4+
- (make-instance 'qd-real :value qdi:+qd-pi/4+)
+ (make-instance 'qd-real :value octi:+qd-pi/4+)
"Quad-double value of pi/4")
(defconstant +2pi+
- (make-instance 'qd-real :value qdi:+qd-2pi+)
+ (make-instance 'qd-real :value octi:+qd-2pi+)
"Quad-double value of 2*pi")
(defconstant +log2+
- (make-instance 'qd-real :value qdi:+qd-log2+)
+ (make-instance 'qd-real :value octi:+qd-log2+)
"Quad-double value of log(2), natural log of 2")
#+cmu
@@ -57,7 +57,7 @@
(defconstant +most-positive-quad-double-float+
(make-instance 'qd-real
- :value (qdi::%make-qd-d most-positive-double-float
+ :value (octi::%make-qd-d most-positive-double-float
(cl:scale-float most-positive-double-float (cl:* 1 -53))
(cl:scale-float most-positive-double-float (cl:* 2 -53))
(cl:scale-float most-positive-double-float (cl:* 3 -53)))))
@@ -79,13 +79,7 @@
(defmethod make-qd ((x cl:rational))
;; We should do something better than this.
- (let ((top (numerator x))
- (bot (denominator x)))
- (if (= bot 1)
- (make-instance 'qd-real :value (qdi::make-float (signum top) (abs top) 0 0 0))
- (make-instance 'qd-real
- :value (div-qd (qdi::make-float (signum top) (abs top) 0 0 0)
- (qdi::make-float (signum bot) (abs bot) 0 0 0))))))
+ (make-instance 'qd-real :value (rational-to-qd x)))
(defmethod add1 ((a number))
@@ -424,7 +418,7 @@
(declaim (inline qd-cssqs))
(defun qd-cssqs (z)
(multiple-value-bind (rho k)
- (qdi::hypot-aux-qd (qd-value (realpart z))
+ (octi::hypot-aux-qd (qd-value (realpart z))
(qd-value (imagpart z)))
(values (make-instance 'qd-real :value rho)
k)))
@@ -811,7 +805,7 @@
(defmethod random ((x qd-real) &optional (state *random-state*))
(* x (make-instance 'qd-real
- :value (qdi:random-qd state))))
+ :value (octi:random-qd state))))
(defmethod float-digits ((x cl:real))
(cl:float-digits x))
--- /project/oct/cvsroot/oct/qd-package.lisp 2007/10/10 15:24:07 1.42
+++ /project/oct/cvsroot/oct/qd-package.lisp 2007/10/15 18:21:47 1.43
@@ -23,9 +23,9 @@
;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
;;;; OTHER DEALINGS IN THE SOFTWARE.
-(defpackage #:quad-double-internal
+(defpackage #:oct-internal
(:use #:cl #+cmu #:extensions)
- (:nicknames #:qdi)
+ (:nicknames #:octi)
(:export #:%quad-double
#:read-qd
#:add-qd
@@ -102,9 +102,9 @@
#:two-prod
#:two-sqr))
-(defpackage #:quad-double
- (:use #:cl #:quad-double-internal)
- (:nicknames #:oct #:qd)
+(defpackage #:net.common-lisp.oct
+ (:use #:cl #:oct-internal)
+ (:nicknames #:oct)
(:shadow #:+
#:-
#:*
@@ -240,6 +240,7 @@
#:float-digits
#:rational
#:rationalize
+ #:make-qd
)
;; Constants
(:export #:+pi+
--- /project/oct/cvsroot/oct/qd-test.lisp 2007/09/18 03:05:56 1.19
+++ /project/oct/cvsroot/oct/qd-test.lisp 2007/10/15 18:21:47 1.20
@@ -38,8 +38,8 @@
(cl:- (log err 2d0)))))
(defun print-result (est true)
- (format t "est: ~/qdi::qd-format/~%" est)
- (format t "tru: ~/qdi::qd-format/~%" true)
+ (format t "est: ~/octi::qd-format/~%" est)
+ (format t "tru: ~/octi::qd-format/~%" true)
(format t "err: ~A~%" (qd-0 (sub-qd est true)))
(format t "bits: ~,1f~%" (bit-accuracy est true)))
@@ -218,24 +218,24 @@
(let* ((arg (div-qd +qd-one+ (sqrt-qd (make-qd-d 3d0))))
(y (div-qd (funcall fun arg) +qd-pi+))
(true (div-qd +qd-one+ (make-qd-d 6d0))))
- (format t "atan(1/sqrt(3))/pi = ~/qdi::qd-format/~%" y)
- (format t "1/6 = ~/qdi::qd-format/~%" true)
+ (format t "atan(1/sqrt(3))/pi = ~/octi::qd-format/~%" y)
+ (format t "1/6 = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
;; atan(sqrt(3)) = %pi/3
(let* ((arg (sqrt-qd (make-qd-d 3d0)))
(y (div-qd (funcall fun arg) +qd-pi+))
(true (div-qd +qd-one+ (make-qd-d 3d0))))
- (format t "atan(sqrt(3))/pi = ~/qdi::qd-format/~%" y)
- (format t "1/3 = ~/qdi::qd-format/~%" true)
+ (format t "atan(sqrt(3))/pi = ~/octi::qd-format/~%" y)
+ (format t "1/3 = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
;; atan(1) = %pi/4
(let* ((arg +qd-one+)
(y (div-qd (funcall fun arg) +qd-pi+))
(true (div-qd +qd-one+ (make-qd-d 4d0))))
- (format t "atan(1)/pi = ~/qdi::qd-format/~%" y)
- (format t "1/4 = ~/qdi::qd-format/~%" true)
+ (format t "atan(1)/pi = ~/octi::qd-format/~%" y)
+ (format t "1/4 = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true))))
@@ -244,22 +244,22 @@
(let* ((arg (div-qd +qd-pi+ (make-qd-d 6d0)))
(y (sin-qd arg))
(true (make-qd-d 0.5d0)))
- (format t "sin(pi/6) = ~/qdi::qd-format/~%" y)
- (format t "1/2 = ~/qdi::qd-format/~%" true)
+ (format t "sin(pi/6) = ~/octi::qd-format/~%" y)
+ (format t "1/2 = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
(let* ((arg (div-qd +qd-pi+ (make-qd-d 4d0)))
(y (sin-qd arg))
(true (sqrt-qd (make-qd-d 0.5d0))))
- (format t "sin(pi/4) = ~/qdi::qd-format/~%" y)
- (format t "1/sqrt(2) = ~/qdi::qd-format/~%" true)
+ (format t "sin(pi/4) = ~/octi::qd-format/~%" y)
+ (format t "1/sqrt(2) = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
(let* ((arg (div-qd +qd-pi+ (make-qd-d 3d0)))
(y (sin-qd arg))
(true (div-qd (sqrt-qd (make-qd-d 3d0)) (make-qd-d 2d0))))
- (format t "sin(pi/3) = ~/qdi::qd-format/~%" y)
- (format t "sqrt(3)/2 = ~/qdi::qd-format/~%" true)
+ (format t "sin(pi/3) = ~/octi::qd-format/~%" y)
+ (format t "sqrt(3)/2 = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
)
@@ -269,22 +269,22 @@
(let* ((arg (div-qd +qd-pi+ (make-qd-d 6d0)))
(y (funcall f arg))
(true (div-qd +qd-one+ (sqrt-qd (make-qd-d 3d0)))))
- (format t "tan(pi/6) = ~/qdi::qd-format/~%" y)
- (format t "1/sqrt(3) = ~/qdi::qd-format/~%" true)
+ (format t "tan(pi/6) = ~/octi::qd-format/~%" y)
+ (format t "1/sqrt(3) = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
(let* ((arg (div-qd +qd-pi+ (make-qd-d 4d0)))
(y (funcall f arg))
(true +qd-one+))
- (format t "tan(pi/4) = ~/qdi::qd-format/~%" y)
- (format t "1 = ~/qdi::qd-format/~%" true)
+ (format t "tan(pi/4) = ~/octi::qd-format/~%" y)
+ (format t "1 = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
(let* ((arg (div-qd +qd-pi+ (make-qd-d 3d0)))
(y (funcall f arg))
(true (sqrt-qd (make-qd-d 3d0))))
- (format t "tan(pi/3) = ~/qdi::qd-format/~%" y)
- (format t "sqrt(3) = ~/qdi::qd-format/~%" true)
+ (format t "tan(pi/3) = ~/octi::qd-format/~%" y)
+ (format t "sqrt(3) = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
)
@@ -294,22 +294,22 @@
(let* ((arg (make-qd-d 0.5d0))
(y (asin-qd arg))
(true (div-qd +qd-pi+ (make-qd-d 6d0))))
- (format t "asin(1/2) = ~/qdi::qd-format/~%" y)
- (format t "pi/6 = ~/qdi::qd-format/~%" true)
+ (format t "asin(1/2) = ~/octi::qd-format/~%" y)
+ (format t "pi/6 = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
(let* ((arg (sqrt-qd (make-qd-d 0.5d0)))
(y (asin-qd arg))
(true (div-qd +qd-pi+ (make-qd-d 4d0))))
- (format t "asin(1/sqrt(2))= ~/qdi::qd-format/~%" y)
- (format t "pi/4 = ~/qdi::qd-format/~%" true)
+ (format t "asin(1/sqrt(2))= ~/octi::qd-format/~%" y)
+ (format t "pi/4 = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
(let* ((arg (div-qd (sqrt-qd (make-qd-d 3d0)) (make-qd-d 2d0)))
(y (asin-qd arg))
(true (div-qd +qd-pi+ (make-qd-d 3d0))))
- (format t "asin(sqrt(3)/2)= ~/qdi::qd-format/~%" y)
- (format t "pi/3 = ~/qdi::qd-format/~%" true)
+ (format t "asin(sqrt(3)/2)= ~/octi::qd-format/~%" y)
+ (format t "pi/3 = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
)
@@ -319,22 +319,22 @@
(let* ((arg (make-qd-d 2d0))
(y (funcall f arg))
(true +qd-log2+))
- (format t "log(2) = ~/qdi::qd-format/~%" y)
- (format t "log(2) = ~/qdi::qd-format/~%" true)
+ (format t "log(2) = ~/octi::qd-format/~%" y)
+ (format t "log(2) = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
(let* ((arg (make-qd-d 10d0))
(y (funcall f arg))
(true +qd-log10+))
- (format t "log(10) = ~/qdi::qd-format/~%" y)
- (format t "log(10) = ~/qdi::qd-format/~%" true)
+ (format t "log(10) = ~/octi::qd-format/~%" y)
+ (format t "log(10) = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
(let* ((arg (add-d-qd 1d0 (scale-float-qd (make-qd-d 1d0) -80)))
(y (funcall f arg))
(true (qd-from-string "8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25")))
- (format t "log(1+2^-80) = ~/qdi::qd-format/~%" y)
- (format t "log(1+2^-80) = ~/qdi::qd-format/~%" true)
+ (format t "log(1+2^-80) = ~/octi::qd-format/~%" y)
+ (format t "log(1+2^-80) = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
)
@@ -344,22 +344,22 @@
(let* ((arg (make-qd-d 1d0))
(y (funcall f arg))
(true +qd-log2+))
- (format t "log1p(1) = ~/qdi::qd-format/~%" y)
- (format t "log(2) = ~/qdi::qd-format/~%" true)
+ (format t "log1p(1) = ~/octi::qd-format/~%" y)
+ (format t "log(2) = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
(let* ((arg (make-qd-d 9d0))
(y (funcall f arg))
(true (qd-from-string "2.3025850929940456840179914546843642076011014886287729760333279009675726096773525q0")))
- (format t "log1p(9) = ~/qdi::qd-format/~%" y)
- (format t "log(10) = ~/qdi::qd-format/~%" true)
+ (format t "log1p(9) = ~/octi::qd-format/~%" y)
+ (format t "log(10) = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
(let* ((arg (scale-float-qd (make-qd-d 1d0) -80))
(y (funcall f arg))
(true (qd-from-string "8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25")))
- (format t "log1p(2^-80) = ~/qdi::qd-format/~%" y)
- (format t "log(1+2^-80) = ~/qdi::qd-format/~%" true)
+ (format t "log1p(2^-80) = ~/octi::qd-format/~%" y)
+ (format t "log(1+2^-80) = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
)
@@ -369,23 +369,23 @@
(let* ((arg +qd-zero+)
(y (funcall f arg))
(true +qd-zero+))
- (format t "exp(0)-1 = ~/qdi::qd-format/~%" y)
- (format t "0 = ~/qdi::qd-format/~%" true)
+ (format t "exp(0)-1 = ~/octi::qd-format/~%" y)
+ (format t "0 = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
(let* ((arg +qd-one+)
(y (funcall f arg))
(true (qd-from-string "1.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952q0")))
- (format t "exp(1)-1 = ~/qdi::qd-format/~%" y)
- (format t "e-1 = ~/qdi::qd-format/~%" true)
+ (format t "exp(1)-1 = ~/octi::qd-format/~%" y)
+ (format t "e-1 = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
(let* ((arg (scale-float-qd +qd-one+ -100))
(y (funcall f arg))
(true (qd-from-string "7.888609052210118054117285652830973804370994921943802079729680186943164342372119432861876389514693341738324702996270767390039172777809233288470357147q-31")))
- (format t "exp(2^-80)-1 = ~/qdi::qd-format/~%" y)
- (format t "exp(2^-80)-1 = ~/qdi::qd-format/~%" true)
+ (format t "exp(2^-80)-1 = ~/octi::qd-format/~%" y)
+ (format t "exp(2^-80)-1 = ~/octi::qd-format/~%" true)
(format t "bits = ~,1f~%"
(bit-accuracy y true)))
--- /project/oct/cvsroot/oct/rt-tests.lisp 2007/10/15 15:45:33 1.4
+++ /project/oct/cvsroot/oct/rt-tests.lisp 2007/10/15 18:21:47 1.5
@@ -50,7 +50,7 @@
;; Pi via Machin's formula
(rt:deftest oct.pi.machin
(let* ((*standard-output* *null*)
- (val (make-instance 'qd-real :value (qdi::test2 nil)))
+ (val (make-instance 'qd-real :value (octi::test2 nil)))
(true oct:+pi+))
(check-accuracy 213 val true))
nil)
@@ -58,7 +58,7 @@
;; Pi via Salamin-Brent algorithm
(rt:deftest oct.pi.salamin-brent
(let* ((*standard-output* *null*)
- (val (make-instance 'qd-real :value (qdi::test3 nil)))
+ (val (make-instance 'qd-real :value (octi::test3 nil)))
(true oct:+pi+))
(check-accuracy 202 val true))
nil)
@@ -66,7 +66,7 @@
;; Pi via Borweign's Quartic formula
(rt:deftest oct.pi.borweign
(let* ((*standard-output* *null*)
- (val (make-instance 'qd-real :value (qdi::test4 nil)))
+ (val (make-instance 'qd-real :value (octi::test4 nil)))
(true oct:+pi+))
(check-accuracy 211 val true))
nil)
@@ -74,16 +74,16 @@
;; e via Taylor series
(rt:deftest oct.e.taylor
(let* ((*standard-output* *null*)
- (val (make-instance 'qd-real :value (qdi::test5 nil)))
- (true (make-instance 'qd-real :value qdi::+qd-e+)))
+ (val (make-instance 'qd-real :value (octi::test5 nil)))
+ (true (make-instance 'qd-real :value octi::+qd-e+)))
(check-accuracy 212 val true))
nil)
;; log(2) via Taylor series
(rt:deftest oct.log2.taylor
(let* ((*standard-output* *null*)
- (val (make-instance 'qd-real :value (qdi::test6 nil)))
- (true (make-instance 'qd-real :value qdi::+qd-log2+)))
+ (val (make-instance 'qd-real :value (octi::test6 nil)))
+ (true (make-instance 'qd-real :value octi::+qd-log2+)))
(check-accuracy 212 val true))
nil)
@@ -126,7 +126,7 @@
(defun atan-qd/duplication (arg)
(make-instance 'qd-real
- :value (qdi::atan-qd/duplication (qd-value arg))))
+ :value (octi::atan-qd/duplication (qd-value arg))))
;;; Tests of atan where we know the analytical result. Same tests,
;;; but using the atan duplication formula.
@@ -139,15 +139,15 @@
(rt:deftest oct.atan/dup.2
(let* ((arg (sqrt #q3))
- (y (/ (atan-qd/duplication arg) +pi+))
- (true (/ #q3)))
+ (y (/ (atan-qd/duplication arg) +pi+))
+ (true (/ #q3)))
(check-accuracy 212 y true))
nil)
(rt:deftest oct.atan/dup.3
(let* ((arg #q1)
- (y (/ (atan-qd/duplication arg) +pi+))
- (true (/ #q4)))
+ (y (/ (atan-qd/duplication arg) +pi+))
+ (true (/ #q4)))
(check-accuracy 212 y true))
nil)
@@ -160,8 +160,8 @@
(rt:deftest oct.atan/dup.5
(let* ((arg #q-1q100)
- (y (/ (atan-qd/duplication arg) +pi+))
- (true #q-.5))
+ (y (/ (atan-qd/duplication arg) +pi+))
+ (true #q-.5))
(check-accuracy 212 y true))
nil)
@@ -169,7 +169,7 @@
;;; but using a CORDIC implementation.
(defun atan-qd/cordic (arg)
(make-instance 'qd-real
- :value (qdi::atan-qd/cordic (qd-value arg))))
+ :value (octi::atan-qd/cordic (qd-value arg))))
(rt:deftest oct.atan/cordic.1
(let* ((arg (/ (sqrt #q3)))
@@ -180,15 +180,15 @@
(rt:deftest oct.atan/cordic.2
(let* ((arg (sqrt #q3))
- (y (/ (atan-qd/cordic arg) +pi+))
- (true (/ #q3)))
+ (y (/ (atan-qd/cordic arg) +pi+))
+ (true (/ #q3)))
(check-accuracy 212 y true))
nil)
(rt:deftest oct.atan/cordic.3
(let* ((arg #q1)
- (y (/ (atan-qd/cordic arg) +pi+))
- (true (/ #q4)))
+ (y (/ (atan-qd/cordic arg) +pi+))
+ (true (/ #q4)))
(check-accuracy 212 y true))
nil)
@@ -201,8 +201,8 @@
(rt:deftest oct.atan/cordic.5
(let* ((arg #q-1q100)
- (y (/ (atan-qd/cordic arg) +pi+))
- (true #q-.5))
+ (y (/ (atan-qd/cordic arg) +pi+))
+ (true #q-.5))
(check-accuracy 212 y true))
nil)
@@ -210,25 +210,39 @@
;;; Tests of sin where we know the analytical result.
(rt:deftest oct.sin.1
(let* ((arg (/ +pi+ 6))
- (y (sin arg))
- (true #q.5))
+ (y (sin arg))
+ (true #q.5))
(check-accuracy 212 y true))
nil)
(rt:deftest oct.sin.2
(let* ((arg (/ +pi+ 4))
- (y (sin arg))
- (true (sqrt #q.5)))
+ (y (sin arg))
+ (true (sqrt #q.5)))
(check-accuracy 212 y true))
nil)
(rt:deftest oct.sin.3
(let* ((arg (/ +pi+ 3))
- (y (sin arg))
- (true (/ (sqrt #q3) 2)))
+ (y (sin arg))
+ (true (/ (sqrt #q3) 2)))
(check-accuracy 212 y true))
nil)
+(rt:deftest oct.big-sin.1
+ (let* ((arg (oct:make-qd (ash 1 120)))
+ (y (sin arg))
+ (true #q3.778201093607520226555484700569229919605866976512306642257987199414885q-1))
+ (check-accuracy 205 y true))
+ nil)
+
+(rt:deftest oct.big-sin.2
+ (let* ((arg (oct:make-qd (ash 1 1023)))
+ (y (sin arg))
+ (true #q5.631277798508840134529434079444683477103854907361251399182750155357133q-1))
+ (check-accuracy 205 y true))
+ nil)
+
;;; Tests of tan where we know the analytical result.
(rt:deftest oct.tan.1
(let* ((arg (/ +pi+ 6))
@@ -256,7 +270,7 @@
(defun tan/cordic (arg)
(make-instance 'qd-real
- :value (qdi::tan-qd/cordic (qd-value arg))))
+ :value (octi::tan-qd/cordic (qd-value arg))))
(rt:deftest oct.tan/cordic.1
(let* ((arg (/ +pi+ 6))
@@ -307,14 +321,14 @@
(rt:deftest oct.log.1
(let* ((arg #q2)
(y (log arg))
- (true (make-instance 'qd-real :value qdi::+qd-log2+)))
+ (true (make-instance 'qd-real :value octi::+qd-log2+)))
(check-accuracy 212 y true))
nil)
(rt:deftest oct.log.2
(let* ((arg #q10)
(y (log arg))
- (true (make-instance 'qd-real :value qdi::+qd-log10+)))
+ (true (make-instance 'qd-real :value octi::+qd-log10+)))
(check-accuracy 207 y true))
nil)
@@ -329,19 +343,19 @@
(defun log/newton (arg)
(make-instance 'qd-real
- :value (qdi::log-qd/newton (qd-value arg))))
+ :value (octi::log-qd/newton (qd-value arg))))
(rt:deftest oct.log/newton.1
(let* ((arg #q2)
(y (log/newton arg))
- (true (make-instance 'qd-real :value qdi::+qd-log2+)))
+ (true (make-instance 'qd-real :value octi::+qd-log2+)))
(check-accuracy 212 y true))
nil)
(rt:deftest oct.log/newton.2
(let* ((arg #q10)
(y (log/newton arg))
- (true (make-instance 'qd-real :value qdi::+qd-log10+)))
+ (true (make-instance 'qd-real :value octi::+qd-log10+)))
(check-accuracy 207 y true))
nil)
@@ -356,19 +370,19 @@
(defun log/agm (arg)
(make-instance 'qd-real
- :value (qdi::log-qd/agm (qd-value arg))))
+ :value (octi::log-qd/agm (qd-value arg))))
(rt:deftest oct.log/agm.1
(let* ((arg #q2)
(y (log/agm arg))
- (true (make-instance 'qd-real :value qdi::+qd-log2+)))
+ (true (make-instance 'qd-real :value octi::+qd-log2+)))
(check-accuracy 203 y true))
nil)
(rt:deftest oct.log/agm.2
(let* ((arg #q10)
(y (log/agm arg))
- (true (make-instance 'qd-real :value qdi::+qd-log10+)))
+ (true (make-instance 'qd-real :value octi::+qd-log10+)))
(check-accuracy 205 y true))
nil)
@@ -383,19 +397,19 @@
(defun log/agm2 (arg)
(make-instance 'qd-real
- :value (qdi::log-qd/agm2 (qd-value arg))))
+ :value (octi::log-qd/agm2 (qd-value arg))))
(rt:deftest oct.log/agm2.1
(let* ((arg #q2)
(y (log/agm2 arg))
- (true (make-instance 'qd-real :value qdi::+qd-log2+)))
+ (true (make-instance 'qd-real :value octi::+qd-log2+)))
(check-accuracy 203 y true))
nil)
(rt:deftest oct.log/agm2.2
(let* ((arg #q10)
(y (log/agm2 arg))
- (true (make-instance 'qd-real :value qdi::+qd-log10+)))
+ (true (make-instance 'qd-real :value octi::+qd-log10+)))
(check-accuracy 205 y true))
nil)
@@ -409,19 +423,19 @@
;;; Tests of log using AGM3, a faster variation of AGM2.
(defun log/agm3 (arg)
(make-instance 'qd-real
- :value (qdi::log-qd/agm3 (qd-value arg))))
+ :value (octi::log-qd/agm3 (qd-value arg))))
(rt:deftest oct.log/agm3.1
(let* ((arg #q2)
(y (log/agm3 arg))
- (true (make-instance 'qd-real :value qdi::+qd-log2+)))
+ (true (make-instance 'qd-real :value octi::+qd-log2+)))
(check-accuracy 203 y true))
nil)
(rt:deftest oct.log/agm3.2
(let* ((arg #q10)
(y (log/agm3 arg))
- (true (make-instance 'qd-real :value qdi::+qd-log10+)))
+ (true (make-instance 'qd-real :value octi::+qd-log10+)))
(check-accuracy 205 y true))
nil)
@@ -474,7 +488,7 @@
(defun log1p/dup (arg)
(make-instance 'qd-real
- :value (qdi::log1p-qd/duplication (qd-value arg))))
+ :value (octi::log1p-qd/duplication (qd-value arg))))
(rt:deftest oct.log1p.1
(let* ((arg #q9)
@@ -495,7 +509,7 @@
(defun expm1/series (arg)
(make-instance 'qd-real
- :value (qdi::expm1-qd/series (qd-value arg))))
+ :value (octi::expm1-qd/series (qd-value arg))))
(rt:deftest oct.expm1/series.1
(let* ((arg #q0)
@@ -522,7 +536,7 @@
(defun expm1/dup (arg)
(make-instance 'qd-real
- :value (qdi::expm1-qd/duplication (qd-value arg))))
+ :value (octi::expm1-qd/duplication (qd-value arg))))
(rt:deftest oct.expm1/dup.1
@@ -551,7 +565,7 @@
;; thing.
(rt:deftest oct.integer-decode.1
(multiple-value-bind (frac exp s)
- (qdi:integer-decode-qd (qdi::%make-qd-d -0.03980126756814893d0
+ (octi:integer-decode-qd (octi::%make-qd-d -0.03980126756814893d0
-2.7419792323327893d-18
0d0 0d0))
(unless (and (eql frac 103329998279901916046530991816704)
--- /project/oct/cvsroot/oct/tests.lisp 2007/10/13 02:14:43 1.2
+++ /project/oct/cvsroot/oct/tests.lisp 2007/10/15 18:21:47 1.3
@@ -43,10 +43,10 @@
(format t "bits: ~,1f~%" (bit-accuracy est true)))
(defconstant +e+
- (make-instance 'qd-real :value qdi::+qd-e+))
+ (make-instance 'qd-real :value octi::+qd-e+))
(defconstant +log2+
- (make-instance 'qd-real :value qdi::+qd-log2+))
+ (make-instance 'qd-real :value octi::+qd-log2+))
(defun test2 ()
;; pi/4 = 4 * arctan(1/5) - arctan(1/239)
@@ -268,7 +268,7 @@
(print-result y true))
(let* ((arg #q10)
(y (log arg))
- (true (make-instance 'qd-real :value qdi::+qd-log10+)))
+ (true (make-instance 'qd-real :value octi::+qd-log10+)))
(format t "log(10)~%")
(print-result y true))
(let* ((arg (+ 1 (scale-float #q1 -80)))
@@ -287,7 +287,7 @@
(destructuring-bind (arg true)
f
(let ((y (sqrt arg)))
- (format t "sqrt(~/qdi::qd-format/)~%" (qd-value arg))
+ (format t "sqrt(~/octi::qd-format/)~%" (qd-value arg))
(print-result y true)))))
(defun all-tests ()
--- /project/oct/cvsroot/oct/timing.lisp 2007/08/27 17:49:19 1.2
+++ /project/oct/cvsroot/oct/timing.lisp 2007/10/15 18:21:47 1.3
@@ -38,9 +38,9 @@
(setf sum (cl:+ sum 1d0)))
sum))
(sum-%qd ()
- (let ((sum (qdi::make-qd-d 0d0))
- (one (qdi::make-qd-d 1d0)))
- (declare (type qdi::%quad-double sum)
+ (let ((sum (octi::make-qd-d 0d0))
+ (one (octi::make-qd-d 1d0)))
+ (declare (type octi::%quad-double sum)
(optimize (speed 3)))
(dotimes (k n)
(declare (fixnum k))
@@ -76,9 +76,9 @@
(setf sum (cl:* sum 1d0)))
sum))
(mul-%qd ()
- (let ((sum (qdi::make-qd-d 0d0))
- (one (qdi::make-qd-d 1d0)))
- (declare (type qdi::%quad-double sum)
+ (let ((sum (octi::make-qd-d 0d0))
+ (one (octi::make-qd-d 1d0)))
+ (declare (type octi::%quad-double sum)
(optimize (speed 3)))
(dotimes (k n)
(declare (fixnum k))
@@ -113,9 +113,9 @@
(setf sum (cl:/ sum 1d0)))
sum))
(div-%qd ()
- (let ((sum (qdi::make-qd-d 7d0))
- (one (qdi::make-qd-d 1d0)))
- (declare (type qdi::%quad-double sum)
+ (let ((sum (octi::make-qd-d 7d0))
+ (one (octi::make-qd-d 1d0)))
+ (declare (type octi::%quad-double sum)
(optimize (speed 3)))
(dotimes (k n)
(declare (fixnum k))
@@ -150,8 +150,8 @@
(setf sum (cl:sqrt sum)))
sum))
(sqrt-%qd ()
- (let ((sum (qdi::make-qd-d 7d0)))
- (declare (type qdi::%quad-double sum)
+ (let ((sum (octi::make-qd-d 7d0)))
+ (declare (type octi::%quad-double sum)
(optimize (speed 3)))
(dotimes (k n)
(declare (fixnum k))
1
0
Update of /project/oct/cvsroot/oct
In directory clnet:/tmp/cvs-serv28447
Modified Files:
qd.lisp rt-tests.lisp
Log Message:
qd.lisp:
o Oops. In INTEGER-DECODE-QD, the signs of the parts were not
computed correctly when combining them into the final integer
result.
rt-tests.lisp:
o Add a test for INTEGER-DECODE-QD.
o Use OCT as the package, not QD.
--- /project/oct/cvsroot/oct/qd.lisp 2007/10/13 15:34:51 1.55
+++ /project/oct/cvsroot/oct/qd.lisp 2007/10/15 15:45:33 1.56
@@ -1135,19 +1135,19 @@
(multiple-value-bind (q1-int q1-exp q1-sign)
(integer-decode-float (qd-1 q))
(setf q0-int (+ (ash q0-int (- q0-exp q1-exp))
- (* q1-sign q1-int)))
+ (* q0-sign q1-sign q1-int)))
(when (zerop (qd-2 q))
(return-from integer-decode-qd (values q0-int q1-exp q0-sign)))
(multiple-value-bind (q2-int q2-exp q2-sign)
(integer-decode-float (qd-2 q))
(setf q0-int (+ (ash q0-int (- q1-exp q2-exp))
- (* q2-sign q2-int)))
+ (* q0-sign q2-sign q2-int)))
(when (zerop (qd-3 q))
(return-from integer-decode-qd (values q0-int q2-exp q0-sign)))
(multiple-value-bind (q3-int q3-exp q3-sign)
(integer-decode-float (qd-3 q))
(values (+ (ash q0-int (- q2-exp q3-exp))
- (* q3-sign q3-int))
+ (* q0-sign q3-sign q3-int))
q3-exp
q0-sign))))))
--- /project/oct/cvsroot/oct/rt-tests.lisp 2007/10/13 02:14:43 1.3
+++ /project/oct/cvsroot/oct/rt-tests.lisp 2007/10/15 15:45:33 1.4
@@ -51,7 +51,7 @@
(rt:deftest oct.pi.machin
(let* ((*standard-output* *null*)
(val (make-instance 'qd-real :value (qdi::test2 nil)))
- (true qd:+pi+))
+ (true oct:+pi+))
(check-accuracy 213 val true))
nil)
@@ -59,7 +59,7 @@
(rt:deftest oct.pi.salamin-brent
(let* ((*standard-output* *null*)
(val (make-instance 'qd-real :value (qdi::test3 nil)))
- (true qd:+pi+))
+ (true oct:+pi+))
(check-accuracy 202 val true))
nil)
@@ -67,7 +67,7 @@
(rt:deftest oct.pi.borweign
(let* ((*standard-output* *null*)
(val (make-instance 'qd-real :value (qdi::test4 nil)))
- (true qd:+pi+))
+ (true oct:+pi+))
(check-accuracy 211 val true))
nil)
@@ -545,3 +545,18 @@
(true #q7.888609052210118054117285652830973804370994921943802079729680186943164342372119432861876389514693341738324702996270767390039172777809233288470357147q-31))
(check-accuracy 211 y true))
nil)
+
+;; If we screw up integer-decode-qd, printing is wrong. Here is one
+;; case where integer-decode-qd was screwed up and printing the wrong
+;; thing.
+(rt:deftest oct.integer-decode.1
+ (multiple-value-bind (frac exp s)
+ (qdi:integer-decode-qd (qdi::%make-qd-d -0.03980126756814893d0
+ -2.7419792323327893d-18
+ 0d0 0d0))
+ (unless (and (eql frac 103329998279901916046530991816704)
+ (eql exp -111)
+ (eql s -1))
+ (list frac exp s)))
+ nil)
+
1
0
Update of /project/oct/cvsroot/oct
In directory clnet:/tmp/cvs-serv29190
Modified Files:
qd-fun.lisp
Log Message:
o Fix some typos accurate-sin-qd and accurate-cos-qd.
o Adjust code in accurate-sin-qd and accurate-cos-qd to handle values
of 0 <= j <= 3, instead of -1 <= j <= 2.
o Add accurate-sincos-qd.
--- /project/oct/cvsroot/oct/qd-fun.lisp 2007/10/13 20:14:45 1.82
+++ /project/oct/cvsroot/oct/qd-fun.lisp 2007/10/15 03:26:38 1.83
@@ -666,66 +666,62 @@
;; cos(s+k*pi/1024) = cos(s)*cos(k*pi/1024)
;; - sin(s)*sin(k*pi/1024)
(when (zerop-qd a)
- (return-from sin-qd a))
+ (return-from accurate-sin-qd a))
(multiple-value-bind (j r)
(rem-pi/2-int a)
- (let* ((j (if (= j 3)
- -1
- j))
- (abs-j (abs j)))
- (multiple-value-bind (k tmp)
- (divrem-qd r +qd-pi/1024+)
- (let* ((k (truncate (qd-0 k)))
- (abs-k (abs k)))
- (assert (<= abs-j 2))
- (assert (<= abs-k 256))
- ;; Compute sin(s) and cos(s)
- (multiple-value-bind (sin-t cos-t)
- (sincos-taylor tmp)
- (multiple-value-bind (s c)
- (cond ((zerop abs-k)
- (values sin-t cos-t))
- (t
- ;; Compute sin(s+k*pi/1024), cos(s+k*pi/1024)
- (let ((u (aref +qd-cos-table+ (cl:1- abs-k)))
- (v (aref +qd-sin-table+ (cl:1- abs-k))))
- (cond ((plusp k)
- ;; sin(s) * cos(k*pi/1024)
- ;; + cos(s) * sin(k*pi/1024)
- ;;
- ;; cos(s) * cos(k*pi/1024)
- ;; - sin(s) * sin(k*pi/1024)
- (values (add-qd (mul-qd u sin-t)
- (mul-qd v cos-t))
- (sub-qd (mul-qd u cos-t)
- (mul-qd v sin-t))))
- (t
- ;; sin(s) * cos(k*pi/1024)
- ;; - cos(s) * sin(|k|*pi/1024)
- ;;
- ;; cos(s) * cos(k*pi/1024)
- ;; + sin(s) * sin(|k|*pi/1024)
- (values (sub-qd (mul-qd u sin-t)
- (mul-qd v cos-t))
- (add-qd (mul-qd u cos-t)
- (mul-qd v sin-t))))))))
- ;;(format t "s = ~/qd::qd-format/~%" s)
- ;;(format t "c = ~/qd::qd-format/~%" c)
- ;; sin(x) = sin(s+k*pi/1024) * cos(j*pi/2)
- ;; + cos(s+k*pi/1024) * sin(j*pi/2)
- (cond ((zerop abs-j)
- ;; cos(j*pi/2) = 1, sin(j*pi/2) = 0
- s)
- ((= j 1)
- ;; cos(j*pi/2) = 0, sin(j*pi/2) = 1
- c)
- ((= j -1)
- ;; cos(j*pi/2) = 0, sin(j*pi/2) = -1
- (neg-qd c))
+ (multiple-value-bind (k tmp)
+ (divrem-qd r +qd-pi/1024+)
+ (let* ((k (truncate (qd-0 k)))
+ (abs-k (abs k)))
+ (assert (<= 0 j 3))
+ (assert (<= abs-k 256))
+ ;; Compute sin(s) and cos(s)
+ (multiple-value-bind (sin-t cos-t)
+ (sincos-taylor tmp)
+ (multiple-value-bind (s c)
+ (cond ((zerop abs-k)
+ (values sin-t cos-t))
(t
- ;; cos(j*pi/2) = -1, sin(j*pi/2) = 0
- (neg-qd s))))))))))
+ ;; Compute sin(s+k*pi/1024), cos(s+k*pi/1024)
+ (let ((u (aref +qd-cos-table+ (cl:1- abs-k)))
+ (v (aref +qd-sin-table+ (cl:1- abs-k))))
+ (cond ((plusp k)
+ ;; sin(s) * cos(k*pi/1024)
+ ;; + cos(s) * sin(k*pi/1024)
+ ;;
+ ;; cos(s) * cos(k*pi/1024)
+ ;; - sin(s) * sin(k*pi/1024)
+ (values (add-qd (mul-qd u sin-t)
+ (mul-qd v cos-t))
+ (sub-qd (mul-qd u cos-t)
+ (mul-qd v sin-t))))
+ (t
+ ;; sin(s) * cos(k*pi/1024)
+ ;; - cos(s) * sin(|k|*pi/1024)
+ ;;
+ ;; cos(s) * cos(k*pi/1024)
+ ;; + sin(s) * sin(|k|*pi/1024)
+ (values (sub-qd (mul-qd u sin-t)
+ (mul-qd v cos-t))
+ (add-qd (mul-qd u cos-t)
+ (mul-qd v sin-t))))))))
+ ;;(format t "s = ~/qd::qd-format/~%" s)
+ ;;(format t "c = ~/qd::qd-format/~%" c)
+ ;; sin(x) = sin(s+k*pi/1024) * cos(j*pi/2)
+ ;; + cos(s+k*pi/1024) * sin(j*pi/2)
+ (cond ((zerop j)
+ ;; cos(j*pi/2) = 1, sin(j*pi/2) = 0
+ s)
+ ((= j 1)
+ ;; cos(j*pi/2) = 0, sin(j*pi/2) = 1
+ c)
+ ((= j 2)
+ ;; cos(j*pi/2) = -1, sin(j*pi/2) = 0
+ (neg-qd s))
+ ((= j 3)
+ ;; cos(j*pi/2) = 0, sin(j*pi/2) = -1
+ (neg-qd c)))))))))
(defun accurate-cos-qd (a)
"Cos(a)"
@@ -748,68 +744,130 @@
;; cos(s+k*pi/1024) = cos(s)*cos(k*pi/1024)
;; - sin(s)*sin(k*pi/1024)
(when (zerop-qd a)
- (return-from cos-qd +qd-one+))
+ (return-from accurate-cos-qd +qd-one+))
(multiple-value-bind (j r)
(rem-pi/2-int a)
- (let* ((j (if (= j 3)
- -1
- j))
- (abs-j (abs j)))
- (multiple-value-bind (k tmp)
- (divrem-qd r +qd-pi/1024+)
- (let* ((k (truncate (qd-0 k)))
- (abs-k (abs k)))
- (assert (<= abs-j 2))
- (assert (<= abs-k 256))
- ;; Compute sin(s) and cos(s)
- (multiple-value-bind (sin-t cos-t)
- (sincos-taylor tmp)
- (multiple-value-bind (s c)
- (cond ((zerop abs-k)
- (values sin-t cos-t))
- (t
- ;; Compute sin(s+k*pi/1024), cos(s+k*pi/1024)
- (let ((u (aref +qd-cos-table+ (cl:1- abs-k)))
- (v (aref +qd-sin-table+ (cl:1- abs-k))))
- (cond ((plusp k)
- ;; sin(s) * cos(k*pi/1024)
- ;; + cos(s) * sin(k*pi/1024)
- ;;
- ;; cos(s) * cos(k*pi/1024)
- ;; - sin(s) * sin(k*pi/1024)
- (values (add-qd (mul-qd u sin-t)
- (mul-qd v cos-t))
- (sub-qd (mul-qd u cos-t)
- (mul-qd v sin-t))))
- (t
- ;; sin(s) * cos(k*pi/1024)
- ;; - cos(s) * sin(|k|*pi/1024)
- ;;
- ;; cos(s) * cos(k*pi/1024)
- ;; + sin(s) * sin(|k|*pi/1024)
- (values (sub-qd (mul-qd u sin-t)
- (mul-qd v cos-t))
- (add-qd (mul-qd u cos-t)
- (mul-qd v sin-t))))))))
- #+nil
- (progn
- (format t "s = ~/qd::qd-format/~%" s)
- (format t "c = ~/qd::qd-format/~%" c))
- ;; sin(x) = sin(s+k*pi/1024) * cos(j*pi/2)
- ;; + cos(s+k*pi/1024) * sin(j*pi/2)
- (cond ((zerop abs-j)
- ;; cos(j*pi/2) = 1, sin(j*pi/2) = 0
- c)
- ((= j 1)
- ;; cos(j*pi/2) = 0, sin(j*pi/2) = 1
- (neg-qd s))
- ((= j -1)
- ;; cos(j*pi/2) = 0, sin(j*pi/2) = -1
- s)
+ (multiple-value-bind (k tmp)
+ (divrem-qd r +qd-pi/1024+)
+ (let* ((k (truncate (qd-0 k)))
+ (abs-k (abs k)))
+ (assert (<= 0 j 3))
+ (assert (<= abs-k 256))
+ ;; Compute sin(s) and cos(s)
+ (multiple-value-bind (sin-t cos-t)
+ (sincos-taylor tmp)
+ (multiple-value-bind (s c)
+ (cond ((zerop abs-k)
+ (values sin-t cos-t))
(t
- ;; cos(j*pi/2) = -1, sin(j*pi/2) = 0
- (neg-qd c))))))))))
+ ;; Compute sin(s+k*pi/1024), cos(s+k*pi/1024)
+ (let ((u (aref +qd-cos-table+ (cl:1- abs-k)))
+ (v (aref +qd-sin-table+ (cl:1- abs-k))))
+ (cond ((plusp k)
+ ;; sin(s) * cos(k*pi/1024)
+ ;; + cos(s) * sin(k*pi/1024)
+ ;;
+ ;; cos(s) * cos(k*pi/1024)
+ ;; - sin(s) * sin(k*pi/1024)
+ (values (add-qd (mul-qd u sin-t)
+ (mul-qd v cos-t))
+ (sub-qd (mul-qd u cos-t)
+ (mul-qd v sin-t))))
+ (t
+ ;; sin(s) * cos(k*pi/1024)
+ ;; - cos(s) * sin(|k|*pi/1024)
+ ;;
+ ;; cos(s) * cos(k*pi/1024)
+ ;; + sin(s) * sin(|k|*pi/1024)
+ (values (sub-qd (mul-qd u sin-t)
+ (mul-qd v cos-t))
+ (add-qd (mul-qd u cos-t)
+ (mul-qd v sin-t))))))))
+ #+nil
+ (progn
+ (format t "s = ~/qd::qd-format/~%" s)
+ (format t "c = ~/qd::qd-format/~%" c))
+ ;; sin(x) = sin(s+k*pi/1024) * cos(j*pi/2)
+ ;; + cos(s+k*pi/1024) * sin(j*pi/2)
+ (cond ((zerop j)
+ ;; cos(j*pi/2) = 1, sin(j*pi/2) = 0
+ c)
+ ((= j 1)
+ ;; cos(j*pi/2) = 0, sin(j*pi/2) = 1
+ (neg-qd s))
+ ((= j 2)
+ ;; cos(j*pi/2) = -1, sin(j*pi/2) = 0
+ (neg-qd c))
+ ((= j 3)
+ ;; cos(j*pi/2) = 0, sin(j*pi/2) = -1
+ s))))))))
+
+(defun accurate-sincos-qd (a)
+ (declare (type %quad-double a))
+ (when (zerop-qd a)
+ (return-from accurate-sincos-qd
+ (values +qd-zero+
+ +qd-one+)))
+
+ (multiple-value-bind (j r)
+ (rem-pi/2-int a)
+ (multiple-value-bind (k tmp)
+ (divrem-qd tmp +qd-pi/1024+)
+ (let* ((k (truncate (qd-0 k)))
+ (abs-j (abs j))
+ (abs-k (abs k)))
+ (assert (<= abs-j 2))
+ (assert (<= abs-k 256))
+ ;; Compute sin(s) and cos(s)
+ (multiple-value-bind (sin-t cos-t)
+ (sincos-taylor tmp)
+ (multiple-value-bind (s c)
+ (cond ((zerop abs-k)
+ (values sin-t cos-t))
+ (t
+ ;; Compute sin(s+k*pi/1024), cos(s+k*pi/1024)
+ (let ((u (aref +qd-cos-table+ (cl:1- abs-k)))
+ (v (aref +qd-sin-table+ (cl:1- abs-k))))
+ (cond ((plusp k)
+ ;; sin(s) * cos(k*pi/1024)
+ ;; + cos(s) * sin(k*pi/1024)
+ ;;
+ ;; cos(s) * cos(k*pi/1024)
+ ;; - sin(s) * sin(k*pi/1024)
+ (values (add-qd (mul-qd u sin-t)
+ (mul-qd v cos-t))
+ (sub-qd (mul-qd u cos-t)
+ (mul-qd v sin-t))))
+ (t
+ ;; sin(s) * cos(k*pi/1024)
+ ;; - cos(s) * sin(|k|*pi/1024)
+ ;;
+ ;; cos(s) * cos(k*pi/1024)
+ ;; + sin(s) * sin(|k|*pi/1024)
+ (values (sub-qd (mul-qd u sin-t)
+ (mul-qd v cos-t))
+ (add-qd (mul-qd u cos-t)
+ (mul-qd v sin-t))))))))
+ #+nil
+ (progn
+ (format t "s = ~/qd::qd-format/~%" s)
+ (format t "c = ~/qd::qd-format/~%" c))
+ ;; sin(x) = sin(s+k*pi/1024) * cos(j*pi/2)
+ ;; + cos(s+k*pi/1024) * sin(j*pi/2)
+ (cond ((zerop abs-j)
+ ;; cos(j*pi/2) = 1, sin(j*pi/2) = 0
+ (values s c))
+ ((= j 1)
+ ;; cos(j*pi/2) = 0, sin(j*pi/2) = 1
+ (values c (neg-qd s)))
+ ((= j 2)
+ ;; cos(j*pi/2) = -1, sin(j*pi/2) = 0
+ (values (neg-qd s)
+ (neg-qd c)))
+ ((= j 3)
+ ;; cos(j*pi/2) = 0, sin(j*pi/2) = -1
+ (values (neg-qd c) s))))))))))
(defun atan2-qd/newton (y x)
(declare (type %quad-double y x)
1
0
Update of /project/oct/cvsroot/oct
In directory clnet:/tmp/cvs-serv15775
Modified Files:
qd-const.lisp
Log Message:
Add +2/pi-bits+, 1584 bits of 2/pi.
--- /project/oct/cvsroot/oct/qd-const.lisp 2007/10/11 17:47:08 1.17
+++ /project/oct/cvsroot/oct/qd-const.lisp 2007/10/14 18:38:14 1.18
@@ -31,6 +31,12 @@
(defconstant +qd-one+
(make-qd-d 1d0))
+;; The bits of 2/pi. Scale these bits by 2^(-1584) and you'll get
+;; 2/pi. These are used for accurate argument reduction for the trig
+;; functions.
+(defconstant +2/pi-bits+
+ #xA2F9836E4E441529FC2757D1F534DDC0DB6295993C439041FE5163ABDEBBC561B7246E3A424DD2E006492EEA09D1921CFE1DEB1CB129A73EE88235F52EBB4484E99C7026B45F7E413991D639835339F49C845F8BBDF9283B1FF897FFDE05980FEF2F118B5A0A6D1F6D367ECF27CB09B74F463F669E5FEA2D7527BAC7EBE5F17B3D0739F78A5292EA6BFB5FB11F8D5D0856033046FC7B6BABF0CFBC209AF4361DA9E391615EE61B086599855F14A068408DFFD8804D73273106061556CA73A8C960E27BC08C6B)
+
;; 3.1415926535897932384626433832795028841971693993751058209749445923078L0
;; #q3.1415926535897932384626433832795028841971693993751058209749445923q0
(defconstant +qd-pi+
1
0
Update of /project/oct/cvsroot/oct
In directory clnet:/tmp/cvs-serv3017
Modified Files:
qd-fun.lisp
Log Message:
o New routine to compute x mod pi/2 accurately, using many bits of
2/pi.
o Implement accurate sin and cos routines to use this new routine.
(Not used yet.)
--- /project/oct/cvsroot/oct/qd-fun.lisp 2007/10/10 15:21:47 1.81
+++ /project/oct/cvsroot/oct/qd-fun.lisp 2007/10/13 20:14:45 1.82
@@ -612,7 +612,205 @@
(values (neg-qd s)
(neg-qd c))))))))))))
-
+;; A more accurate implementation of sin and cos. We use a more
+;; accurate argument reduction using 1584 bits of 2/pi. This makes
+;; sin and cos more expensive, of course.
+
+(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)
+ (integer-decode-qd qd)
+ (let* ((exp (- e 1584))
+ (prod (* (* s i) +2/pi-bits+))
+ (int (ash prod exp))
+ (mod (logand int 3))
+ (frac (ldb (byte 212 (- (- exp) 212)) prod))
+ (f (mul-qd (scale-float-qd (rational-to-qd frac) -212)
+ +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 (<= (abs (qd-0 f)) (/ pi 4))
+ (values mod f)
+ (values (mod (1+ mod) 4) (sub-qd f +qd-pi/2+))))))
+
+(defun accurate-sin-qd (a)
+ "Sin(a)"
+ (declare (type %quad-double a))
+ ;; To compute sin(x), choose integers a, b so that
+ ;;
+ ;; x = s + a * (pi/2) + b*(pi/1024)
+ ;;
+ ;; with |x| <= pi/2048. Using a precomputed table of sin(k*pi/1024)
+ ;; and cos(k*pi/1024), we can compute sin(x) from sin(s) and cos(s).
+ ;;
+ ;; sin(x) = sin(s+k*(pi/1024) + j*pi/2)
+ ;; = sin(s+k*(pi/1024))*cos(j*pi/2)
+ ;; + cos(s+k*(pi/1024))*sin(j*pi/2)
+ ;;
+ ;; sin(s+k*pi/1024) = sin(s)*cos(k*pi/1024)
+ ;; + cos(s)*sin(k*pi/1024)
+ ;;
+ ;; cos(s+k*pi/1024) = cos(s)*cos(k*pi/1024)
+ ;; - sin(s)*sin(k*pi/1024)
+ (when (zerop-qd a)
+ (return-from sin-qd a))
+
+ (multiple-value-bind (j r)
+ (rem-pi/2-int a)
+ (let* ((j (if (= j 3)
+ -1
+ j))
+ (abs-j (abs j)))
+ (multiple-value-bind (k tmp)
+ (divrem-qd r +qd-pi/1024+)
+ (let* ((k (truncate (qd-0 k)))
+ (abs-k (abs k)))
+ (assert (<= abs-j 2))
+ (assert (<= abs-k 256))
+ ;; Compute sin(s) and cos(s)
+ (multiple-value-bind (sin-t cos-t)
+ (sincos-taylor tmp)
+ (multiple-value-bind (s c)
+ (cond ((zerop abs-k)
+ (values sin-t cos-t))
+ (t
+ ;; Compute sin(s+k*pi/1024), cos(s+k*pi/1024)
+ (let ((u (aref +qd-cos-table+ (cl:1- abs-k)))
+ (v (aref +qd-sin-table+ (cl:1- abs-k))))
+ (cond ((plusp k)
+ ;; sin(s) * cos(k*pi/1024)
+ ;; + cos(s) * sin(k*pi/1024)
+ ;;
+ ;; cos(s) * cos(k*pi/1024)
+ ;; - sin(s) * sin(k*pi/1024)
+ (values (add-qd (mul-qd u sin-t)
+ (mul-qd v cos-t))
+ (sub-qd (mul-qd u cos-t)
+ (mul-qd v sin-t))))
+ (t
+ ;; sin(s) * cos(k*pi/1024)
+ ;; - cos(s) * sin(|k|*pi/1024)
+ ;;
+ ;; cos(s) * cos(k*pi/1024)
+ ;; + sin(s) * sin(|k|*pi/1024)
+ (values (sub-qd (mul-qd u sin-t)
+ (mul-qd v cos-t))
+ (add-qd (mul-qd u cos-t)
+ (mul-qd v sin-t))))))))
+ ;;(format t "s = ~/qd::qd-format/~%" s)
+ ;;(format t "c = ~/qd::qd-format/~%" c)
+ ;; sin(x) = sin(s+k*pi/1024) * cos(j*pi/2)
+ ;; + cos(s+k*pi/1024) * sin(j*pi/2)
+ (cond ((zerop abs-j)
+ ;; cos(j*pi/2) = 1, sin(j*pi/2) = 0
+ s)
+ ((= j 1)
+ ;; cos(j*pi/2) = 0, sin(j*pi/2) = 1
+ c)
+ ((= j -1)
+ ;; cos(j*pi/2) = 0, sin(j*pi/2) = -1
+ (neg-qd c))
+ (t
+ ;; cos(j*pi/2) = -1, sin(j*pi/2) = 0
+ (neg-qd s))))))))))
+
+(defun accurate-cos-qd (a)
+ "Cos(a)"
+ ;; Just like sin-qd, but for cos.
+ (declare (type %quad-double a))
+ ;; To compute sin(x), choose integers a, b so that
+ ;;
+ ;; x = s + a * (pi/2) + b*(pi/1024)
+ ;;
+ ;; with |x| <= pi/2048. Using a precomputed table of sin(k*pi/1024)
+ ;; and cos(k*pi/1024), we can compute sin(x) from sin(s) and cos(s).
+ ;;
+ ;; sin(x) = sin(s+k*(pi/1024) + j*pi/2)
+ ;; = sin(s+k*(pi/1024))*cos(j*pi/2)
+ ;; + cos(s+k*(pi/1024))*sin(j*pi/2)
+ ;;
+ ;; sin(s+k*pi/1024) = sin(s)*cos(k*pi/1024)
+ ;; + cos(s)*sin(k*pi/1024)
+ ;;
+ ;; cos(s+k*pi/1024) = cos(s)*cos(k*pi/1024)
+ ;; - sin(s)*sin(k*pi/1024)
+ (when (zerop-qd a)
+ (return-from cos-qd +qd-one+))
+
+ (multiple-value-bind (j r)
+ (rem-pi/2-int a)
+ (let* ((j (if (= j 3)
+ -1
+ j))
+ (abs-j (abs j)))
+ (multiple-value-bind (k tmp)
+ (divrem-qd r +qd-pi/1024+)
+ (let* ((k (truncate (qd-0 k)))
+ (abs-k (abs k)))
+ (assert (<= abs-j 2))
+ (assert (<= abs-k 256))
+ ;; Compute sin(s) and cos(s)
+ (multiple-value-bind (sin-t cos-t)
+ (sincos-taylor tmp)
+ (multiple-value-bind (s c)
+ (cond ((zerop abs-k)
+ (values sin-t cos-t))
+ (t
+ ;; Compute sin(s+k*pi/1024), cos(s+k*pi/1024)
+ (let ((u (aref +qd-cos-table+ (cl:1- abs-k)))
+ (v (aref +qd-sin-table+ (cl:1- abs-k))))
+ (cond ((plusp k)
+ ;; sin(s) * cos(k*pi/1024)
+ ;; + cos(s) * sin(k*pi/1024)
+ ;;
+ ;; cos(s) * cos(k*pi/1024)
+ ;; - sin(s) * sin(k*pi/1024)
+ (values (add-qd (mul-qd u sin-t)
+ (mul-qd v cos-t))
+ (sub-qd (mul-qd u cos-t)
+ (mul-qd v sin-t))))
+ (t
+ ;; sin(s) * cos(k*pi/1024)
+ ;; - cos(s) * sin(|k|*pi/1024)
+ ;;
+ ;; cos(s) * cos(k*pi/1024)
+ ;; + sin(s) * sin(|k|*pi/1024)
+ (values (sub-qd (mul-qd u sin-t)
+ (mul-qd v cos-t))
+ (add-qd (mul-qd u cos-t)
+ (mul-qd v sin-t))))))))
+ #+nil
+ (progn
+ (format t "s = ~/qd::qd-format/~%" s)
+ (format t "c = ~/qd::qd-format/~%" c))
+ ;; sin(x) = sin(s+k*pi/1024) * cos(j*pi/2)
+ ;; + cos(s+k*pi/1024) * sin(j*pi/2)
+ (cond ((zerop abs-j)
+ ;; cos(j*pi/2) = 1, sin(j*pi/2) = 0
+ c)
+ ((= j 1)
+ ;; cos(j*pi/2) = 0, sin(j*pi/2) = 1
+ (neg-qd s))
+ ((= j -1)
+ ;; cos(j*pi/2) = 0, sin(j*pi/2) = -1
+ s)
+ (t
+ ;; cos(j*pi/2) = -1, sin(j*pi/2) = 0
+ (neg-qd c))))))))))
+
(defun atan2-qd/newton (y x)
(declare (type %quad-double y x)
#+nil
1
0
Update of /project/oct/cvsroot/oct
In directory clnet:/tmp/cvs-serv23423
Modified Files:
qd.lisp
Log Message:
Redo implementation of INTEGER-DECODE-QD. It used to return way too
many digits if one of the components was 0. This causes problems
because the resulting integer can't even be coerced back to a
quad-double.
--- /project/oct/cvsroot/oct/qd.lisp 2007/09/18 11:20:16 1.54
+++ /project/oct/cvsroot/oct/qd.lisp 2007/10/13 15:34:51 1.55
@@ -1064,6 +1064,7 @@
lo-exp
sign)))))
+#+(or)
(defun integer-decode-qd (q)
(declare (type %quad-double q))
;; Integer decode each component and then create the appropriate
@@ -1093,6 +1094,63 @@
q3-exp
q0-sign))))))
+#+(or)
+(defun integer-decode-qd (q)
+ (declare (type %quad-double q))
+ ;; Integer decode each component and then create the appropriate
+ ;; integer by shifting and adding all the parts together. If any
+ ;; component is zero, we stop.
+ (multiple-value-bind (q0-int q0-exp q0-sign)
+ (integer-decode-float (qd-0 q))
+ (if (zerop (qd-1 q))
+ (values q0-int q0-exp q0-sign)
+ (multiple-value-bind (q1-int q1-exp q1-sign)
+ (integer-decode-float (qd-1 q))
+ (setf q0-int (+ (ash q0-int (- q0-exp q1-exp))
+ (* q1-sign q1-int)))
+ (if (zerop (qd-2 q))
+ (values q0-int q1-exp q0-sign)
+ (multiple-value-bind (q2-int q2-exp q2-sign)
+ (integer-decode-float (qd-2 q))
+ (setf q0-int (+ (ash q0-int (- q1-exp q2-exp))
+ (* q2-sign q2-int)))
+ (if (zerop (qd-3 q))
+ (values q0-int q2-exp q0-sign)
+ (multiple-value-bind (q3-int q3-exp q3-sign)
+ (integer-decode-float (qd-3 q))
+ (values (+ (ash q0-int (- q2-exp q3-exp))
+ (* q3-sign q3-int))
+ q3-exp
+ q0-sign)))))))))
+
+(defun integer-decode-qd (q)
+ (declare (type %quad-double q))
+ ;; Integer decode each component and then create the appropriate
+ ;; integer by shifting and adding all the parts together. If any
+ ;; component is zero, we stop.
+ (multiple-value-bind (q0-int q0-exp q0-sign)
+ (integer-decode-float (qd-0 q))
+ (when (zerop (qd-1 q))
+ (return-from integer-decode-qd (values q0-int q0-exp q0-sign)))
+ (multiple-value-bind (q1-int q1-exp q1-sign)
+ (integer-decode-float (qd-1 q))
+ (setf q0-int (+ (ash q0-int (- q0-exp q1-exp))
+ (* q1-sign q1-int)))
+ (when (zerop (qd-2 q))
+ (return-from integer-decode-qd (values q0-int q1-exp q0-sign)))
+ (multiple-value-bind (q2-int q2-exp q2-sign)
+ (integer-decode-float (qd-2 q))
+ (setf q0-int (+ (ash q0-int (- q1-exp q2-exp))
+ (* q2-sign q2-int)))
+ (when (zerop (qd-3 q))
+ (return-from integer-decode-qd (values q0-int q2-exp q0-sign)))
+ (multiple-value-bind (q3-int q3-exp q3-sign)
+ (integer-decode-float (qd-3 q))
+ (values (+ (ash q0-int (- q2-exp q3-exp))
+ (* q3-sign q3-int))
+ q3-exp
+ q0-sign))))))
+
(declaim (inline scale-float-qd))
#+(or)
(defun scale-float-qd (qd k)
1
0
[oct-cvs] Oct commit: oct branch-test.lisp qd-class.lisp qd-complex.lisp qd-format.lisp qd-methods.lisp rt-tests.lisp tests.lisp timing2.lisp
by rtoy 13 Oct '07
by rtoy 13 Oct '07
13 Oct '07
Update of /project/oct/cvsroot/oct
In directory clnet:/tmp/cvs-serv32041
Modified Files:
branch-test.lisp qd-class.lisp qd-complex.lisp qd-format.lisp
qd-methods.lisp rt-tests.lisp tests.lisp timing2.lisp
Log Message:
Use package #:oct, not #:qd.
--- /project/oct/cvsroot/oct/branch-test.lisp 2007/08/25 18:49:27 1.1
+++ /project/oct/cvsroot/oct/branch-test.lisp 2007/10/13 02:14:43 1.2
@@ -30,7 +30,7 @@
;;; computing the values correctly for the branch cuts. We need to
;;; fix this.
-(in-package #:qd)
+(in-package #:oct)
(defun check-signs (fun arg real-sign imag-sign)
(let* ((z (funcall fun arg))
--- /project/oct/cvsroot/oct/qd-class.lisp 2007/09/19 17:30:04 1.26
+++ /project/oct/cvsroot/oct/qd-class.lisp 2007/10/13 02:14:43 1.27
@@ -23,7 +23,7 @@
;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
;;;; OTHER DEALINGS IN THE SOFTWARE.
-(in-package #:qd)
+(in-package #:oct)
(define-symbol-macro * cl:*)
(define-symbol-macro - cl:-)
--- /project/oct/cvsroot/oct/qd-complex.lisp 2007/09/18 12:46:36 1.37
+++ /project/oct/cvsroot/oct/qd-complex.lisp 2007/10/13 02:14:43 1.38
@@ -23,7 +23,7 @@
;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
;;;; OTHER DEALINGS IN THE SOFTWARE.
-(in-package #:qd)
+(in-package #:oct)
(defmethod add1 ((a qd-complex))
(make-instance 'qd-complex
--- /project/oct/cvsroot/oct/qd-format.lisp 2007/08/29 14:37:42 1.6
+++ /project/oct/cvsroot/oct/qd-format.lisp 2007/10/13 02:14:43 1.7
@@ -23,7 +23,7 @@
;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
;;;; OTHER DEALINGS IN THE SOFTWARE.
-(in-package #:qd)
+(in-package #:oct)
(defun qd-scale-exponent (original-x)
(let* ((x original-x))
--- /project/oct/cvsroot/oct/qd-methods.lisp 2007/10/10 15:24:07 1.62
+++ /project/oct/cvsroot/oct/qd-methods.lisp 2007/10/13 02:14:43 1.63
@@ -23,7 +23,7 @@
;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
;;;; OTHER DEALINGS IN THE SOFTWARE.
-(in-package #:qd)
+(in-package #:oct)
(defconstant +pi+
(make-instance 'qd-real :value qdi:+qd-pi+)
--- /project/oct/cvsroot/oct/rt-tests.lisp 2007/09/18 03:05:56 1.2
+++ /project/oct/cvsroot/oct/rt-tests.lisp 2007/10/13 02:14:43 1.3
@@ -23,7 +23,7 @@
;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
;;;; OTHER DEALINGS IN THE SOFTWARE.
-(in-package #:qd)
+(in-package #:oct)
;; Compute how many bits are the same for two numbers EST and TRUE.
;; Return T if they are identical.
--- /project/oct/cvsroot/oct/tests.lisp 2007/08/25 18:49:27 1.1
+++ /project/oct/cvsroot/oct/tests.lisp 2007/10/13 02:14:43 1.2
@@ -24,7 +24,7 @@
;;;; OTHER DEALINGS IN THE SOFTWARE.
-(in-package #:qd)
+(in-package #:oct)
(defun bit-accuracy (est true)
(let* ((diff (abs (- est true)))
--- /project/oct/cvsroot/oct/timing2.lisp 2007/10/11 18:59:05 1.2
+++ /project/oct/cvsroot/oct/timing2.lisp 2007/10/13 02:14:43 1.3
@@ -27,7 +27,7 @@
;;; test code. I've tried to make these versions time the same
;;; operations as Yozo's.
-(in-package #:qd)
+(in-package #:oct)
(defun time-add (&optional (n 100000))
(declare (fixnum n)
1
0
Update of /project/oct/cvsroot/oct
In directory clnet:/tmp/cvs-serv8913
Modified Files:
timing2.lisp
Log Message:
Correct the timing info. I accidentally ran the C++ version on a
different (slower) machine. The results make some sense now.
--- /project/oct/cvsroot/oct/timing2.lisp 2007/10/11 16:29:08 1.1
+++ /project/oct/cvsroot/oct/timing2.lisp 2007/10/11 18:59:05 1.2
@@ -121,27 +121,25 @@
Some test results. These were all run on a Sun Blade 1500 using a 1.5
GHz Ultrasparc III. I used the default configuration when compiling
-qd, and used Sun's C++ compiler. For the Lisp timing, I used CMUCL.
+qd, and used Sun's C++ compiler with -O. For the Lisp timing, I used
+CMUCL.
Executive summary:
Test Time
qd oct
---- -----------
-add 0.036 0.09
-mul 0.117 0.13
-div 0.388 0.29
-sqrt 0.142 0.11
-sin 0.136 0.14
-log 0.231 0.12
+add 0.023 0.09
+mul 0.075 0.13
+div 0.299 0.29
+sqrt 0.105 0.11
+sin 0.115 0.14
+log 0.194 0.12
Times are in sec for the test. The default number of iterations were
-used. Some of the results are a bit surprising. I was expecting the
-C++ code to be faster, and that's the case for add and mul. But oct
-with CMUCL was faster in div, sqrt, and log. Note that oct uses a
-different algorithm for log than qd, so that could explain the
-difference. Note also that CLOS dispatch is included in these timing
-results. CMUCL's CLOS dispatch is ok, but not great.
+used. Most of the timings match my expectations, including the log
+test. Oct uses a different algorithm (Halley's method) which is
+faster (in Lisp) than the algorithm used in qd (Newtwon iteration).
-------------------------------------------------------------------------------
@@ -149,45 +147,41 @@
The output from qd_timer -qd -v:
-Timing qd_real
---------------
-
Timing addition...
-n = 100000 t = 0.0362288
+n = 100000 t = 0.0231462
b = 1.428571e+04
-100000 operations in 0.0362288 s.
- 0.362288 us
+100000 operations in 0.0231462 s.
+ 0.231462 us
Timing multiplication ...
-n = 100000 t = 0.11686
+n = 100000 t = 0.0749929
b = 2.718268e+00
-100000 operations in 0.11686 s.
- 1.168602 us
+100000 operations in 0.0749929 s.
+ 0.749929 us
Timing division ...
-n = 100000 t = 0.388279
+n = 100000 t = 0.298858
b = 0.367881
-100000 operations in 0.388279 s.
- 3.882788 us
+100000 operations in 0.298858 s.
+ 2.988580 us
Timing square root ...
-n = 10000 t = 0.141866
+n = 10000 t = 0.105049
a = 2.821980
-10000 operations in 0.141866 s.
- 14.186575 us
+10000 operations in 0.105049 s.
+ 10.504860 us
Timing sin ...
-n = 2000 t = 0.136080
+n = 2000 t = 0.114943
a = 3.141593
-2000 operations in 0.136080 s.
- 68.039791 us
+2000 operations in 0.114943 s.
+ 57.471350 us
Timing log ...
-n = 1000 t = 0.230506
+n = 1000 t = 0.193698
a = -50.100000
-1000 operations in 0.230506 s.
-230.506166 us
-
+1000 operations in 0.193698 s.
+193.697800 us
The output from CMUCL:
QD> (time-add)
1
0
Update of /project/oct/cvsroot/oct
In directory clnet:/tmp/cvs-serv32694
Modified Files:
qd-const.lisp qd-extra.lisp
Log Message:
qd-const.lisp:
o Add code for clisp to create the constants we need. (Clisp has
arbitrary length long-float numbers.)
o Update the constants accordingly.
o Move the atan table stuff for CORDIC to qd-extra.lisp, since we
don't need them in oct itself.
qd-extra.lisp:
o atan table stuff for CORDIC moved here.
--- /project/oct/cvsroot/oct/qd-const.lisp 2007/08/27 17:49:19 1.16
+++ /project/oct/cvsroot/oct/qd-const.lisp 2007/10/11 17:47:08 1.17
@@ -31,50 +31,78 @@
(defconstant +qd-one+
(make-qd-d 1d0))
-(defconstant +qd-2pi+
- (%make-qd-d (scale-float (float 7074237752028440 1.0d0) -50)
- (scale-float (float 4967757600021511 1.0d0) -104)
- (scale-float (float -8753721960665020 1.0d0) -160)
- (scale-float (float 5857755168774013 1.0d0) -214)))
-
-;; 3.1415926535897932384626433832795028841971693993751058209749445923078L0
+;; 3.1415926535897932384626433832795028841971693993751058209749445923078L0
+;; #q3.1415926535897932384626433832795028841971693993751058209749445923q0
(defconstant +qd-pi+
- (%make-qd-d (scale-float (float 7074237752028440 1.0d0) -51)
- (scale-float (float 4967757600021511 1.0d0) -105)
- (scale-float (float -8753721960665020 1.0d0) -161)
- (scale-float (float 5857755168774013 1.0d0) -215)))
+ (multiple-value-call #'%make-qd-d
+ (renorm-5 (scale-float (float 7074237752028440 1.0d0) -51)
+ (scale-float (float 4967757600021511 1.0d0) -105)
+ (scale-float (float -8753721960665020 1.0d0) -161)
+ (scale-float (float 5857755168774013 1.0d0) -215)
+ (scale-float (float 5380502254059520 1.0d0) -269))))
+
+;; 6.2831853071795864769252867665590057683943387987502116419498891846156328125724L0
+;; #q6.2831853071795864769252867665590057683943387987502116419498891846q0
+(defconstant +qd-2pi+
+ (multiple-value-call #'%make-qd-d
+ (renorm-5 (scale-float (float 7074237752028440 1.0d0) -50)
+ (scale-float (float 4967757600021511 1.0d0) -104)
+ (scale-float (float -8753721960665020 1.0d0) -160)
+ (scale-float (float 5857755168774013 1.0d0) -214)
+ (scale-float (float 5380502254059520 1.0d0) -268))))
+;; 1.5707963267948966192313216916397514420985846996875529104874722961539082031431L0
+;; #q1.57079632679489661923132169163975144209858469968755291048747229615q0
(defconstant +qd-pi/2+
- (%make-qd-d (scale-float (float 7074237752028440 1.0d0) -52)
- (scale-float (float 4967757600021511 1.0d0) -106)
- (scale-float (float -8753721960665020 1.0d0) -162)
- (scale-float (float 5857755168774013 1.0d0) -216)))
+ (multiple-value-call #'%make-qd-d
+ (renorm-5 (scale-float (float 7074237752028440 1.0d0) -52)
+ (scale-float (float 4967757600021511 1.0d0) -106)
+ (scale-float (float -8753721960665020 1.0d0) -162)
+ (scale-float (float 5857755168774013 1.0d0) -216)
+ (scale-float (float 5380502254059520 1.0d0) -270))))
+;; 0.78539816339744830961566084581987572104929234984377645524373614807695410157155L0
+;; #q0.785398163397448309615660845819875721049292349843776455243736148076q0
(defconstant +qd-pi/4+
- (%make-qd-d (scale-float (float 7074237752028440 1.0d0) -53)
- (scale-float (float 4967757600021511 1.0d0) -107)
- (scale-float (float -8753721960665020 1.0d0) -163)
- (scale-float (float 5857755168774013 1.0d0) -217)))
+ (multiple-value-call #'%make-qd-d
+ (renorm-5 (scale-float (float 7074237752028440 1.0d0) -53)
+ (scale-float (float 4967757600021511 1.0d0) -107)
+ (scale-float (float -8753721960665020 1.0d0) -163)
+ (scale-float (float 5857755168774013 1.0d0) -217)
+ (scale-float (float 5380502254059520 1.0d0) -271))))
+;; 2.35619449019234492884698253745962716314787704953132936573120844423086230471467L0
+;; #q2.35619449019234492884698253745962716314787704953132936573120844423q0
(defconstant +qd-3pi/4+
- (%make-qd-d (scale-float (float 5305678314021330 1.0d0) -51)
- (scale-float (float 7451636400032266 1.0d0) -106)
- (scale-float (float 5724553519491610 1.0d0) -160)
- (scale-float (float -6810541066450736 1.0d0) -214)))
+ (multiple-value-call #'%make-qd-d
+ (renorm-5 (scale-float (float 5305678314021330 1.0d0) -51)
+ (scale-float (float 7451636400032266 1.0d0) -106)
+ (scale-float (float 5724553519491610 1.0d0) -160)
+ (scale-float (float -6810541066450737 1.0d0) -214)
+ (scale-float (float -7491566988951552 1.0d0) -273))))
+;; 0.00306796157577128245943617517898388953534879824157725177829584432842560195926387L0
+;; #q0.00306796157577128245943617517898388953534879824157725177829584432842q0
(defconstant +qd-pi/1024+
- (%make-qd-d (scale-float (float 7074237752028440 1.0d0) -61)
- (scale-float (float 4967757600021511 1.0d0) -115)
- (scale-float (float -8753721960665020 1.0d0) -171)
- (scale-float (float 5857755168774013 1.0d0) -225)))
+ (multiple-value-call #'%make-qd-d
+ (renorm-5 (scale-float (float 7074237752028440 1.0d0) -61)
+ (scale-float (float 4967757600021511 1.0d0) -115)
+ (scale-float (float -8753721960665020 1.0d0) -171)
+ (scale-float (float 5857755168774013 1.0d0) -225)
+ (scale-float (float 5380502254059520 1.0d0) -279))))
+;; 2.71828182845904523536028747135266249775724709369995957496696762772407663035355L0
+;; #q2.71828182845904523536028747135266249775724709369995957496696762773q0
(defconstant +qd-e+
- (%make-qd-d (scale-float (float 6121026514868073 1.0d0) -51)
- (scale-float (float 5864240480059706 1.0d0) -105)
- (scale-float (float -6219324074349538 1.0d0) -161)
- (scale-float (float 7980724272743020 1.0d0) -215)))
+ (multiple-value-call #'%make-qd-d
+ (renorm-5 (scale-float (float 6121026514868073 1.0d0) -51)
+ (scale-float (float 5864240480059706 1.0d0) -105)
+ (scale-float (float -6219324074349538 1.0d0) -161)
+ (scale-float (float 7980724272743021 1.0d0) -215)
+ (scale-float (float -8855251465666560 1.0d0) -269))))
-;; 0.6931471805599453094172321214581765680755001343602552541206800094934L0
+;; 0.693147180559945309417232121458176568075500134360255254120680009493393621969696L0
+;; #q0.693147180559945309417232121458176568075500134360255254120680009495q0
(defconstant +qd-log2+
#+nil
(make-qd-d 6.931471805599452862d-01
@@ -82,2510 +110,3198 @@
5.707708438416212066d-34
-3.582432210601811423d-50)
(%make-qd-d (scale-float (float 6243314768165359 1.0d0) -53)
- (scale-float (float 7525737178955839 1.0d0) -108)
- (scale-float (float 6673460182522164 1.0d0) -163)
- (scale-float (float -7545482916914656 1.0d0) -217)))
+ (scale-float (float 7525737178955839 1.0d0) -108)
+ (scale-float (float 6673460182522164 1.0d0) -163)
+ (scale-float (float -7545482916914641 1.0d0) -217)))
;; The rest of log(2) such that (+ +qd-log2+ +qd-log2-extra+) is
;; log(2) to twice the precision of a quad-double.
(defconstant +qd-log2-extra+
- #+nil
- (make-qd-d (scale-float (float 4141960528156623 1d0) (- -53 212))
- (scale-float (float 3973120087747366 1d0) (- -106 212))
- (scale-float (float 752798645508048 1d0) (- -159 212))
- (scale-float (float 7618435247650241 1d0) (- -212 212)))
- (%make-qd-d (scale-float (float 8283921056313247 1.0d0) -266)
- (scale-float (float -8487672633970079 1.0d0) -322)
- (scale-float (float 6075158146775579 1.0d0) -376)
- (scale-float (float 4764384374407424 1.0d0) -432)))
+ (%make-qd-d (scale-float (float -5130503840205860 1.0d0) -271)
+ (scale-float (float 8312425932334613 1.0d0) -326)
+ (scale-float (float 7130537800999345 1.0d0) -380)
+ (scale-float (float 8345111905183492 1.0d0) -437)))
+;; Log(10)
+;; 2.30258509299404568401799145468436420760110148862877297603332790096757260967737L0
+;; #q2.30258509299404568401799145468436420760110148862877297603332790095q0
(defconstant +qd-log10+
- (%make-qd-d (scale-float (float 5184960683398422 1.0d0) -51)
- (scale-float (float -8805633374462953 1.0d0) -105)
- (scale-float (float -7296007962371596 1.0d0) -159)
- (scale-float (float -5296362421624049 1.0d0) -213)))
+ (multiple-value-call #'%make-qd-d
+ (renorm-5 (scale-float (float 5184960683398422 1.0d0) -51)
+ (scale-float (float -8805633374462953 1.0d0) -105)
+ (scale-float (float -7296007962371596 1.0d0) -159)
+ (scale-float (float -5296362421624049 1.0d0) -213)
+ (scale-float (float 4574234754834432 1.0d0) -267))))
(defconstant +qd-eps+
(scale-float 1d0 -209))
(defconstant +qd-sin-table+
- (make-array 256
- :initial-contents
+ (make-array 256 :initial-contents
(list
- (%make-qd-d (scale-float (float 7074226654454970 1.0d0) -61)
- (scale-float (float 5271335698347442 1.0d0) -115)
- (scale-float (float 7913822574154586 1.0d0) -170)
- (scale-float (float -4804352564403217 1.0d0) -224))
- (%make-qd-d (scale-float (float 7074193361797233 1.0d0) -60)
- (scale-float (float 7522205695703707 1.0d0) -116)
- (scale-float (float 7787047147324868 1.0d0) -175)
- (scale-float (float -6728650631839117 1.0d0) -232))
- (%make-qd-d (scale-float (float 5305603405682435 1.0d0) -59)
- (scale-float (float -5041342953893321 1.0d0) -115)
- (scale-float (float 8335474310793984 1.0d0) -170)
- (scale-float (float 6743025479788607 1.0d0) -225))
- (%make-qd-d (scale-float (float 7074060192106372 1.0d0) -59)
- (scale-float (float 7185921569156509 1.0d0) -113)
- (scale-float (float -6016799991653399 1.0d0) -170)
- (scale-float (float -5577663458559350 1.0d0) -224))
- (%make-qd-d (scale-float (float 8842450394781643 1.0d0) -59)
- (scale-float (float -8771095658131507 1.0d0) -113)
- (scale-float (float 8705573170983202 1.0d0) -167)
- (scale-float (float -4692271881407162 1.0d0) -221))
- (%make-qd-d (scale-float (float 5305378684473085 1.0d0) -58)
- (scale-float (float 7704904742673764 1.0d0) -113)
- (scale-float (float 7308626221243388 1.0d0) -167)
- (scale-float (float 4905810005485122 1.0d0) -223))
- (%make-qd-d (scale-float (float 6189482235310630 1.0d0) -58)
- (scale-float (float -4715432182777116 1.0d0) -113)
- (scale-float (float -8243063834437948 1.0d0) -168)
- (scale-float (float 8566398131092182 1.0d0) -222))
- (%make-qd-d (scale-float (float 7073527528384126 1.0d0) -58)
- (scale-float (float -7632135562854704 1.0d0) -116)
- (scale-float (float 7289253486557291 1.0d0) -170)
- (scale-float (float -7590289326531904 1.0d0) -224))
- (%make-qd-d (scale-float (float 7957506242722589 1.0d0) -58)
- (scale-float (float -8272553670927212 1.0d0) -112)
- (scale-float (float -6591995918745221 1.0d0) -166)
- (scale-float (float -4636857096009409 1.0d0) -220))
- (%make-qd-d (scale-float (float 8841410057981697 1.0d0) -58)
- (scale-float (float -5627969559234246 1.0d0) -118)
- (scale-float (float -5998146560009404 1.0d0) -171)
- (scale-float (float -8771421402740511 1.0d0) -225))
- (%make-qd-d (scale-float (float 4862615327261055 1.0d0) -57)
- (scale-float (float -5217239170857332 1.0d0) -111)
- (scale-float (float -6336246586437174 1.0d0) -165)
- (scale-float (float 5522987872330073 1.0d0) -219))
- (%make-qd-d (scale-float (float 5304479856743885 1.0d0) -57)
- (scale-float (float 6340842145528512 1.0d0) -113)
- (scale-float (float -7566831660699656 1.0d0) -167)
- (scale-float (float -7114412056121165 1.0d0) -221))
- (%make-qd-d (scale-float (float 5746294458442105 1.0d0) -57)
- (scale-float (float 4845186982751657 1.0d0) -113)
- (scale-float (float 6382779118838290 1.0d0) -167)
- (scale-float (float 8090523197847278 1.0d0) -221))
- (%make-qd-d (scale-float (float 6188054973828419 1.0d0) -57)
- (scale-float (float 7360584602593505 1.0d0) -111)
- (scale-float (float 7946504576809704 1.0d0) -165)
- (scale-float (float 5731224328660037 1.0d0) -219))
- (%make-qd-d (scale-float (float 6629757244884614 1.0d0) -57)
- (scale-float (float -5806448969106076 1.0d0) -112)
- (scale-float (float 7037190739826080 1.0d0) -166)
- (scale-float (float 5552037724818900 1.0d0) -223))
- (%make-qd-d (scale-float (float 7071397114140692 1.0d0) -57)
- (scale-float (float -7057477599664784 1.0d0) -113)
- (scale-float (float -8290828586102809 1.0d0) -167)
- (scale-float (float -6697971486035361 1.0d0) -222))
- (%make-qd-d (scale-float (float 7512970424714007 1.0d0) -57)
- (scale-float (float -6294023145184360 1.0d0) -111)
- (scale-float (float -6395720725273163 1.0d0) -165)
- (scale-float (float -7071774105018818 1.0d0) -219))
- (%make-qd-d (scale-float (float 7954473020348387 1.0d0) -57)
- (scale-float (float -6926679705961747 1.0d0) -112)
- (scale-float (float -6427721451859560 1.0d0) -167)
- (scale-float (float 8028310015339695 1.0d0) -222))
- (%make-qd-d (scale-float (float 8395900745453257 1.0d0) -57)
- (scale-float (float 4839201044457661 1.0d0) -114)
- (scale-float (float 5799570435052711 1.0d0) -171)
- (scale-float (float -5529353034270453 1.0d0) -226))
- (%make-qd-d (scale-float (float 8837249445142752 1.0d0) -57)
- (scale-float (float -5314952832401406 1.0d0) -113)
- (scale-float (float -7992910057771694 1.0d0) -167)
- (scale-float (float 8886165424019995 1.0d0) -221))
- (%make-qd-d (scale-float (float 4639257482637412 1.0d0) -56)
- (scale-float (float -5494228531443205 1.0d0) -110)
- (scale-float (float 6222043711193090 1.0d0) -167)
- (scale-float (float 4966112393074763 1.0d0) -221))
- (%make-qd-d (scale-float (float 4859846576245171 1.0d0) -56)
- (scale-float (float -8985502920124077 1.0d0) -110)
- (scale-float (float 7440463379940699 1.0d0) -165)
- (scale-float (float -6641325637598524 1.0d0) -219))
- (%make-qd-d (scale-float (float 5080389927126093 1.0d0) -56)
- (scale-float (float -8898661047761268 1.0d0) -110)
- (scale-float (float -4667727006484474 1.0d0) -164)
- (scale-float (float 8478536563323652 1.0d0) -218))
- (%make-qd-d (scale-float (float 5300885459442166 1.0d0) -56)
- (scale-float (float -7213383224879916 1.0d0) -111)
- (scale-float (float -8534269322822802 1.0d0) -166)
- (scale-float (float -6601733372219099 1.0d0) -221))
- (%make-qd-d (scale-float (float 5521331097805465 1.0d0) -56)
- (scale-float (float 4829604598179156 1.0d0) -114)
- (scale-float (float -7893599535446481 1.0d0) -172)
- (scale-float (float -5384211130874634 1.0d0) -226))
- (%make-qd-d (scale-float (float 5741724767297686 1.0d0) -56)
- (scale-float (float -5824155811343436 1.0d0) -110)
- (scale-float (float 6673985610265747 1.0d0) -164)
- (scale-float (float 4800782739249759 1.0d0) -220))
- (%make-qd-d (scale-float (float 5962064393489674 1.0d0) -56)
- (scale-float (float 7651360099479755 1.0d0) -112)
- (scale-float (float 6974817744235260 1.0d0) -167)
- (scale-float (float 7824970858345265 1.0d0) -224))
- (%make-qd-d (scale-float (float 6182347902460953 1.0d0) -56)
- (scale-float (float -8796242544913562 1.0d0) -111)
- (scale-float (float -7546272905440146 1.0d0) -165)
- (scale-float (float 6512061516709508 1.0d0) -219))
- (%make-qd-d (scale-float (float 6402573220819241 1.0d0) -56)
- (scale-float (float -4868008820827392 1.0d0) -110)
- (scale-float (float 8779206749884591 1.0d0) -164)
- (scale-float (float -8648470131210884 1.0d0) -218))
- (%make-qd-d (scale-float (float 6622738275719969 1.0d0) -56)
- (scale-float (float 6182934508221337 1.0d0) -110)
- (scale-float (float 7353282250945404 1.0d0) -165)
- (scale-float (float -8166192353391047 1.0d0) -220))
- (%make-qd-d (scale-float (float 6842840994885793 1.0d0) -56)
- (scale-float (float -8552477024466766 1.0d0) -110)
- (scale-float (float -4997427595980002 1.0d0) -164)
- (scale-float (float 5821398332046138 1.0d0) -218))
- (%make-qd-d (scale-float (float 7062879306626092 1.0d0) -56)
- (scale-float (float -8487236864497288 1.0d0) -112)
- (scale-float (float -4942162982072151 1.0d0) -168)
- (scale-float (float -4811145377091453 1.0d0) -223))
- (%make-qd-d (scale-float (float 7282851139856476 1.0d0) -56)
- (scale-float (float 8609951212389606 1.0d0) -111)
- (scale-float (float -5614142183842944 1.0d0) -165)
- (scale-float (float 6068416796043208 1.0d0) -219))
- (%make-qd-d (scale-float (float 7502754424118275 1.0d0) -56)
- (scale-float (float 8536169017599928 1.0d0) -110)
- (scale-float (float 7981550951145323 1.0d0) -165)
- (scale-float (float -8340622458345952 1.0d0) -220))
- (%make-qd-d (scale-float (float 7722587089598028 1.0d0) -56)
- (scale-float (float 8362719068102409 1.0d0) -110)
- (scale-float (float -7853008551658302 1.0d0) -166)
- (scale-float (float -6852498985157079 1.0d0) -220))
- (%make-qd-d (scale-float (float 7942347067146965 1.0d0) -56)
- (scale-float (float -5897359216343841 1.0d0) -113)
- (scale-float (float 7767443765766550 1.0d0) -169)
- (scale-float (float 8203371449514667 1.0d0) -225))
- (%make-qd-d (scale-float (float 8162032288300481 1.0d0) -56)
- (scale-float (float 7035687121204020 1.0d0) -111)
- (scale-float (float 5733124558862890 1.0d0) -168)
- (scale-float (float 6681272325774503 1.0d0) -223))
- (%make-qd-d (scale-float (float 8381640685297609 1.0d0) -56)
- (scale-float (float 5345425417737710 1.0d0) -112)
- (scale-float (float -8790124191829010 1.0d0) -166)
- (scale-float (float 7297964455931531 1.0d0) -225))
- (%make-qd-d (scale-float (float 8601170191100479 1.0d0) -56)
- (scale-float (float -5127404751534987 1.0d0) -110)
- (scale-float (float 8258618371098839 1.0d0) -164)
- (scale-float (float 6354475203379338 1.0d0) -221))
- (%make-qd-d (scale-float (float 8820618739413774 1.0d0) -56)
- (scale-float (float 7361249450583405 1.0d0) -111)
- (scale-float (float 8489175872343797 1.0d0) -165)
- (scale-float (float -4838736209868180 1.0d0) -220))
- (%make-qd-d (scale-float (float 4519992132352091 1.0d0) -55)
- (scale-float (float 6319901705801489 1.0d0) -110)
- (scale-float (float 5600764727998866 1.0d0) -166)
- (scale-float (float -5650573825026412 1.0d0) -220))
- (%make-qd-d (scale-float (float 4629632351109917 1.0d0) -55)
- (scale-float (float 4958462282897610 1.0d0) -110)
- (scale-float (float -8711719640720475 1.0d0) -165)
- (scale-float (float -8191834096705166 1.0d0) -221))
- (%make-qd-d (scale-float (float 4739228994004870 1.0d0) -55)
- (scale-float (float -6495525500297301 1.0d0) -110)
- (scale-float (float -6071898187587615 1.0d0) -165)
- (scale-float (float -7897922259441451 1.0d0) -220))
- (%make-qd-d (scale-float (float 4848781029471607 1.0d0) -55)
- (scale-float (float -5949746474335482 1.0d0) -109)
- (scale-float (float 7443833788929433 1.0d0) -165)
- (scale-float (float 6779921506403360 1.0d0) -220))
- (%make-qd-d (scale-float (float 4958287426364647 1.0d0) -55)
- (scale-float (float 8600164307618932 1.0d0) -110)
- (scale-float (float -5552907197025195 1.0d0) -164)
- (scale-float (float -6315628760329704 1.0d0) -222))
- (%make-qd-d (scale-float (float 5067747153968079 1.0d0) -55)
- (scale-float (float -5139980116898898 1.0d0) -109)
- (scale-float (float 7128908501078503 1.0d0) -163)
- (scale-float (float 5174334784143035 1.0d0) -217))
- (%make-qd-d (scale-float (float 5177159182005257 1.0d0) -55)
- (scale-float (float 7446222959753664 1.0d0) -109)
- (scale-float (float -4853116546479197 1.0d0) -166)
- (scale-float (float -7114171320789543 1.0d0) -220))
- (%make-qd-d (scale-float (float 5286522480648506 1.0d0) -55)
- (scale-float (float 4837853990883808 1.0d0) -110)
- (scale-float (float 8734494569006220 1.0d0) -164)
- (scale-float (float -6888573428631769 1.0d0) -220))
- (%make-qd-d (scale-float (float 5395836020528807 1.0d0) -55)
- (scale-float (float 5245006079192873 1.0d0) -109)
- (scale-float (float 6070092190452118 1.0d0) -165)
- (scale-float (float 8372804672547503 1.0d0) -219))
- (%make-qd-d (scale-float (float 5505098772745492 1.0d0) -55)
- (scale-float (float -4953034127800088 1.0d0) -109)
- (scale-float (float 6748009167290918 1.0d0) -163)
- (scale-float (float -7947104079128991 1.0d0) -217))
- (%make-qd-d (scale-float (float 5614309708875923 1.0d0) -55)
- (scale-float (float 7879649848150358 1.0d0) -111)
- (scale-float (float -5133702133285397 1.0d0) -165)
- (scale-float (float 6633380945160774 1.0d0) -219))
- (%make-qd-d (scale-float (float 5723467800985178 1.0d0) -55)
- (scale-float (float -5213481504208817 1.0d0) -110)
- (scale-float (float -8683946243665639 1.0d0) -166)
- (scale-float (float -4813600752545885 1.0d0) -220))
- (%make-qd-d (scale-float (float 5832572021635720 1.0d0) -55)
- (scale-float (float 7691426989666512 1.0d0) -109)
- (scale-float (float -5897567218587937 1.0d0) -163)
- (scale-float (float 6426162084210735 1.0d0) -217))
- (%make-qd-d (scale-float (float 5941621343897074 1.0d0) -55)
- (scale-float (float -7311303147276965 1.0d0) -113)
- (scale-float (float 6212875260931578 1.0d0) -167)
- (scale-float (float 5841094814764128 1.0d0) -222))
- (%make-qd-d (scale-float (float 6050614741355486 1.0d0) -55)
- (scale-float (float 7046535347736856 1.0d0) -110)
- (scale-float (float -7778307984467600 1.0d0) -164)
- (scale-float (float -7808429293088315 1.0d0) -218))
- (%make-qd-d (scale-float (float 6159551188123590 1.0d0) -55)
- (scale-float (float 5965947804179142 1.0d0) -109)
- (scale-float (float -7914176807889465 1.0d0) -163)
- (scale-float (float -8743204545259556 1.0d0) -219))
- (%make-qd-d (scale-float (float 6268429658850061 1.0d0) -55)
- (scale-float (float 7548560474328400 1.0d0) -110)
- (scale-float (float -7834123143654772 1.0d0) -165)
- (scale-float (float -5643244224302487 1.0d0) -219))
- (%make-qd-d (scale-float (float 6377249128729266 1.0d0) -55)
- (scale-float (float 8739842904414420 1.0d0) -110)
- (scale-float (float 6557520883328920 1.0d0) -164)
- (scale-float (float 6198578405040918 1.0d0) -220))
- (%make-qd-d (scale-float (float 6486008573510911 1.0d0) -55)
- (scale-float (float 5172944262567044 1.0d0) -109)
- (scale-float (float -8273960648700810 1.0d0) -163)
- (scale-float (float 8265967084369073 1.0d0) -218))
[5397 lines skipped]
--- /project/oct/cvsroot/oct/qd-extra.lisp 2007/08/25 17:08:48 1.2
+++ /project/oct/cvsroot/oct/qd-extra.lisp 2007/10/11 17:47:08 1.3
@@ -477,6 +477,416 @@
(cl:* s y))))
||#
+#||
+;; Here is a function for clisp that can be used to create the atan2 table
+;; that we need.
+
+(defun make-atan-table-data ()
+ (let ((scale 1l0))
+ (dotimes (k 67)
+ (let* ((x (scale-float 1L0 (- 2 k)))
+ (p (atan x)))
+ (setf scale (* scale (cos p)))
+ (multiple-value-bind (int exp sign)
+ (integer-decode-float p)
+ (let* ((len (integer-length int))
+ (wanted (ldb (byte 212 (- len 212)) int))
+ (bit (ldb (byte 1 (- len (* 4 53) 1)) int))
+ (roundp (not (zerop (ldb (byte (- len (* 4 53) 2) 0) int)))))
+ ;;(format t "~&~v,'0b~%" len int)
+ ;;(format t "~b~a~%" wanted (make-string (- len 212) :initial-element #\-))
+ ;;(format t "~v,'-b~%" len (ash bit (- len 212 1)))
+ ;;(format t "~v,'-b~%" len (ldb (byte (- len (* 4 53) 2) 0) int))
+ ;; See if we need to round up the answer.
+ (when (= bit 1)
+ ;; Round to even
+ (cond (roundp
+ (incf wanted))
+ (t
+ ;; Round to even
+ (when (oddp wanted)
+ (incf wanted)))))
+ ;;(format t "~b~a~%" wanted (make-string (- len 212) :initial-element #\-))
+
+ (let* ((i0 (ldb (byte 53 (* 3 53)) wanted))
+ (i1 (ldb (byte 53 (* 2 53)) wanted))
+ (i2 (ldb (byte 53 (* 1 53)) wanted))
+ (i3 (ldb (byte 53 0) wanted)))
+ (write `(make-qd-d
+ (scale-float (float ,i0 1d0) ,(+ exp (- len (* 1 53))))
+ (scale-float (float ,i1 1d0) ,(+ exp (- len (* 2 53))))
+ (scale-float (float ,i2 1d0) ,(+ exp (- len (* 3 53))))
+ (scale-float (float ,i3 1d0) ,(+ exp (- len (* 4 53)))))
+ :case :downcase))))))
+ scale))
+||#
+
+
+#+nil
+(defconstant +atan-table+
+ (make-array 66
+ :initial-contents
+ (list
+ +qd-pi/4+
+ +qd-pi/4+
+ +qd-pi/4+
+ ;; Do we need to make these values more accurate? (The
+ ;; reader has quite a bit of roundoff.)
+ #.(qd-from-string "0.46364760900080611621425623146121440202853705428612026381093308872018q0")
+ #.(qd-from-string "0.24497866312686415417208248121127581091414409838118406712737591466738q0")
+ #.(qd-from-string "0.12435499454676143503135484916387102557317019176980408991511411911572q0")
+ #.(qd-from-string "0.062418809995957348473979112985505113606273887797499194607527816898697q0")
+ #.(qd-from-string "0.031239833430268276253711744892490977032495663725400040255315586255793q0")
+ #.(qd-from-string "0.0156237286204768308028015212565703189111141398009054178814105073966645q0")
+ #.(qd-from-string "0.0078123410601011112964633918421992816212228117250147235574539022483893q0")
+ #.(qd-from-string "0.003906230131966971827628665311424387140357490115202856215213095149011q0")
+ #.(qd-from-string "0.00195312251647881868512148262507671393161074677723351033905753396043094q0")
+ #.(qd-from-string "9.7656218955931943040343019971729085163419701581008759004900725226773q-4")
+ #.(qd-from-string "4.8828121119489827546923962564484866619236113313500303710940335348752q-4")
+ #.(qd-from-string "2.4414062014936176401672294325965998621241779097061761180790046091019q-4")
+ #.(qd-from-string "1.22070311893670204239058646117956300930829409015787498451939837846645q-4")
+ #.(qd-from-string "6.1035156174208775021662569173829153785143536833346179337671134316588q-5")
+ #.(qd-from-string "3.0517578115526096861825953438536019750949675119437837531021156883611q-5")
+ #.(qd-from-string "1.5258789061315762107231935812697885137429238144575874846241186407446q-5")
+ #.(qd-from-string "7.6293945311019702633884823401050905863507439184680771577638306965336q-6")
+ #.(qd-from-string "3.8146972656064962829230756163729937228052573039688663101874392503939q-6")
+ #.(qd-from-string "1.9073486328101870353653693059172441687143421654501533666700577234671q-6")
+ #.(qd-from-string "9.53674316405960879420670689923112390019634124498790160133611802076q-7")
+ #.(qd-from-string "4.7683715820308885992758382144924707587049404378664196740053215887142q-7")
+ #.(qd-from-string "2.3841857910155798249094797721893269783096898769063155913766911372218q-7")
+ #.(qd-from-string "1.19209289550780685311368497137922112645967587664586735576738225215437q-7")
+ #.(qd-from-string "5.9604644775390554413921062141788874250030195782366297314294565710003q-8")
+ #.(qd-from-string "2.9802322387695303676740132767709503349043907067445107249258477840843q-8")
+ #.(qd-from-string "1.4901161193847655147092516595963247108248930025964720012170057805491q-8")
+ #.(qd-from-string "7.4505805969238279871365645744953921132066925545665870075947601416172q-9")
+ #.(qd-from-string "3.725290298461914045267070571811923583671948328737040524231998269239q-9")
+ #.(qd-from-string "1.8626451492309570290958838214764904345065282835738863513491050124951q-9")
+ #.(qd-from-string "9.3132257461547851535573547768456130389292649614929067394376854242196q-10")
+ #.(qd-from-string "4.6566128730773925777884193471057016297347863891561617421323492554414q-10")
+ #.(qd-from-string "2.32830643653869628902042741838821270371274293204981860525486662280605q-10")
+ #.(qd-from-string "1.16415321826934814452599092729852658796396457380014290026584979170883q-10")
+ #.(qd-from-string "5.8207660913467407226496761591231582349549156257795272423976206167147q-11")
+ #.(qd-from-string "2.9103830456733703613273032698903947793693632003639830495829934525029q-11")
+ #.(qd-from-string "1.4551915228366851806639597837362993474211703608936710732067270213307q-11")
+ #.(qd-from-string "7.2759576141834259033201841046703741842764629388821429640111752890838q-12")
+ #.(qd-from-string "3.6379788070917129516601402005837967730345578669779258118296083646486q-12")
+ #.(qd-from-string "1.81898940354585647583007611882297459662931973336029253714520765350336q-12")
+ #.(qd-from-string "9.094947017729282379150388117278718245786649666696631862264792881855q-13")
+ #.(qd-from-string "4.5474735088646411895751949990348397807233312083369623012466392138249q-13")
+ #.(qd-from-string "2.2737367544323205947875976170668549725904164010421166413578155299654q-13")
+ #.(qd-from-string "1.1368683772161602973937988232271068715738020501302644662229139921281q-13")
+ #.(qd-from-string "5.6843418860808014869689941345026335894672525626628305471702634435609q-14")
+ #.(qd-from-string "2.8421709430404007434844970695472041986834065703328538172835210852389q-14")
+ #.(qd-from-string "1.42108547152020037174224853506058802483542582129160672712566632799217q-14")
+ #.(qd-from-string "7.1054273576010018587112426756616725310442822766145084088962160950957q-15")
+ #.(qd-from-string "3.5527136788005009293556213378756778163805352845768135511116874239215q-15")
+ #.(qd-from-string "1.7763568394002504646778106689434441020475669105721016938889503158663q-15")
+ #.(qd-from-string "8.881784197001252323389053344724227002559458638215127117361184578544q-16")
+ #.(qd-from-string "4.440892098500626161694526672362989312819932329776890889670147968684q-16")
+ #.(qd-from-string "2.22044604925031308084726333618160413285249154122211136120876849284695q-16")
+ #.(qd-from-string "1.11022302462515654042363166809081575098156144265276392015109606150467q-16")
+ #.(qd-from-string "5.5511151231257827021181583404540958606019518033159549001888700768492q-17")
+ #.(qd-from-string "2.7755575615628913510590791702270500685127439754144943625236087596052q-17")
+ #.(qd-from-string "1.3877787807814456755295395851135253015328429969268117953154510949506q-17")
+ #.(qd-from-string "6.9388939039072283776476979255676268417598037461585147441443138686883q-18")
+ #.(qd-from-string "3.4694469519536141888238489627838134626418504682698143430180392335861q-18")
+ #.(qd-from-string "1.7347234759768070944119244813919067365411688085337267928772549041983q-18")
+ #.(qd-from-string "8.673617379884035472059622406959533689231148510667158491096568630248q-19")
+ #.(qd-from-string "4.336808689942017736029811203479766845431237313833394811387071078781q-19")
+ #.(qd-from-string "2.16840434497100886801490560173988342281757653922917435142338388484765q-19")
+ #.(qd-from-string "1.08420217248550443400745280086994171142153300490364679392792298560597q-19")
+
+ ))
+ "Table of atan(2^(-k)) for k = 1 to 64. But the first three entries are 1")
+
+(defconstant +atan-table+
+ (make-array 67
+ :initial-contents
+ (list
+ (%make-qd-d (scale-float (float 5970951936056572 1.0d0) -52)
+ (scale-float (float 5427585433121543 1.0d0) -105)
+ (scale-float (float 5608515294538868 1.0d0) -158)
+ (scale-float (float 445395631680583 1.0d0) -211))
+ (%make-qd-d (scale-float (float 4986154552901188 1.0d0) -52)
+ (scale-float (float 3814906810089799 1.0d0) -105)
+ (scale-float (float 1896417689773139 1.0d0) -158)
+ (scale-float (float 3393132800284032 1.0d0) -211))
+ (%make-qd-d (scale-float (float 7074237752028440 1.0d0) -53)
+ (scale-float (float 2483878800010755 1.0d0) -106)
+ (scale-float (float 3956492004828932 1.0d0) -159)
+ (scale-float (float 2434854662709436 1.0d0) -212))
+ (%make-qd-d (scale-float (float 8352332796509007 1.0d0) -54)
+ (scale-float (float 3683087214424816 1.0d0) -107)
+ (scale-float (float 8240297260223171 1.0d0) -160)
+ (scale-float (float 5174086704442609 1.0d0) -213))
+ (%make-qd-d (scale-float (float 8826286527774941 1.0d0) -55)
+ (scale-float (float 3471944699336670 1.0d0) -108)
+ (scale-float (float 4798212191802497 1.0d0) -161)
+ (scale-float (float 6908472993489831 1.0d0) -214))
+ (%make-qd-d (scale-float (float 8960721713639277 1.0d0) -56)
+ (scale-float (float 6978747913895162 1.0d0) -109)
+ (scale-float (float 1204496828771308 1.0d0) -162)
+ (scale-float (float 6150314016033077 1.0d0) -215))
+ (%make-qd-d (scale-float (float 8995498542038505 1.0d0) -57)
+ (scale-float (float 6996384121843768 1.0d0) -110)
+ (scale-float (float 6481245652257127 1.0d0) -163)
+ (scale-float (float 6083920726820778 1.0d0) -216))
+ (%make-qd-d (scale-float (float 9004268940523044 1.0d0) -58)
+ (scale-float (float 5921825575778154 1.0d0) -111)
+ (scale-float (float 1742767809528138 1.0d0) -164)
+ (scale-float (float 3392785816514584 1.0d0) -217))
+ (%make-qd-d (scale-float (float 9006466354344602 1.0d0) -59)
+ (scale-float (float 6455912199422039 1.0d0) -112)
+ (scale-float (float 7793493312778976 1.0d0) -165)
+ (scale-float (float 4748718880757240 1.0d0) -218))
+ (%make-qd-d (scale-float (float 9007016009513623 1.0d0) -60)
+ (scale-float (float 1583402193514233 1.0d0) -113)
+ (scale-float (float 4599960241393675 1.0d0) -166)
+ (scale-float (float 4964226307734805 1.0d0) -219))
+ (%make-qd-d (scale-float (float 9007153442175927 1.0d0) -61)
+ (scale-float (float 1458797116501429 1.0d0) -114)
+ (scale-float (float 2180379843517813 1.0d0) -167)
+ (scale-float (float 7244224576758923 1.0d0) -220))
+ (%make-qd-d (scale-float (float 9007187801521083 1.0d0) -62)
+ (scale-float (float 5961909987006481 1.0d0) -115)
+ (scale-float (float 1439161705865198 1.0d0) -168)
+ (scale-float (float 1250151122136839 1.0d0) -221))
+ (%make-qd-d (scale-float (float 9007196391431099 1.0d0) -63)
+ (scale-float (float 6595226783193595 1.0d0) -116)
+ (scale-float (float 7270788700276565 1.0d0) -169)
+ (scale-float (float 5212528258452836 1.0d0) -222))
+ (%make-qd-d (scale-float (float 9007198538913211 1.0d0) -64)
+ (scale-float (float 6605122380416172 1.0d0) -117)
+ (scale-float (float 2579496809882929 1.0d0) -170)
+ (scale-float (float 2545695100421145 1.0d0) -223))
+ (%make-qd-d (scale-float (float 9007199075784027 1.0d0) -65)
+ (scale-float (float 6605276999209814 1.0d0) -118)
+ (scale-float (float 8635423593413256 1.0d0) -171)
+ (scale-float (float 6747877897971029 1.0d0) -224))
+ (%make-qd-d (scale-float (float 9007199210001749 1.0d0) -66)
+ (scale-float (float 6605279415128805 1.0d0) -119)
+ (scale-float (float 5633073770825222 1.0d0) -172)
+ (scale-float (float 744251135568860 1.0d0) -225))
+ (%make-qd-d (scale-float (float 9007199243556181 1.0d0) -67)
+ (scale-float (float 3227579732349669 1.0d0) -120)
+ (scale-float (float 1645511649516378 1.0d0) -173)
+ (scale-float (float 7212311609477561 1.0d0) -226))
+ (%make-qd-d (scale-float (float 9007199251944789 1.0d0) -68)
+ (scale-float (float 3016473500406501 1.0d0) -121)
+ (scale-float (float 1629935234837168 1.0d0) -174)
+ (scale-float (float 1206159191623029 1.0d0) -227))
+ (%make-qd-d (scale-float (float 9007199254041941 1.0d0) -69)
+ (scale-float (float 3003279360882405 1.0d0) -122)
+ (scale-float (float 1629874389467187 1.0d0) -175)
+ (scale-float (float 8712158240272416 1.0d0) -228))
+ (%make-qd-d (scale-float (float 9007199254566229 1.0d0) -70)
+ (scale-float (float 3002454727161717 1.0d0) -123)
+ (scale-float (float 1629874151789961 1.0d0) -176)
+ (scale-float (float 3116377062563786 1.0d0) -229))
+ (%make-qd-d (scale-float (float 9007199254697301 1.0d0) -71)
+ (scale-float (float 3002403187554167 1.0d0) -124)
+ (scale-float (float 3881673964546782 1.0d0) -177)
+ (scale-float (float 6119176246102625 1.0d0) -230))
+ (%make-qd-d (scale-float (float 9007199254730069 1.0d0) -72)
+ (scale-float (float 3002399966328695 1.0d0) -125)
+ (scale-float (float 4198333313342644 1.0d0) -178)
+ (scale-float (float 114377133012236 1.0d0) -231))
+ (%make-qd-d (scale-float (float 9007199254738261 1.0d0) -73)
+ (scale-float (float 3002399765002103 1.0d0) -126)
+ (scale-float (float 4203281115667621 1.0d0) -179)
+ (scale-float (float 7620376512343991 1.0d0) -232))
+ (%make-qd-d (scale-float (float 9007199254740309 1.0d0) -74)
+ (scale-float (float 3002399752419191 1.0d0) -127)
+ (scale-float (float 4203358425078949 1.0d0) -180)
+ (scale-float (float 7121931241085909 1.0d0) -233))
+ (%make-qd-d (scale-float (float 9007199254740821 1.0d0) -75)
+ (scale-float (float 3002399751632759 1.0d0) -128)
+ (scale-float (float 4203359633038501 1.0d0) -181)
+ (scale-float (float 7119984189245056 1.0d0) -234))
+ (%make-qd-d (scale-float (float 9007199254740949 1.0d0) -76)
+ (scale-float (float 3002399751583607 1.0d0) -129)
+ (scale-float (float 4203359651912869 1.0d0) -182)
+ (scale-float (float 7119976583573803 1.0d0) -235))
+ (%make-qd-d (scale-float (float 9007199254740981 1.0d0) -77)
+ (scale-float (float 3002399751580535 1.0d0) -130)
+ (scale-float (float 4203359652207781 1.0d0) -183)
+ (scale-float (float 7119976553864150 1.0d0) -236))
+ (%make-qd-d (scale-float (float 9007199254740989 1.0d0) -78)
+ (scale-float (float 3002399751580343 1.0d0) -131)
+ (scale-float (float 4203359652212389 1.0d0) -184)
+ (scale-float (float 7119976553748096 1.0d0) -237))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -79)
+ (scale-float (float 3002399751580331 1.0d0) -132)
+ (scale-float (float 4203359652212461 1.0d0) -185)
+ (scale-float (float 7119976553747643 1.0d0) -238))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -80)
+ (scale-float (float 7505999378950826 1.0d0) -133)
+ (scale-float (float 6455159465897710 1.0d0) -186)
+ (scale-float (float 8245876460590265 1.0d0) -239))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -81)
+ (scale-float (float 8631899285793450 1.0d0) -134)
+ (scale-float (float 6032947000831726 1.0d0) -187)
+ (scale-float (float 8404206134990009 1.0d0) -240))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -82)
+ (scale-float (float 8913374262504106 1.0d0) -135)
+ (scale-float (float 6006558721765102 1.0d0) -188)
+ (scale-float (float 8406680036152505 1.0d0) -241))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -83)
+ (scale-float (float 8983743006681770 1.0d0) -136)
+ (scale-float (float 6004909454323438 1.0d0) -189)
+ (scale-float (float 8406718690858169 1.0d0) -242))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -84)
+ (scale-float (float 9001335192726186 1.0d0) -137)
+ (scale-float (float 6004806375108334 1.0d0) -190)
+ (scale-float (float 8406719294837945 1.0d0) -243))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -85)
+ (scale-float (float 9005733239237290 1.0d0) -138)
+ (scale-float (float 6004799932657390 1.0d0) -191)
+ (scale-float (float 8406719304275129 1.0d0) -244))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -86)
+ (scale-float (float 9006832750865066 1.0d0) -139)
+ (scale-float (float 6004799530004206 1.0d0) -192)
+ (scale-float (float 8406719304422585 1.0d0) -245))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -87)
+ (scale-float (float 9007107628772010 1.0d0) -140)
+ (scale-float (float 6004799504838382 1.0d0) -193)
+ (scale-float (float 8406719304424889 1.0d0) -246))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -88)
+ (scale-float (float 9007176348248746 1.0d0) -141)
+ (scale-float (float 6004799503265518 1.0d0) -194)
+ (scale-float (float 8406719304424925 1.0d0) -247))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -89)
+ (scale-float (float 9007193528117930 1.0d0) -142)
+ (scale-float (float 6004799503167214 1.0d0) -195)
+ (scale-float (float 8406719304424926 1.0d0) -248))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -90)
+ (scale-float (float 9007197823085226 1.0d0) -143)
+ (scale-float (float 6004799503161070 1.0d0) -196)
+ (scale-float (float 8406719304424926 1.0d0) -249))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -91)
+ (scale-float (float 9007198896827050 1.0d0) -144)
+ (scale-float (float 6004799503160686 1.0d0) -197)
+ (scale-float (float 8406719304424926 1.0d0) -250))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -92)
+ (scale-float (float 9007199165262506 1.0d0) -145)
+ (scale-float (float 6004799503160662 1.0d0) -198)
+ (scale-float (float 8406719304424926 1.0d0) -251))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -93)
+ (scale-float (float 9007199232371370 1.0d0) -146)
+ (scale-float (float 6004799503160661 1.0d0) -199)
+ (scale-float (float 3903119677054430 1.0d0) -252))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -94)
+ (scale-float (float 9007199249148586 1.0d0) -147)
+ (scale-float (float 6004799503160661 1.0d0) -200)
+ (scale-float (float 3058694746922462 1.0d0) -253))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -95)
+ (scale-float (float 9007199253342890 1.0d0) -148)
+ (scale-float (float 6004799503160661 1.0d0) -201)
+ (scale-float (float 3005918188789214 1.0d0) -254))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -96)
+ (scale-float (float 9007199254391466 1.0d0) -149)
+ (scale-float (float 6004799503160661 1.0d0) -202)
+ (scale-float (float 3002619653905886 1.0d0) -255))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -97)
+ (scale-float (float 9007199254653610 1.0d0) -150)
+ (scale-float (float 6004799503160661 1.0d0) -203)
+ (scale-float (float 3002413495475678 1.0d0) -256))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -98)
+ (scale-float (float 9007199254719146 1.0d0) -151)
+ (scale-float (float 6004799503160661 1.0d0) -204)
+ (scale-float (float 3002400610573790 1.0d0) -257))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -99)
+ (scale-float (float 9007199254735530 1.0d0) -152)
+ (scale-float (float 6004799503160661 1.0d0) -205)
+ (scale-float (float 3002399805267422 1.0d0) -258))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -100)
+ (scale-float (float 9007199254739626 1.0d0) -153)
+ (scale-float (float 6004799503160661 1.0d0) -206)
+ (scale-float (float 3002399754935774 1.0d0) -259))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -101)
+ (scale-float (float 9007199254740650 1.0d0) -154)
+ (scale-float (float 6004799503160661 1.0d0) -207)
+ (scale-float (float 3002399751790046 1.0d0) -260))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -102)
+ (scale-float (float 9007199254740906 1.0d0) -155)
+ (scale-float (float 6004799503160661 1.0d0) -208)
+ (scale-float (float 3002399751593438 1.0d0) -261))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -103)
+ (scale-float (float 9007199254740970 1.0d0) -156)
+ (scale-float (float 6004799503160661 1.0d0) -209)
+ (scale-float (float 3002399751581150 1.0d0) -262))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -104)
+ (scale-float (float 9007199254740986 1.0d0) -157)
+ (scale-float (float 6004799503160661 1.0d0) -210)
+ (scale-float (float 3002399751580382 1.0d0) -263))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -105)
+ (scale-float (float 9007199254740990 1.0d0) -158)
+ (scale-float (float 6004799503160661 1.0d0) -211)
+ (scale-float (float 3002399751580334 1.0d0) -264))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -106)
+ (scale-float (float 9007199254740991 1.0d0) -159)
+ (scale-float (float 6004799503160661 1.0d0) -212)
+ (scale-float (float 3002399751580331 1.0d0) -265))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -107)
+ (scale-float (float 9007199254740991 1.0d0) -160)
+ (scale-float (float 8256599316845909 1.0d0) -213)
+ (scale-float (float 3002399751580331 1.0d0) -266))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -108)
+ (scale-float (float 9007199254740991 1.0d0) -161)
+ (scale-float (float 8819549270267221 1.0d0) -214)
+ (scale-float (float 3002399751580331 1.0d0) -267))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -109)
+ (scale-float (float 9007199254740991 1.0d0) -162)
+ (scale-float (float 8960286758622549 1.0d0) -215)
+ (scale-float (float 3002399751580331 1.0d0) -268))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -110)
+ (scale-float (float 9007199254740991 1.0d0) -163)
+ (scale-float (float 8995471130711381 1.0d0) -216)
+ (scale-float (float 3002399751580331 1.0d0) -269))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -111)
+ (scale-float (float 9007199254740991 1.0d0) -164)
+ (scale-float (float 9004267223733589 1.0d0) -217)
+ (scale-float (float 3002399751580331 1.0d0) -270))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -112)
+ (scale-float (float 9007199254740991 1.0d0) -165)
+ (scale-float (float 9006466246989141 1.0d0) -218)
+ (scale-float (float 3002399751580331 1.0d0) -271))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -113)
+ (scale-float (float 9007199254740991 1.0d0) -166)
+ (scale-float (float 9007016002803029 1.0d0) -219)
+ (scale-float (float 3002399751580331 1.0d0) -272))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -114)
+ (scale-float (float 9007199254740991 1.0d0) -167)
+ (scale-float (float 9007153441756501 1.0d0) -220)
+ (scale-float (float 3002399751580331 1.0d0) -273))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -115)
+ (scale-float (float 9007199254740991 1.0d0) -168)
+ (scale-float (float 9007187801494869 1.0d0) -221)
+ (scale-float (float 3002399751580331 1.0d0) -274))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -116)
+ (scale-float (float 9007199254740991 1.0d0) -169)
+ (scale-float (float 9007196391429461 1.0d0) -222)
+ (scale-float (float 3002399751580331 1.0d0) -275))
+ (%make-qd-d (scale-float (float 9007199254740991 1.0d0) -117)
+ (scale-float (float 9007199254740991 1.0d0) -170)
+ (scale-float (float 9007198538913109 1.0d0) -223)
+ (scale-float (float 3002399751580331 1.0d0) -276))
[35 lines skipped]
1
0