[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 6ab5226ca3e32b443d87934ec138ff0efc8aaecc

[Raymond Toy]

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  6ab5226ca3e32b443d87934ec138ff0efc8aaecc (commit)
       via  d795ba718dc53f591c82994811f50250aceec1d7 (commit)
       via  c0e12ddf6b61f571555d46c6f168e6bebabc80b1 (commit)
      from  7c5a3186070096ee93e16f2ddf51b2c84e7c5895 (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 6ab5226ca3e32b443d87934ec138ff0efc8aaecc
Author: Raymond Toy <toy.raymond at gmail.com>
Date:   Wed Apr 11 20:32:23 2012 -0700

    First cut at Bessel Y.  Not working yet.

diff --git a/qd-bessel.lisp b/qd-bessel.lisp
index e7ca034..588b5af 100644
--- a/qd-bessel.lisp
+++ b/qd-bessel.lisp
@@ -376,6 +376,7 @@
 	  (format t " term = ~S~%" term)
 	  (format t " sum  = ~S~%" sum))))))
 
+;; 
 ;; TODO:
 ;;  o For |z| <= 1 use the series.
 ;;  o Currently accuracy is not good for large z and half-integer
@@ -402,6 +403,41 @@
 		   (+ (/ -1 z)
 		      (sum-ab big-n v z)
 		      (sum-big-ia big-n v z)))))))))
+
+;; Bessel Y
+;;
+;; bessel_y(v, z) = 1/(2*%pi*%i)*(exp(-%i*v*%pi/2)*I(%i*v,z) - exp(%i*v*%pi/2)*I(-%i*z, v))
+;;                   + z/v/%pi*((1-cos(v*%pi)/z) + S(N,z,v)*cos(v*%pi)-S(N,z,-v))
+;;
+;; where
+;;
+;;   S(N,z,v) = sum(alpha[n](z)*a[n](0,v) + beta[n](z)*sum(exp(-k*z)*a[n](k,v),k,1,N),n,0,inf)
+;;               + sum(A[n](v)*I[n](N+1/2,z,v),n,0,inf)
+;;
+(defun bessel-y (v z)
+  (flet ((ipart (v z)
+	   (let* ((iz (* #c(0 1) z))
+		  (c+ (exp (* v (float-pi z) 1/2)))
+		  (c- (exp (* v (float-pi z) -1/2)))
+		  (i+ (exp-arc-i-2 iz v))
+		  (i- (exp-arc-i-2 (- iz) v)))
+	     (/ (- (* c- i+) (* c+ i-))
+		(* #c(0 2) (float-pi z)))))
+	 (s (big-n z v)
+	   (+ (sum-ab big-n v z)
+	      (sum-big-ia big-n v z))))
+    (let* ((big-n 100)
+	   (vpi (* v (float-pi z)))
+	   (c (cos vpi)))
+      (+ (ipart v z)
+	 (* (/ z vpi)
+	    (+ (/ (- 1 c)
+		  z)
+	       (* c
+		  (s big-n z v))
+	       (- (s big-n z (- v)))))))))
+	   
+  
 
 (defun paris-series (v z n)
   (labels ((pochhammer (a k)

commit d795ba718dc53f591c82994811f50250aceec1d7
Author: Raymond Toy <toy.raymond at gmail.com>
Date:   Wed Apr 11 20:15:47 2012 -0700

    Add TODO list for bessel-j.

diff --git a/qd-bessel.lisp b/qd-bessel.lisp
index 3958ff2..e7ca034 100644
--- a/qd-bessel.lisp
+++ b/qd-bessel.lisp
@@ -375,13 +375,25 @@
 	  (format t " f    = ~S~%" f)
 	  (format t " term = ~S~%" term)
 	  (format t " sum  = ~S~%" sum))))))
-  
+
+;; TODO:
+;;  o For |z| <= 1 use the series.
+;;  o Currently accuracy is not good for large z and half-integer
+;;    order.
+;;  o For real v and z, return a real number instead of complex.
+;;  o Handle the case of Re(z) < 0. (The formulas are for Re(z) > 0:
+;;    bessel_j(v,z*exp(m*%pi*%i)) = exp(m*v*%pi*%i)*bessel_j(v, z)
+;;  o The paper suggests using
+;;      bessel_i(v,z) = exp(-v*%pi*%i/2)*bessel_j(v, %i*z)
+;;    when Im(z) >> Re(z)
+;; 
 (defun bessel-j (v z)
   (let ((vv (ftruncate v)))
     (cond ((= vv v)
 	   ;; v is an integer
 	   (integer-bessel-j-exp-arc v z))
 	  (t
+	   ;; Need to fine-tune the value of big-n.
 	   (let ((big-n 100)
 		 (vpi (* v (float-pi (realpart z)))))
 	     (+ (integer-bessel-j-exp-arc v z)

commit c0e12ddf6b61f571555d46c6f168e6bebabc80b1
Author: Raymond Toy <toy.raymond at gmail.com>
Date:   Wed Apr 11 20:06:38 2012 -0700

    Correct the computation of alpha and beta.

diff --git a/qd-bessel.lisp b/qd-bessel.lisp
index 0197e7f..3958ff2 100644
--- a/qd-bessel.lisp
+++ b/qd-bessel.lisp
@@ -213,13 +213,13 @@
 
 (defun alpha (n z)
   (let ((n (float n (realpart z))))
-    (/ (cf-incomplete-gamma (1+ n) (/ z 2))
+    (/ (incomplete-gamma (1+ n) (/ z 2))
        (expt z (1+ n)))))
 
 (defun beta (n z)
   (let ((n (float n (realpart z))))
-    (/ (- (cf-incomplete-gamma (1+ n) (/ z 2))
-	  (cf-incomplete-gamma (1+ n) (/ z -2)))
+    (/ (- (incomplete-gamma (1+ n) (/ z 2))
+	  (incomplete-gamma (1+ n) (/ z -2)))
        (expt z (1+ n)))))
 
 ;; a[0](k,v) := (k+sqrt(k^2+1))^(-v);

-----------------------------------------------------------------------

Summary of changes:
 qd-bessel.lisp |   56 ++++++++++++++++++++++++++++++++++++++++++++++++++++----
 1 files changed, 52 insertions(+), 4 deletions(-)


hooks/post-receive
-- 
OCT:  A portable Lisp implementation for quad-double precision floats