Switch to using atan-qd/taylor as the default atan implementation.
authorRaymond Toy <toy.raymond@gmail.com>
Tue, 26 Nov 2013 17:55:46 +0000 (09:55 -0800)
committerRaymond Toy <toy.raymond@gmail.com>
Tue, 26 Nov 2013 17:55:59 +0000 (09:55 -0800)
qd-fun.lisp:
o Use atan-qd/taylor and atan2-qd/taylor as the atan-qd and atan2-qd
  implementation.
o Fix bug in atan-qd/taylor. We need a special case for x small where
  we don't want to use the table. (Found by our rt-test.  Yay!)

rt-tests.lisp:
o Add comment for oct.atan.4 for why it works.

qd-fun.lisp
rt-tests.lisp

index 9564dcc..207a51b 100644 (file)
@@ -1188,7 +1188,10 @@ is the cosine of A"
      ;; atan(x) = pi/2-atan(1/x), for x > 1
      (sub-qd +qd-pi/2+
             (atan-qd/taylor (div-qd +qd-one+ y))))
+    ((qd-<= y (aref octi::+qd-atan-partition+ 0))
+     (atan-taylor y))
     (t
+     ;; The general case where the table is needed.
      (let* ((i (find-atan-partition y))
            (atan-xi (mul-qd +qd-pi/2+
                             (rational-to-qd (/ (1- (+ i 2))
@@ -1244,6 +1247,8 @@ is the cosine of A"
          (neg-qd +qd-pi/4+))))
 
   (let ((arg (atan-qd/taylor (div-qd y (abs-qd x)))))
+    ;; arg = atan2(y,|x|), so |arg| <= pi/2.  Handle the case when x
+    ;; is negative.
     (cond ((minusp-qd x)
           (if (minusp-qd y)
               (- (neg-qd +qd-pi+) arg)
@@ -1254,12 +1259,12 @@ is the cosine of A"
 (defun atan2-qd (y x)
   "atan2(y, x) = atan(y/x), but carefully handling the quadrant"
   (declare (type %quad-double y x))
-  (atan2-qd/newton y x))
+  (atan2-qd/taylor y x))
 
 (defun atan-qd (y)
   "Return the arc tangent of the %QUAD-DOUBLE number Y"
   (declare (type %quad-double y))
-  (atan-qd/newton y))
+  (atan-qd/taylor y))
 
 (defun asin-qd (a)
   "Return the arc sine of the %QUAD-DOUBLE number A"
index 69f6ee7..289131e 100644 (file)
   nil)
 
 (rt:deftest oct.atan.4
+    ;; We know atan(10^100) is pi/2 because atan(10^100) =
+    ;; pi/2-atan(10^-100) and atan(10^-100) is approx 10^-100, which
+    ;; is too small to affect pi/2.
     (let* ((arg #q1q100)
           (y (/ (atan arg) +pi+))
           (true #q.5))