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

[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  b609574d330c7b0c9a6de588d3884a00f9aad13b (commit)
       via  5079d741cabce7dbaa84410fdfe565dfff03a51c (commit)
      from  5da3f9d5ccf7ddc9fe3f9334391984426acaeb9d (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 b609574d330c7b0c9a6de588d3884a00f9aad13b
Author: Raymond Toy <toy.raymond at gmail.com>
Date:   Sun Dec 4 21:29:17 2011 -0800

    Signal error for gamma of negative integers.

diff --git a/qd-gamma.lisp b/qd-gamma.lisp
index 8cab1bb..05f85d4 100644
--- a/qd-gamma.lisp
+++ b/qd-gamma.lisp
@@ -134,14 +134,18 @@
 
 (defun gamma-aux (z limit nterms)
   (let ((precision (float-contagion z)))
-    (cond ((minusp (realpart z))
+    (cond ((<= (realpart z) 0)
 	   ;; Use reflection formula if realpart(z) < 0:
 	   ;;  gamma(-z) = -pi*csc(pi*z)/gamma(z+1)
 	   ;; or
 	   ;;  gamma(z) = pi*csc(pi*z)/gamma(1-z)
-	   (/ (float-pi z)
-	      (sin (* (float-pi z) z))
-	      (gamma-aux (- 1 z) limit nterms)))
+	   (if (and (realp z)
+		    (= (truncate z) z))
+	       ;; Gamma of a negative integer is infinity.  Signal an error
+	       (error "Gamma of non-positive integer ~S" z)
+	       (/ (float-pi z)
+		  (sin (* (float-pi z) z))
+		  (gamma-aux (- 1 z) limit nterms))))
 	  ((and (zerop (imagpart z))
 		(= z (truncate z)))
 	   ;; We have gamma(n) where an integer value n and is small
@@ -623,4 +627,30 @@
     (if (and (realp z) (plusp z))
 	(realpart (ci z))
 	(ci z))))
-  
\ No newline at end of file
+
+;;; Exponential integral e defined by
+;;;
+;;; E(v,z) = z^(v-1) * integrate(t^(-v)*exp(-t), t, z, inf);
+;;;
+;;; for |arg(z)| < pi.
+;;;
+;;;
+;;; We use the continued fraction
+;;;
+;;; E(v,z) = exp(-z)/cf(z)
+;;;
+;; where the continued fraction cf(z) is
+;;
+;; a[k] = -k*(k+v-1)
+;; b[k] = v + 2*k + z
+;;
+;; for k = 1, inf
+
+(defun expintegral-e (v z)
+  (let ((z+v (+ z v)))
+    (/ (exp (- z))
+       (lentz #'(lambda (k)
+		  (+ z+v (* 2 k)))
+	      #'(lambda (k)
+		  (* (- k)
+		     (1- (+ k v))))))))

commit 5079d741cabce7dbaa84410fdfe565dfff03a51c
Author: Raymond Toy <toy.raymond at gmail.com>
Date:   Sun Dec 4 21:28:27 2011 -0800

    Fix bug in CEILING and FCEILING; add tests.

diff --git a/qd-methods.lisp b/qd-methods.lisp
index e4df131..bb851c4 100644
--- a/qd-methods.lisp
+++ b/qd-methods.lisp
@@ -652,7 +652,7 @@ underlying floating-point format"
   (multiple-value-bind (f rem)
       (floor x y)
     (if (zerop rem)
-	(values (+ f 1)
+	(values f
 		rem)
 	(values (+ f 1)
 		(- rem 1)))))
@@ -661,7 +661,7 @@ underlying floating-point format"
   (multiple-value-bind (f rem)
       (ffloor x y)
     (if (zerop rem)
-	(values (+ f 1)
+	(values f
 		rem)
 	(values (+ f 1)
 		(- rem 1)))))
diff --git a/rt-tests.lisp b/rt-tests.lisp
index acc8b74..736f7f8 100644
--- a/rt-tests.lisp
+++ b/rt-tests.lisp
@@ -55,6 +55,74 @@
 
 ;;; Some simple tests from the Yozo Hida's qd package.
 
+(rt:deftest ceiling-d.1
+    (multiple-value-list (ceiling -50d0))
+  (-50 0d0))
+
+(rt:deftest ceiling-d.2
+    (let ((z -50.1d0))
+      (multiple-value-bind (res rem)
+	  (ceiling -50.1d0)
+	(list res (= z (+ res rem)))))
+  (-50 t))
+
+(rt:deftest ceiling-q.1
+    (multiple-value-bind (res rem)
+	(ceiling #q-50q0)
+      (list res (zerop rem)))
+  (-50 t))
+
+(rt:deftest ceiling-q.2
+    (let ((z #q-50.1q0))
+      (multiple-value-bind (res rem)
+	  (ceiling z)
+	(list res (= z (+ res rem)))))
+  (-50 t))
+
+(rt:deftest truncate-d.1
+    (multiple-value-list (truncate -50d0))
+  (-50 0d0))
+
+(rt:deftest truncate-q.1
+    (multiple-value-bind (res rem)
+	(truncate #q-50q0)
+      (list res (zerop rem)))
+  (-50 t))
+
+(rt:deftest fceiling-d.1
+    (multiple-value-list (fceiling -50d0))
+  (-50d0 0d0))
+
+(rt:deftest fceiling-d.2
+    (let ((z -50.1d0))
+      (multiple-value-bind (res rem)
+	  (fceiling -50.1d0)
+	(list res (= z (+ res rem)))))
+  (-50d0 t))
+
+(rt:deftest fceiling-q.1
+    (multiple-value-bind (res rem)
+	(fceiling #q-50q0)
+      (list (= res -50) (zerop rem)))
+  (t t))
+
+(rt:deftest fceiling-q.2
+    (let ((z #q-50.1q0))
+      (multiple-value-bind (res rem)
+	  (fceiling z)
+	(list (= res -50) (= z (+ res rem)))))
+  (t t))
+
+(rt:deftest ftruncate-d.1
+    (multiple-value-list (ftruncate -50d0))
+  (-50d0 0d0))
+
+(rt:deftest ftruncate-q.1
+    (multiple-value-bind (res rem)
+	(ftruncate #q-50q0)
+      (list (= res -50) (zerop rem)))
+  (t t))
+
 ;; Pi via Machin's formula
 (rt:deftest oct.pi.machin
     (let* ((*standard-output* *null*)
@@ -1161,7 +1229,7 @@
        for y = (+ 1 (random 100d0))
        for g = (abs (gamma (complex 0 y)))
        for true = (sqrt (/ pi y (sinh (* pi y))))
-       for result = (check-accuracy 45 g true)
+       for result = (check-accuracy 44 g true)
        when result
        append (list (list (list k y) result)))
   nil)
@@ -1294,3 +1362,4 @@
 	   (true (fresnel-s-series z)))
       (check-accuracy 212 s true))
   nil)
+

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

Summary of changes:
 qd-gamma.lisp   |   40 +++++++++++++++++++++++++++----
 qd-methods.lisp |    4 +-
 rt-tests.lisp   |   71 ++++++++++++++++++++++++++++++++++++++++++++++++++++++-
 3 files changed, 107 insertions(+), 8 deletions(-)


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