While playing around with some elliptic integrals, I needed to
evaluate some gamma functions. The numerical gamma functions included
with maxima seem to have about 5-6 digits of accuracy.
I've included an updated version that is supposed to accurate to about
15 digits and is valid for the entire complex plane. (The current
code doesn't support numerically evaluating the gamma function at
complex values.)
To use this, just replace gammafloat in csimp2 with the following
code. You can probably also delete gamma and do-gammafloat in
csimp2.lisp.
Ray
;; This implementation is based on Lanczos convergent formula for the
;; gamma function for Re(z) > 0. We can use the reflection formula
;;
;; -z*Gamma(z)*Gamma(-z) = pi/sin(pi*z)
;;
;; to handle the case of Re(z) <= 0.
;;
;; See http://winnie.fit.edu/~gabdo/gamma.m for some matlab code to
;; compute this and http://winnie.fit.edu/~gabdo/gamma.txt for a nice
;; discussion of Lanczos method and an improvement of Lanczos method.
;;
;;
;; The document says this should give about 15 digits of accuracy for
;; double-precision IEEE floats. The document also indicates how to
;; compute a new set of coefficients if you need more range or
;; accuracy.
(defun gamma-lanczos (z)
(declare (type (complex double-float) z)
(let ((g 607/128)
(c (make-array 15 :element-type 'double-float
:initial-contents
'(0.99999999999999709182d0
57.156235665862923517d0
-59.597960355475491248d0
14.136097974741747174d0
-0.49191381609762019978d0
.33994649984811888699d-4
.46523628927048575665d-4
-.98374475304879564677d-4
.15808870322491248884d-3
-.21026444172410488319d-3
.21743961811521264320d-3
-.16431810653676389022d-3
.84418223983852743293d-4
-.26190838401581408670d-4
.36899182659531622704d-5))))
(declare (type (rational 607/128 607/128) g)
(type (simple-array double-float (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? What
(/ pi (* (- z) (sin (* pi z)) (gamma-lanczos (- z))))
(let* ((z (- z 1))
(zh (+ z 1/2))
(zgh (+ zh 607/128))
(zp (expt zgh (/ zh 2)))
(ss
(loop for pp of-type fixnum from (1- (length c)) downto 1
sum (/ (aref c pp) (+ z pp)))))
(* (sqrt (* 2 pi))
(+ ss (aref c 0))
(* zp (exp (- zgh)) zp))))))
(defun gammafloat (a)
(realpart (gamma-lanczos (complex a 0d0))))