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))