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

[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  e9cd1a46fcf2a5c0101b3473648cb242b55987e8 (commit)
      from  f203389e5e78d1f001f68447ac2a9dd86dcfbbf6 (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 e9cd1a46fcf2a5c0101b3473648cb242b55987e8
Author: Raymond Toy <rtoy at google.com>
Date:   Tue Apr 10 09:51:55 2012 -0700

    Fix bug in s-bessel-j.  Microoptimize integer-bessel-j-exp-arc for the
    case where v is an integer.

diff --git a/qd-bessel.lisp b/qd-bessel.lisp
index 9759044..367b847 100644
--- a/qd-bessel.lisp
+++ b/qd-bessel.lisp
@@ -171,18 +171,25 @@
 	(format t " sum   - ~S~%" sum)))))
 
 
-;; This currently only works for v an integer.
 ;;
 (defun integer-bessel-j-exp-arc (v z)
   (let* ((iz (* #c(0 1) z))
-	 (i+ (exp-arc-i-2 iz v))
-	 (i- (exp-arc-i-2 (- iz ) v)))
-    (/ (+ (* (cis (* v (float-pi i+) -1/2))
-	     i+)
-	  (* (cis (* v (float-pi i+) 1/2))
-	     i-))
-       (float-pi i+)
-       2)))
+	 (i+ (exp-arc-i-2 iz v)))
+    (cond ((= v (ftruncate v))
+	   ;; We can simplify the result
+	   (let ((c (cis (* v (float-pi i+) -1/2))))
+	     (/ (+ (* c i+)
+		   (* (conjugate c) (conjugate i+)))
+		(float-pi i+)
+		2)))
+	  (t
+	   (let ((i- (exp-arc-i-2 (- iz ) v)))
+	     (/ (+ (* (cis (* v (float-pi i+) -1/2))
+		      i+)
+		   (* (cis (* v (float-pi i+) 1/2))
+		      i-))
+		(float-pi i+)
+		2))))))
 
 ;; alpha[n](z) = integrate(exp(-z*s)*s^n, s, 0, 1/2)
 ;; beta[n](z)  = integrate(exp(-z*s)*s^n, s, -1/2, 1/2)
@@ -349,15 +356,18 @@
 	  (eps (epsilon z)))
       (do* ((k 0 (+ 1 k))
 	    (f (gamma (+ v 1))
-	       (* f (* k (+ v k))))
+	       (* k (+ v k)))
 	    (term (/ f)
 		  (/ (* (- term) z2/4) f))
 	    (sum term (+ sum term)))
 	   ((<= (abs term) (* eps (abs sum)))
 	    (* sum (expt (* z 1/2) v)))
-	(format t "k = ~D~%" k)
-	(format t " term = ~S~%" term)
-	(format t " sum  = ~S~%" sum)))))
+	#+nil
+	(progn
+	  (format t "k = ~D~%" k)
+	  (format t " f    = ~S~%" f)
+	  (format t " term = ~S~%" term)
+	  (format t " sum  = ~S~%" sum))))))
   
 (defun bessel-j (v z)
   (let ((vv (ftruncate v)))

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

Summary of changes:
 qd-bessel.lisp |   36 +++++++++++++++++++++++-------------
 1 files changed, 23 insertions(+), 13 deletions(-)


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