October 2007 - oct-cvs - mailman3.common-lisp.net

[oct-cvs] Oct commit: oct qd-class.lisp qd-fun.lisp qd-package.lisp qd.lisp
by rtoy 16 Oct '07

16 Oct '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv5153 Modified Files: qd-class.lisp qd-fun.lisp qd-package.lisp qd.lisp Log Message: qd-package.lisp: o Don't :USE CMUCL's EXTENSIONS package anymore. Import just the symbols we need. qd-class.lisp: o No need to use package qualifiers. qd-fun.lisp: o Need package qualifier for MAYBE-INLINE. qd.lisp: o Need package qualifier for *INLINE-EXPANSION-LIMIT*

1 0

[oct-cvs] Oct commit: oct qd-extra.lisp
by rtoy 15 Oct '07

15 Oct '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv23713 Modified Files: qd-extra.lisp Log Message: o Oops. Fix up a few IN-PACKAGE's for the new package names. --- /project/oct/cvsroot/oct/qd-extra.lisp 2007/10/11 17:47:08 1.3 +++ /project/oct/cvsroot/oct/qd-extra.lisp 2007/10/15 18:54:02 1.4 @@ -32,7 +32,7 @@ ;;; The end of the file also includes some tests of the different ;;; implementations for speed. -(in-package #:qdi) +(in-package #:octi) ;; This works but seems rather slow, so we don't even compile it. #+(or)

1 0

[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

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

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

[oct-cvs] Oct commit: oct qd.lisp rt-tests.lisp
by rtoy 15 Oct '07

15 Oct '07

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

[oct-cvs] Oct commit: oct qd-fun.lisp
by rtoy 15 Oct '07

15 Oct '07

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

[oct-cvs] Oct commit: oct qd-const.lisp
by rtoy 14 Oct '07

14 Oct '07

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

[oct-cvs] Oct commit: oct qd-fun.lisp
by rtoy 13 Oct '07

13 Oct '07

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

[oct-cvs] Oct commit: oct qd.lisp
by rtoy 13 Oct '07

13 Oct '07

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

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