[Maxima] Complex Bigfloats for the Gamma function

Dieter Kaiser drdieterkaiser at web.de
Thu Sep 4 16:47:51 CDT 2008

```I have added support for Complex Bigfloat values to the Gamma function. With
this extension the Gamma function would have complete numerical support.

With Complex Bigfloats for the Gamma function it is possible to support complete
Complex Bigfloat evaluation for the Exponential Integral E(n,z) too.

I have rewritten the functions bbfac and cbfac in Lisp and called them
bfloat-factorial and complex-bfloat-factorial. What do you think? Should we use
the rewritten Lisp functions and not the Maxima functions?

I have tested the Gamma function with Bigfloats up to a precision of 64 digits
against functions.wolfram.com and against the original Maxima functions bffac
and cbffac.

Both GCL 2.6.8 and CLISP 2.46 have no problems. The testsuite runs without
errors.

Dieter Kaiser

Index: csimp2.lisp
===================================================================
RCS file: /cvsroot/maxima/maxima/src/csimp2.lisp,v
retrieving revision 1.21
diff -u -r1.21 csimp2.lisp
--- csimp2.lisp   14 Feb 2008 01:31:37 -0000   1.21
+++ csimp2.lisp   4 Sep 2008 21:18:38 -0000
@@ -180,12 +180,17 @@
(jr (\$realpart j))
(ji (\$imagpart j)))
(cond ((floatp j) (gammafloat j))
-     ((\$bfloatp j) (mfuncall '\$bffac (m+ j -1) \$fpprec))
+     ((\$bfloatp j) (bfloat-factorial (add j -1)))
((and (numberp jr)
(numberp ji)
(or \$numer (floatp jr) (floatp ji)))
(complexify (gamma-lanczos (complex (float jr)
(float ji)))))
+          ((and (mnump jr)
+                (mnump ji)
+                (or \$numer (\$bfloatp jr) (\$bfloatp ji)))
+           (complex-bfloat-factorial
+             (add (\$bfloat (sub jr 1)) (mul '\$%i (\$bfloat ji)))))
((or (not (mnump j))
(ratgreaterp (simpabs (list '(%abs) j) 1 t) \$gammalim))
(eqtest (list '(%gamma) j) x))
@@ -306,6 +311,74 @@
(defun gammafloat (a)
(realpart (gamma-lanczos (complex a 0.0))))

+;;; The numerical routine for Bigfloat evaluation of the Factorial function
+;;; is a translation of the routine bbfac written by Bill Gosper.
+
+(defun bfloat-factorial (z)
+  (let ((\$ratprint nil)
+        (bigfloat%pi  (\$bfloat '\$%pi)))
+    (cond
+      ((eq (\$sign z) '\$neg)
+       ;; for a negative argument we use the reflection formula
+       (div (div (mul bigfloat%pi z) (\$sin (mul bigfloat%pi z)))
+            (bfloat-factorial (mul -1 z))))
+      (t
+       (let* ((bigfloatone  (\$bfloat bigfloatone))
+              (bigfloathalf (div bigfloatone bigfloattwo))
+              (k (* 2 (+ 1 (\$entier (* 0.41 \$fpprec)))))
+              (m bigfloatone)
+              (y (power z+k bigfloattwo))
+              (x (\$bfloat bigfloatzero))
+              (ii))
+       (dotimes (i (/ k 2))
+         (setq ii (+ i 1))
+         (setq m (mul m (add z (* 2 ii) -1) (add z (* 2 ii))))
+         (setq x (div (add x
+                        (div (\$bern (+ k (* -2 ii) 2))
+                             (* (+ k (* -2 ii) 1) (+ k (* -2 ii) 2))))
+                      y)))
+       (div
+         (mul (power (mul bigfloattwo bigfloat%pi z+k) bigfloathalf)
+              (power (\$bfloat '\$%e) (mul z+k (add (\$log z+k) x -1))))
+       m))))))
+
+;;; The numerical routine for Complex Bigfloat evaluation of the Factorial
+;;; function is a translation of the routine cbbfac written by Bill Gosper.
+
+(defun complex-bfloat-factorial (z)
+  (let ((\$ratprint nil)
+        (bigfloat%pi  (\$bfloat '\$%pi)))
+  (cond
+    ((eq (\$sign (\$realpart z)) '\$neg)
+     ;; for a negative argument we use the reflection formula
+     (\$rectform (div (div (mul bigfloat%pi z) (\$sin (mul bigfloat%pi z)))
+                     (complex-bfloat-factorial (mul -1 z)))))
+    (t
+     (let* ((bigfloatone  (\$bfloat bigfloatone))
+            (bigfloathalf (div bigfloatone bigfloattwo))
+            (k (* 2 (+ 1 (\$entier (* 0.41 \$fpprec)))))
+            (m bigfloatone)
+            (y (\$rectform (power z+k bigfloattwo)))
+            (x (\$bfloat bigfloatzero))
+            (ii))
+       (dotimes (i (/ k 2))
+         (setq ii (+ i 1))
+         (setq m (\$rectform (mul m (add z (* 2 ii) -1) (add z (* 2 ii)))))
+         (setq x (\$rectform
+                          (div (\$bern (+ k (* -2 ii) 2))
+                               (* (+ k (* -2 ii) 1) (+ k (* -2 ii) 2))))
+                        y))))
+       (\$rectform
+         (div
+           (mul (power (mul bigfloattwo bigfloat%pi z+k) bigfloathalf)
+                (power (\$bfloat '\$%e) (mul z+k (add (\$log z+k) x -1))))
+           m)))))))
+
(declare-top (special \$numer \$trigsign))

(defmfun simperf (x vestigial z &aux y)

Success, CVS operation completed

```