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

[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  0d5870201359817c679921a2d740fdd1697469b2 (commit)
      from  405df618a38d3b8ddaae691f865bbf068e931ac5 (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 0d5870201359817c679921a2d740fdd1697469b2
Author: Raymond Toy <toy.raymond at gmail.com>
Date:   Wed Mar 21 18:45:44 2012 -0700

    Implement psi and fix exp-integral-e for integral values of v.  Needs
    some more work.

diff --git a/qd-gamma.lisp b/qd-gamma.lisp
index 48f8345..02b5547 100644
--- a/qd-gamma.lisp
+++ b/qd-gamma.lisp
@@ -479,17 +479,44 @@
 		  (* (- k)
 		     (+ k v -1)))))))
 
+
+;; For v not an integer:
+;;
+;; E(v,z) = gamma(1-v)*z^(v-1) - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf)
+;;
+;; For v an integer:
+;;
+;; E(v,z) = (-z)^(v-1)/(v-1)!*(psi(v)-log(z))
+;;          - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf, k != n-1)
+;;
 (defun s-exp-integral-e (v z)
   ;; E(v,z) = gamma(1-v)*z^(v-1) - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf)
   (let ((-z (- z))
 	(-v (- v))
 	(eps (epsilon z)))
-    (loop for k from 0
-	  for term = 1 then (* term (/ -z k))
-	  for sum = (/ (- 1 v)) then (+ sum (/ term (+ k 1 -v)))
-	  when (< (abs term) (* (abs sum) eps))
-	  return (- (* (gamma (- 1 v)) (expt z (- v 1)))
-		    sum))))
+    (if (and (realp v)
+	     (= v (ftruncate v)))
+	;; v is an integer
+	(let ((n (truncate v)))
+	  (- (* (/ (expt -z (- v 1))
+		   (gamma v))
+		(- (psi v) (log z)))
+	     (loop for k from 0 below n
+		   for term = 1 then (* term (/ -z k))
+		   for sum = (/ (- 1 v)) then (+ sum (/ term (+ k 1 -v)))
+		   when (< (abs term) (* (abs sum) eps))
+		     return sum)
+	     (loop for k from n
+		   for term = 1 then (* term (/ -z k))
+		   for sum = 0 then (+ sum (/ term (+ k 1 -v)))
+		   when (< (abs term) (* (abs sum) eps))
+		     return sum)))
+	(loop for k from 0
+	      for term = 1 then (* term (/ -z k))
+	      for sum = (/ (- 1 v)) then (+ sum (/ term (+ k 1 -v)))
+	      when (< (abs term) (* (abs sum) eps))
+		return (- (* (gamma (- 1 v)) (expt z (- v 1)))
+				    sum)))))
 
 (defun exp-integral-e (v z)
   "Exponential integral E:
@@ -679,3 +706,109 @@
     (if (and (realp z) (plusp z))
 	(realpart (ci z))
 	(ci z))))
+
+(defconstant bern-values
+  (make-array 55
+	      :initial-contents
+	      '(1
+		-1/2
+		1/6
+		0
+		-1/30
+		0
+		1/42
+		0
+		-1/30
+		0
+		5/66
+		0
+		-691/2730
+		0
+		7/6
+		0
+		-3617/510
+		0
+		43867/798
+		0
+		-174611/330
+		0
+		854513/138
+		0
+		-236364091/2730
+		0
+		8553103/6
+		0
+		-23749461029/870
+		0
+		8615841276005/14322
+		0
+		-7709321041217/510
+		0
+		2577687858367/6
+		0
+		-26315271553053477373/1919190
+		0
+		2929993913841559/6
+		0
+		-261082718496449122051/13530
+		0
+		1520097643918070802691/1806
+		0
+		-27833269579301024235023/690
+		0
+		596451111593912163277961/282
+		0
+		-5609403368997817686249127547/46410
+		0
+		495057205241079648212477525/66
+		0
+		-801165718135489957347924991853/1590
+		0
+		29149963634884862421418123812691/798
+		)))
+		
+(defun bern (k)
+  (aref bern-values k))
+
+(defun psi (z)
+  "Digamma function defined by
+
+  - %gamma + sum(1/k-1/(k+z-1), k, 1, inf)
+
+  where %gamma is Euler's constant"
+
+  ;; A&S 6.3.7:  Reflection formula
+  ;;
+  ;;   psi(1-z) = psi(z) + %pi*cot(%pi*z)
+  ;;
+  ;; A&S 6.3.6:  Recurrence formula
+  ;;
+  ;;   psi(n+z) = 1/(z+n-1)+1/(z+n-2)+...+1/(z+2)+1/(1+z)+psi(1+z)
+  ;;
+  ;; A&S 6.3.8:  Asymptotic formula
+  ;;
+  ;;   psi(z) ~ log(z) - sum(bern(2*n)/(2*n*z^(2*n)), n, 1, inf)
+  ;;
+  ;; So use reflection formula if Re(z) < 0.  For z > 0, use the recurrence
+  ;; formula to increase the argument and then apply the asymptotic formula.
+
+  (cond ((minusp (realpart z))
+	 (- (psi (- 1 z))
+	    (* +pi+ (/ (tan (* +pi+ z))))))
+	(t
+	 (let* ((k (* 2 (1+ (floor (* .41 (- (log (epsilon z) 10)))))))
+		(m 0)
+		(y (expt (+ z k) 2))
+		(x 0))
+	   (loop for i from 1 upto (floor k 2) do
+	     (progn
+	       (incf m (+ (/ (+ z i i -1))
+			  (/ (+ z i i -2))))
+	       (setf x (/ (+ x (/ (bern (+ k 2 (* -2 i)))
+				  (- k i i -2)))
+			  y))))
+	   (- (log (+ z k))
+	      (/ (* 2 (+ z k)))
+	      x
+	      m)))))
+

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

Summary of changes:
 qd-gamma.lisp |  145 ++++++++++++++++++++++++++++++++++++++++++++++++++++++--
 1 files changed, 139 insertions(+), 6 deletions(-)


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