Commit 4d09f5a3 authored by Raymond Toy's avatar Raymond Toy

Fix atan/taylor bug and add atan2/taylor.

qd-fun.lisp:
o Fix issue with atan-qd/taylor for large args. We use the
  relationship atan(x) = pi/2-atan(1/x) for x > 1.
o Add atan2-qd/taylor.  This appears to be much faster than
  atan2-qd/newton. 

timing.lisp:
o Add timing test for atan2-qd.
parent b9a21f1f
......@@ -1182,21 +1182,76 @@ is the cosine of A"
;; atan(x) = (2*i-1)*pi/(4*s) + atan(e)
;;
;; Note that qd-atan-partition array starts with X[2].
(if (minusp-qd y)
(neg-qd (atan-qd/taylor (neg-qd y)))
(let* ((i (find-atan-partition y))
(atan-xi (mul-qd +qd-pi/2+
(rational-to-qd (/ (1- (+ i 2))
octi::+qd-atan-partition-size+))))
(1/xi (div-qd +qd-one+ (aref octi::+qd-atan-nodes+ i)))
(e1 (div-qd (add-qd-d (sqr-qd 1/xi) 1d0)
(add-qd 1/xi y)))
(e (sub-qd 1/xi
e1)))
;;(format t "partition ~D~%" i)
;;(format t "atan(xi) = ~S~%" atan-xi)
;;(format t "e = ~S~%" e)
(add-qd atan-xi (atan-taylor e)))))
(cond
((minusp-qd y)
;; atan(-x) = -atan(x)
(neg-qd (atan-qd/taylor (neg-qd y))))
((qd-> y +qd-one+)
;; atan(x) = pi/2-atan(1/x), for x > 1
(sub-qd +qd-pi/2+
(atan-qd/taylor (div-qd +qd-one+ y))))
(t
(let* ((i (find-atan-partition y))
(atan-xi (mul-qd +qd-pi/2+
(rational-to-qd (/ (1- (+ i 2))
octi::+qd-atan-partition-size+))))
(1/xi (div-qd +qd-one+ (aref octi::+qd-atan-nodes+ i)))
(e1 (div-qd (add-qd-d (sqr-qd 1/xi) 1d0)
(add-qd 1/xi y)))
(e (sub-qd 1/xi
e1)))
;;(format t "partition ~D~%" i)
;;(format t "atan(xi) = ~S~%" atan-xi)
;;(format t "e = ~S~%" e)
(add-qd atan-xi (atan-taylor e))))))
(defun atan2-qd/taylor (y x)
(declare (type %quad-double y x)
#+nil
(optimize (speed 3) (space 0)))
(cond ((zerop-qd x)
;; x = 0
(cond ((zerop-qd y)
;; Both x and y are zero. Use the signs of x and y to
;; determine the result
(error "atan2(0,0)"))
(t
;; x = 0, but y is not. Use the sign of y.
(return-from atan2-qd/taylor
(cond ((plusp (float-sign (qd-0 y)))
+qd-pi/2+)
(t
(neg-qd +qd-pi/2+)))))))
((zerop-qd y)
;; y = 0.
(return-from atan2-qd/taylor
;; Use the sign of x and y to figure out the result.
(cond ((plusp (float-sign (qd-0 x)))
+qd-zero+)
((plusp (float-sign (qd-0 y)))
+qd-pi+)
(t
(neg-qd +qd-pi+))))))
(when (qd-= x y)
(return-from atan2-qd/taylor
(if (plusp-qd y)
+qd-pi/4+
+qd-3pi/4+)))
(when (qd-= x (neg-qd y))
(return-from atan2-qd/taylor
(if (plusp-qd y)
+qd-3pi/4+
(neg-qd +qd-pi/4+))))
(let ((arg (atan-qd/taylor (div-qd y (abs-qd x)))))
(cond ((minusp-qd x)
(if (minusp-qd y)
(- (neg-qd +qd-pi+) arg)
(- +qd-pi+ arg)))
(t
arg))))
(defun atan2-qd (y x)
"atan2(y, x) = atan(y/x), but carefully handling the quadrant"
......
......@@ -231,3 +231,57 @@
;; We see that the taylor series is 21 times faster (!) and conses 11
;; times less. That's a pretty nice gain, at the expense of two
;; 1024-element tables.
(defun time-atan2 (&optional (n 10000))
(declare (fixnum n))
(flet ((time-atan/newton ()
(let ((sum (octi::make-qd-d -1000d0)))
(declare (type octi::%quad-double sum)
(optimize (speed 3)))
(dotimes (k n)
(declare (fixnum k))
(setf sum (octi::atan2-qd/newton sum (mul-qd-d sum .1d0))))
sum))
(time-atan/taylor ()
(let ((sum (octi::make-qd-d -1000d0)))
(declare (type octi::%quad-double sum)
(optimize (speed 3)))
(dotimes (k n)
(declare (fixnum k))
(setf sum (octi::atan2-qd/taylor sum (mul-qd-d sum .1d0))))
sum)))
(format t "atan2-qd/newton ~d times~%" n)
#+cmu (ext:gc :full t)
(format t "sum = ~A~%" (time (time-atan/newton)))
(format t "atan2-qd/taylor ~d times~%" n)
#+cmu (ext:gc :full t)
(format t "sum = ~A~%" (time (time-atan/taylor)))))
;; Some timing results on an iMac, 3.06 GHz Core i3:
; atan2-qd/newton 10000 times
; Evaluation took:
; 6.37 seconds of real time
; 6.321544 seconds of user run time
; 0.026819 seconds of system run time
; 19,474,336,804 CPU cycles
; [Run times include 0.27 seconds GC run time]
; 0 page faults and
; 934,872,144 bytes consed.
;
; sum = #C(-1.670464979286058652105921398771024w0 1.39712341620073544419882362664172w-33)
;
; atan2-qd/taylor 10000 times
; Evaluation took:
; 0.32 seconds of real time
; 0.320343 seconds of user run time
; 0.001606 seconds of system run time
; 989,558,772 CPU cycles
; [Run times include 0.02 seconds GC run time]
; 0 page faults and
; 88,919,384 bytes consed.
;
; sum = #C(-1.670464979286058652105921398771027w0 4.47861132722031280908838833477765w-33)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment