# [oct-cvs] Oct commit: oct qd-extra.lisp

rtoy rtoy at common-lisp.net
Fri Oct 26 15:48:15 UTC 2007

```Update of /project/oct/cvsroot/oct
In directory clnet:/tmp/cvs-serv1229

Modified Files:
qd-extra.lisp
Log Message:

Experiments show that it is 25% faster than exp-qd/reduce.  However,
it loses 3 bits of accuracy compared with exp-qd/reduce.

--- /project/oct/cvsroot/oct/qd-extra.lisp	2007/10/15 18:54:02	1.4
+++ /project/oct/cvsroot/oct/qd-extra.lisp	2007/10/26 15:48:15	1.5
@@ -1018,6 +1018,71 @@
(mul-qd +cordic-scale+ y))))

+;; Evaluate the polynomial poly at the point q.  The polynomial is
+;; a list of the coefficients arranged in descending powers.
+(defun poly-eval-qd (q poly)
+  (let ((sum (%make-qd-d 0d0 0d0 0d0 0d0)))
+    (dolist (c poly)
+      (setf sum (add-qd (mul-qd q sum) c)))
+    sum))
+
+;; Like poly-eval-qd, except the polynomial is a list of double-floats
+(defun poly-eval-qd-d (q poly)
+  (let ((sum (%make-qd-d 0d0 0d0 0d0 0d0)))
+    (dolist (c poly)
+      (setf sum (add-qd-d (mul-qd q sum) c)))
+    sum))
+
+
+;; This idea is from Richard Fateman.  I've added a few additional
+;; notes, but these are mostly from Richard.
+;;
+;; We can evaluate exp(x) for |x|< 0.01, say, by a pade approximation,
+;; computed using Maxima:
+;;
+;; taylor(exp(x), x, 0, 20)\$
+;;
+;;
+;; [(x^10+110*x^9+5940*x^8+205920*x^7+5045040*x^6+90810720*x^5+1210809600*x^4
+;;  +11762150400*x^3+79394515200*x^2+335221286400*x^+670442572800)
+;;  /(x^10-110*x^9+
+;;  5940*x^8-205920*x^7+5045040*x^6-90810720*x^5+1210809600*x^4-11762150400*x^3+
+;;  79394515200*x^2-335221286400*x^+670442572800)]
+;;
+;; The numerator and denominator have the same coefficients with
+;; different signs on odd terms.  so we note that num and den here can
+;; be evaluated as f(x^2)+x*g(x^2) and f(x^2)-x*g(x^2) respectively using
+;; half the number of multiplies.
+
+  '(1.0d0 5940.0d0 5045040.0d0  1.2108096d+9
+    7.93945152d+10 6.704425728d+11))
+
+  '(110.0d0  205920.0d0  9.081072d+7
+    1.17621504d+10  3.352212864d+11 ))
+
+  ;; exp(x) = (f(x^2) + x*g(x^2))/(f(x^2) - x*g(x^2))
+  (let* ((x^2 (sqr-qd x))
+	 (x*gx (mul-qd x (poly-eval-qd-d x^2 g-exp-pade))))
+	    (sub-qd fx x*gx))))
+
+  ;; Same idea as in exp-qd/reduce, except we use our Pade
+  ;; approximation instead of a Taylor series.
+  (let* ((k 256)
+	 (z (truncate (qd-0 (nint-qd (div-qd a +qd-log2+)))))
+	 (r (div-qd-d (sub-qd a (mul-qd-d +qd-log2+ (float z 1d0)))
+		      (float k 1d0)))
+
+    (scale-float-qd (npow e k) z)))
+
+
;; Some timing and consing tests.
;;
;; The tests are run using the following:
@@ -1035,11 +1100,13 @@
;; exp-qd/reduce	2.06	 3.18	10.46	2.76		 6.12
;; expm1-qd/series	8.81	12.24	18.87	3.26		29.0
;; expm1-qd/dup		5.68	 4.34	18.47	3.64		18.78
;;
;; Consing (MB)		Sparc
;; exp-qd/reduce	 45   	 45   	 638   	44.4   		 45
;; expm1-qd/series	519   	519   	1201  	14.8   		519
;; expm1-qd/dup		 32   	 32   	1224   	32.0   		 32
;;
;; Speeds seem to vary quite a bit between architectures.
;;
@@ -1066,17 +1133,22 @@
(fixnum n))
(let ((y +qd-zero+))
-    #+cmu (gc :full t)
+    #+cmu (ext:gc :full t)
(format t "exp-qd/reduce~%")
(time (dotimes (k n)
(declare (fixnum k))
(setf y (exp-qd/reduce x))))
-    #+cmu (gc :full t)
+    #+cmu (ext:gc :full t)
+    (time (dotimes (k n)
+	    (declare (fixnum k))
+    #+cmu (ext:gc :full t)
(format t "expm1-qd/series~%")
(time (dotimes (k n)
(declare (fixnum k))
(setf y (expm1-qd/series x))))
-    #+cmu (gc :full t)
+    #+cmu (ext:gc :full t)
(format t "expm1-qd/duplication~%")
(time (dotimes (k n)
(declare (fixnum k))
@@ -1084,6 +1156,87 @@

))

+#||
+
+;; This bit of code is meant to be run with Clisp, which has arbitrary
+;; precision long floats.
+;;
+;; We use this to compare the accuracy of exp-qd/reduce and
+;; exp-qd/pade against Clisp's exp implementation.
+;;
+;; Some results:
+;;
+;; k = 159000
+;; red min    =       3.838L-69
+;; red max    =       4.318L-62
+;; red sumsq  =      6.076L-120
+;;
+;; So, worst case error is:
+;;
+;; exp-qd/reduce: 4.318L-62 (203.8 bits)
+;;
+;; We lose 3 bits for Pade.  exp-qd/reduce is "missing" 8 bits of
+;; accuracy.  Bummer.
+
+(in-package #:oct)
+
+;; 240 bit fractions for long-floats.  I hope this is accurate enough
+;; compared against our 212 bit fractions for quad-doubles.
+(setf (ext:long-float-digits) 240)
+(defun compare-exp (n &optional verbose)
+	(red-max -1L0)
+	(red-min 1L100)
+	(red-sum 0L0))
+    (dotimes (k n)
+      (let* ((f (random (make-qd 1/256)))
+	     (red-exp (octi::exp-qd/reduce (qd-value f))))
+	(flet ((qd-to-lf (qd)
+		 (octi:with-qd-parts (q0 q1 q2 q3)
+		     qd
+		   (+ (float q0 1L0)
+		      (float q1 1L0)
+		      (float q2 1L0)
+		      (float q3 1L0)))))
+	  (let* ((lf (qd-to-lf (qd-value f)))
+		 (l-exp (exp lf))
+		 (red-exp-lf (qd-to-lf red-exp))
+		 (red-diff (abs (- l-exp red-exp-lf))))
+	    (when verbose
+	      (format t "f     = ~S~%" f)
+	      (format t "l-exp =   ~A~%" l-exp)
+	      (format t "diff  = ~A~%" pade-diff))
+	    (when (and (plusp k) (zerop (mod k 1000)))
+	      (format t "k = ~A~%" k)
+	      (format t "red min    = ~15,3g~%" red-min)
+	      (format t "red max    = ~15,3g~%" red-max)
+	      (format t "red sumsq  = ~15,3g~%" red-sum))
+	    (setf red-max (max red-max red-diff))
+	    (setf red-min (min red-min red-diff))
+	    (incf red-sum (* red-diff red-diff)))
+	  )))
+	  red-max red-min red-sum
+	  )))
+
+||#
+
;; (time-log #c(3w0 0) 50000)
;;
;; Time (s)		Sparc	PPC	x86	PPC (fma)	Sparc2
@@ -1147,33 +1300,33 @@
(fixnum n))
(let ((y +qd-zero+))
-    #+cmu (gc :full t)
+    #+cmu (ext:gc :full t)
(format t "log-qd/newton~%")
(time (dotimes (k n)
(declare (fixnum k))
(setf y (log-qd/newton x))))
-    #+cmu (gc :full t)
+    #+cmu (ext:gc :full t)
(format t "log1p-qd/duplication~%")
(time (dotimes (k n)
(declare (fixnum k))
(setf y (log1p-qd/duplication x))))
-    #+cmu (gc :full t)
+    #+cmu (ext:gc :full t)
(format t "log-qd/agm~%")
(time (dotimes (k n)
(declare (fixnum k))
(setf y (log-qd/agm x))))

-    #+cmu (gc :full t)
+    #+cmu (ext:gc :full t)
(format t "log-qd/agm2~%")
(time (dotimes (k n)
(declare (fixnum k))
(setf y (log-qd/agm2 x))))
-    #+cmu (gc :full t)
+    #+cmu (ext:gc :full t)
(format t "log-qd/agm3~%")
(time (dotimes (k n)
(declare (fixnum k))
(setf y (log-qd/agm3 x))))
-    #+cmu (gc :full t)
+    #+cmu (ext:gc :full t)
(format t "log-qd/halley~%")
(time (dotimes (k n)
(declare (fixnum k))
@@ -1217,17 +1370,17 @@
(fixnum n))
(let ((y +qd-zero+)
(one +qd-one+))
-    #+cmu (gc :full t)
+    #+cmu (ext:gc :full t)
(format t "atan2-qd/newton~%")
(time (dotimes (k n)
(declare (fixnum k))
(setf y (atan2-qd/newton x one))))
-    #+cmu (gc :full t)
+    #+cmu (ext:gc :full t)
(format t "atan2-qd/cordic~%")
(time (dotimes (k n)
(declare (fixnum k))
(setf y (atan2-qd/cordic x one))))
-    #+cmu (gc :full t)
+    #+cmu (ext:gc :full t)
(format t "atan-qd/duplication~%")
(time (dotimes (k n)
(declare (fixnum k))
@@ -1260,12 +1413,12 @@
(fixnum n))
(let ((y +qd-zero+))
-    #+cmu (gc :full t)
+    #+cmu (ext:gc :full t)
(format t "tan-qd/cordic~%")
(time (dotimes (k n)
(declare (fixnum k))
(setf y (tan-qd/cordic x))))
-    #+cmu (gc :full t)
+    #+cmu (ext:gc :full t)
(format t "tan-qd/sincos~%")
(time (dotimes (k n)
(declare (fixnum k))

```