Commit bc4055f3 authored by Raymond Toy's avatar Raymond Toy
Browse files

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.
parent 89bea403
......@@ -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))))
......
......@@ -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)
......
......@@ -646,7 +646,52 @@ that we can always return an integer"
(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 @@ that we can always return an integer"
;; 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 @@ that we can always return an integer"
;; 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 @@ that we can always return an integer"
;; 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 @@ that we can always return an integer"
+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 @@ that we can always return an integer"
(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))))
......
......@@ -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))))))))
......
......@@ -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 @@ underlying floating-point format"
(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 @@ underlying floating-point format"
(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))
......
......@@ -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+
......
......@@ -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)))
......
......@@ -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