[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 64c424b7f41413424bbab93d4f5dd4f8461b784d
Raymond Toy
rtoy at common-lisp.net
Tue Apr 17 04:45:13 UTC 2012
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "OCT: A portable Lisp implementation for quad-double precision floats".
The branch, master has been updated
via 64c424b7f41413424bbab93d4f5dd4f8461b784d (commit)
via 8c5195a87137fd61952a175a5dae676c70b480ef (commit)
from 78801d6381aaaf4f21967582680f26889582db60 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 64c424b7f41413424bbab93d4f5dd4f8461b784d
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Mon Apr 16 20:45:12 2012 -0700
Add iterative versions for some functions.
diff --git a/qd-bessel.lisp b/qd-bessel.lisp
index 14743af..37a4c8e 100644
--- a/qd-bessel.lisp
+++ b/qd-bessel.lisp
@@ -78,6 +78,21 @@
(* (ash n -1) p)
(- (* (+ a (ash n -1)) p))))))))
+;; Use the recursion
+(defun bk-iter (k p old-bk)
+ (with-floating-point-contagion (p old-bk)
+ (if (zerop k)
+ (* (sqrt (/ (float-pi p) 8))
+ (let ((rp (sqrt p)))
+ (/ (erf rp)
+ rp)))
+ (- (* (- k 1/2)
+ (/ old-bk (* 2 p)))
+ (/ (exp (- p))
+ p
+ (ash 1 (+ k 1))
+ (sqrt (float 2 (realpart p))))))))
+
;; exp-arc I function, as given in the Laguerre paper
;;
;; I(p, q) = 4*exp(p) * sum(g[k](-2*%i*q)/(2*k)!*B[k](p), k, 0, inf)
@@ -178,6 +193,35 @@
(format t " term = ~S~%" term)
(format t " sum - ~S~%" sum)))))
+(defun exp-arc-i-3 (p q)
+ (let* ((v (* #c(0 -2) q))
+ (v2 (expt v 2))
+ (eps (epsilon (realpart p))))
+ (do* ((k 0 (1+ k))
+ (bk (bk 0 p)
+ (bk-iter k p bk))
+ ;; Compute g[k](p)/(2*k)!, not r[2*k+1](p)/(2*k)!
+ (ratio 1
+ (* ratio (/ (+ v2 (expt (1- (* 2 k)) 2))
+ (* 2 k (1- (* 2 k))))))
+ (term (* ratio bk)
+ (* ratio bk))
+ (sum term (+ sum term)))
+ ((< (abs term) (* (abs sum) eps))
+ (when *debug-exparc*
+ (format t "Final k= ~D~%" k)
+ (format t " bk = ~S~%" bk)
+ (format t " ratio = ~S~%" ratio)
+ (format t " term = ~S~%" term)
+ (format t " sum - ~S~%" sum))
+ (* sum 4 (exp p)))
+ (when *debug-exparc*
+ (format t "k = ~D~%" k)
+ (format t " bk = ~S~%" bk)
+ (format t " ratio = ~S~%" ratio)
+ (format t " term = ~S~%" term)
+ (format t " sum - ~S~%" sum)))))
+
;; Not really just for Bessel J for integer orders, but in that case,
;; this is all that's needed to compute Bessel J. For other values,
@@ -212,6 +256,7 @@
;;
;; alpha[n](z) = - exp(-z/2)/2^n/z + n/z*alpha[n-1](z)
;; beta[n]z) = ((-1)^n*exp(z/2)-exp(-z/2))/2^n/z + n/z*beta[n-1](z)
+;; = (-1)^n/(2^n)*2*sinh(z/2)/z + n/z*beta[n-1](z)
;;
;; We also note that
;;
@@ -223,12 +268,34 @@
(/ (incomplete-gamma (1+ n) (/ z 2))
(expt z (1+ n)))))
+(defun alpha-iter (n z alpha-old)
+ (if (zerop n)
+ ;; (1- exp(-z/2))/z.
+ (/ (- 1 (exp (* z -1/2)))
+ z)
+ (- (* (/ n z) alpha-old)
+ (/ (exp (- (* z 1/2)))
+ z
+ (ash 1 n)))))
+
(defun beta (n z)
(let ((n (float n (realpart z))))
(/ (- (incomplete-gamma (1+ n) (/ z 2))
(incomplete-gamma (1+ n) (/ z -2)))
(expt z (1+ n)))))
+(defun beta-iter (n z old-beta)
+ (if (zerop n)
+ ;; integrate(exp(-z*s),s,-1/2,1/2)
+ ;; = (exp(z/2)-exp(-z/2)/z
+ ;; = 2*sinh(z/2)/z
+ ;; = sinh(z/2)/(z/2)
+ (* 2 (/ (sinh (* 1/2 z)) z))
+ (+ (* n (/ old-beta z))
+ (* (/ (sinh (* 1/2 z)) (* 1/2 z))
+ (scale-float (float (if (evenp n) 1 -1) (realpart z)) (- n))))))
+
+
;; a[0](k,v) := (k+sqrt(k^2+1))^(-v);
;; a[1](k,v) := -v*a[0](k,v)/sqrt(k^2+1);
;; a[n](k,v) := 1/(k^2+1)/(n-1)/n*((v^2-(n-2)^2)*a[n-2](k,v)-k*(n-1)*(2*n-3)*a[n-1](k,v));
@@ -305,6 +372,26 @@
(format t " term = ~S~%" term)
(format t " sum = ~S~%" sum)))))
+(defun sum-ab-2 (big-n v z)
+ (let ((eps (epsilon (realpart z))))
+ (an-clrhash)
+ (do* ((n 0 (+ 1 n))
+ (alphan (alpha-iter 0 z 0)
+ (alpha-iter n z alphan))
+ (betan (beta-iter 0 z 0)
+ (beta-iter n z betan))
+ (term (+ (* alphan (an n 0 v))
+ (* betan (sum-an big-n n v z)))
+ (+ (* alphan (an n 0 v))
+ (* betan (sum-an big-n n v z))))
+ (sum term (+ sum term)))
+ ((<= (abs term) (* eps (abs sum)))
+ sum)
+ (when nil
+ (format t "n = ~D~%" n)
+ (format t " term = ~S~%" term)
+ (format t " sum = ~S~%" sum)))))
+
;; Convert to iteration instead of this quick-and-dirty memoization?
(let ((hash (make-hash-table :test 'equal)))
(defun %big-a-clrhash ()
commit 8c5195a87137fd61952a175a5dae676c70b480ef
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Mon Apr 16 19:30:12 2012 -0700
Make big-n a defvar so we can change it easily.
diff --git a/qd-bessel.lisp b/qd-bessel.lisp
index 48a0dac..14743af 100644
--- a/qd-bessel.lisp
+++ b/qd-bessel.lisp
@@ -410,6 +410,7 @@
;; bessel_i(v,z) = exp(-v*%pi*%i/2)*bessel_j(v, %i*z)
;; when Im(z) >> Re(z)
;;
+(defvar *big-n* 100)
(defun bessel-j (v z)
(let ((vv (ftruncate v)))
;; Clear the caches for now.
@@ -420,7 +421,7 @@
(integer-bessel-j-exp-arc v z))
(t
;; Need to fine-tune the value of big-n.
- (let ((big-n 10)
+ (let ((big-n *big-n*)
(vpi (* v (float-pi (realpart z)))))
(+ (integer-bessel-j-exp-arc v z)
(if (= vv v)
-----------------------------------------------------------------------
Summary of changes:
qd-bessel.lisp | 90 +++++++++++++++++++++++++++++++++++++++++++++++++++++++-
1 files changed, 89 insertions(+), 1 deletions(-)
hooks/post-receive
--
OCT: A portable Lisp implementation for quad-double precision floats
More information about the oct-cvs
mailing list