[Maxima-commits] [git] Maxima CAS branch, master, updated. branch-5_31-base-205-g0fb4df3

[Raymond T. <rt...@us...>]

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 "Maxima CAS".

The branch, master has been updated
       via  0fb4df38c199843705ae24a63514be25fe1407ab (commit)
      from  29fe517aad136c878d4e6b33e661b5dbf6460fc9 (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 0fb4df38c199843705ae24a63514be25fe1407ab
Author: Raymond Toy <rt...@go...>
Date:   Sat Dec 7 15:06:00 2013 -0800

    Fix Bug 2668
    
    src/csimp2.lisp:
    o When the arg to gamma is small use the relationship gamma(z) =
      gamma(z+1)/z.
    
    tests/rtest_gamma.mac:
    o Add some tests for small z.

diff --git a/src/csimp2.lisp b/src/csimp2.lisp
index 0e5b204..919c6ca 100644
--- a/src/csimp2.lisp
+++ b/src/csimp2.lisp
@@ -378,21 +378,40 @@
            ;; Adding 4 digits in the call to bffac. For $fpprec up to about 256
            ;; and an argument up to about 500.0 the accuracy of the result is
            ;; better than 10^(-$fpprec).
-	   (let ((result (mfuncall '$bffac (m+ ($bfloat j) -1) (+ $fpprec 4))))
-	     ;; bigfloatp will round the result to the correct fpprec
-	     (bigfloatp result)))
+	   (let ((z (bigfloat:to ($bfloat j))))
+	     (cond
+	       ((bigfloat:<= (bigfloat:abs z) (bigfloat:sqrt (bigfloat:epsilon z)))
+		;; For small z, use gamma(z) = gamma(z+1)/z = z!/z
+		(div (mfuncall '$bffac
+			       ($bfloat j)
+			       (+ $fpprec 4))
+		     ($bfloat j)))
+	       (t
+		(let ((result (mfuncall '$bffac (m+ ($bfloat j) -1) (+ $fpprec 4))))
+		  ;; bigfloatp will round the result to the correct fpprec
+		  (bigfloatp result))))))
 	  ((complex-float-numerical-eval-p j)
            (complexify (gamma-lanczos (complex ($float ($realpart j))
                                                ($float ($imagpart j))))))
           ((complex-bigfloat-numerical-eval-p j)
-           ;; Adding 4 digits in the call to cbffac. See comment above.
-	   (let ((result
-		  (mfuncall '$cbffac 
-			    (add -1 ($bfloat ($realpart j)) 
-				 (mul '$%i ($bfloat ($imagpart j))))
-			    (+ $fpprec 4))))
-	     (add (bigfloatp ($realpart result))
-		  (mul '$%i (bigfloatp ($imagpart result))))))
+	   (let ((z (bigfloat:to ($bfloat j))))
+	     (cond
+	       ((bigfloat:<= (bigfloat:abs z)
+			     (bigfloat:sqrt (bigfloat:epsilon z)))
+		;; For small z, use gamma(z) = gamma(z+1)/z = z!/z
+		(to (bigfloat:/ (bigfloat:to (mfuncall '$cbffac 
+						       (to z)
+						       (+ $fpprec 4)))
+				z)))
+	       (t
+		;; Adding 4 digits in the call to cbffac. See comment above.
+		(let ((result
+		       (mfuncall '$cbffac 
+				 (add -1 ($bfloat ($realpart j)) 
+				      (mul '$%i ($bfloat ($imagpart j))))
+				 (+ $fpprec 4))))
+		  (add (bigfloatp ($realpart result))
+		       (mul '$%i (bigfloatp ($imagpart result)))))))))
           ((taylorize (mop x) (cadr x)))
           ((eq j '$inf) '$inf) ; Simplify to $inf to be more consistent.
           ((and $gamma_expand
@@ -510,51 +529,63 @@
 			 -.26190838401581408670e-4
 			 .36899182659531622704e-5))))
     (declare (type (simple-array flonum (15)) c))
-    (if (minusp (realpart z))
-	;; Use the reflection formula
-	;; -z*Gamma(z)*Gamma(-z) = pi/sin(pi*z)
-	;; or
-	;; Gamma(z) = pi/z/sin(pi*z)/Gamma(-z)
-	;;
-	;; If z is a negative integer, Gamma(z) is infinity.  Should
-	;; we test for this?  Throw an error?
-        ;; The test must be done by the calling routine.
-	(/ (float pi)
-	   (* (- z) (sin (* (float pi) z))
-	      (gamma-lanczos (- z))))
-	(let* ((z (- z 1))
-	       (zh (+ z 1/2))
-	       (zgh (+ zh 607/128))
-	       (ss 
-		(do ((sum 0.0)
-		     (pp (1- (length c)) (1- pp)))
-		    ((< pp 1)
-		     sum)
-		  (incf sum (/ (aref c pp) (+ z pp))))))
-          (let ((result 
-                 ;; We check for an overflow. The last positive value in 
-                 ;; double-float precicsion for which Maxima can calculate 
-                 ;; gamma is ~171.6243 (CLISP 2.46 and GCL 2.6.8)
-                 (ignore-errors
-		   (let ((zp (expt zgh (/ zh 2))))
-		     (* (sqrt (float (* 2 pi)))
-			(+ ss (aref c 0))
-			(* (/ zp (exp zgh)) zp))))))
-            (cond ((null result)
-                   ;; No result. Overflow.
-                   (merror (intl:gettext "gamma: overflow in GAMMA-LANCZOS.")))
-                  ((or (float-nan-p (realpart result))
-                       (float-inf-p (realpart result)))
-                   ;; Result, but beyond extreme values. Overflow.
-                   (merror (intl:gettext "gamma: overflow in GAMMA-LANCZOS.")))
-                  (t result)))))))
+    (cond
+      ((minusp (realpart z))
+       ;; Use the reflection formula
+       ;; -z*Gamma(z)*Gamma(-z) = pi/sin(pi*z)
+       ;; or
+       ;; Gamma(z) = pi/z/sin(pi*z)/Gamma(-z)
+       ;;
+       ;; If z is a negative integer, Gamma(z) is infinity.  Should
+       ;; we test for this?  Throw an error?
+       ;; The test must be done by the calling routine.
+       (/ (float pi)
+	  (* (- z) (sin (* (float pi) z))
+	     (gamma-lanczos (- z)))))
+      ((<= (abs z) (sqrt flonum-epsilon))
+       ;; For |z| small, use Gamma(z) = Gamma(z+1)/z
+       (/ (gamma-lanczos (+ 1 z))
+	  z))
+      (t
+       (let* ((z (- z 1))
+	      (zh (+ z 1/2))
+	      (zgh (+ zh 607/128))
+	      (ss 
+	       (do ((sum 0.0)
+		    (pp (1- (length c)) (1- pp)))
+		   ((< pp 1)
+		    sum)
+		 (incf sum (/ (aref c pp) (+ z pp))))))
+	 (let ((result 
+		;; We check for an overflow. The last positive value in 
+		;; double-float precicsion for which Maxima can calculate 
+		;; gamma is ~171.6243 (CLISP 2.46 and GCL 2.6.8)
+		(ignore-errors
+		  (let ((zp (expt zgh (/ zh 2))))
+		    (* (sqrt (float (* 2 pi)))
+		       (+ ss (aref c 0))
+		       (* (/ zp (exp zgh)) zp))))))
+	   (cond ((null result)
+		  ;; No result. Overflow.
+		  (merror (intl:gettext "gamma: overflow in GAMMA-LANCZOS.")))
+		 ((or (float-nan-p (realpart result))
+		      (float-inf-p (realpart result)))
+		  ;; Result, but beyond extreme values. Overflow.
+		  (merror (intl:gettext "gamma: overflow in GAMMA-LANCZOS.")))
+		 (t result))))))))
 
 (defun gammafloat (a)
   (let ((a (float a)))
     (cond ((minusp a)
+	   ;; Reflection formula to make it positive: gamma(x) =
+	   ;; %pi/sin(%pi*x)/x/gamma(-x)
 	   (/ (float (- pi))
 	      (* a (sin (* (float pi) a)))
 	      (gammafloat (- a))))
+	  ((<= a (sqrt flonum-epsilon))
+	   ;; Use gamma(x) = gamma(1+x)/x when x is very small
+	   (/ (gammafloat (+ 1 a))
+	      a))
 	  ((< a 10)
 	   (slatec::dgamma a))
 	  (t
diff --git a/tests/rtest_gamma.mac b/tests/rtest_gamma.mac
index e2c88f6..f8c5572 100755
--- a/tests/rtest_gamma.mac
+++ b/tests/rtest_gamma.mac
@@ -3695,6 +3695,29 @@ relerror(inverse_erf(- .9763545580611441 ),
  1d-15);
 true;
 
+/* Bug #2668 Bigfloat Gamma inaccurae for small inputs
+ *
+ * For these small values, gamma(z) = gamma(z+1)/z = 1/z
+ * since 1+z = 1 and gamma(1) = 1.
+ */
+relerror(gamma(2.0^-256), 2.0^256, 1d-14);
+true;
+
+relerror(gamma(2b0^-256), 2b0^256, 10^(-fpprec+1));
+true;
+
+closeto(
+  gamma(2.0^-256 + 1d-200*%i),
+  rectform(1/(2.0^-256 + 1d-200*%i)),
+  2d-15);
+true;
+
+closeto(
+  gamma(2b0^-256 + 1b-500*%i),
+  rectform(1/(2b0^-256 + 1b-500*%i)),
+  2d-15);
+true;
+
 (fpprec:oldfpprec,done);
 done;
 

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

Summary of changes:
 src/csimp2.lisp       |  129 ++++++++++++++++++++++++++++++-------------------
 tests/rtest_gamma.mac |   23 +++++++++
 2 files changed, 103 insertions(+), 49 deletions(-)


hooks/post-receive
-- 
Maxima CAS