Complex special functions

Algorithms for calculating gamma, log-gamma and digamma in the complex field

Author: bdeket

math-doc/math/scribblings/math-special-functions.scrbl

@@ -19,7 +19,7 @@
 The term ``special function'' has no formal definition. However, for the purposes of the
 @racketmodname[math] library, a @deftech{special function} is one that is not @tech{elementary}.
 
-The special functions are split into two groups: @secref{real-functions} and
+The special functions are split into three groups: @secref{complex-functions}, @secref{real-functions} and
 @secref{flonum-functions}. Functions that accept real arguments are usually defined
 in terms of their flonum counterparts, but are different in two crucial ways:
 @itemlist[
@@ -28,13 +28,13 @@ in terms of their flonum counterparts, but are different in two crucial ways:
        @racket[exn:fail:contract] instead of returning @racket[+nan.0].}
 ]
 
-Currently, @racketmodname[math/special-functions] does not export any functions that accept
-or return complex numbers. Mathematically, some of them could return complex numbers given
+Only functions in @secref{complex-functions} support accepting or returning complex arguments.
+Mathematically, some of the other functions could return complex numbers given
 real numbers, such @racket[hurwitz-zeta] when given a negative second argument. In
 these cases, they raise an @racket[exn:fail:contract] (for an exact argument) or return
 @racket[+nan.0] (for an inexact argument).
 
-Most real functions have more than one type, but they are documented as having only
+Most real and complex functions have more than one type, but they are documented as having only
 one. The documented type is the most general type, which is used to generate a contract for
 uses in untyped code. Use @racket[:print-type] to see all of a function's types.
 
@@ -56,11 +56,11 @@ The most general type @racket[Real -> (U Zero Flonum)] is used to generate
 @racket[lambert]'s contract when it is used in untyped code. Except for this discussion,
 this the only type documented for @racket[lambert].
 
-@section[#:tag "real-functions"]{Real Functions}
+@section[#:tag "complex-functions"]{Complex Functions}
 
-@defproc[(gamma [x Real]) (U Positive-Integer Flonum)]{
+@defproc[(gamma [x Number]) Number]{
 Computes the @hyperlink["http://en.wikipedia.org/wiki/Gamma_function"]{gamma function},
-a generalization of the factorial function to the entire real line, except nonpositive integers.
+a generalization of the factorial function to the entire complex plane, except nonpositive exact integers.
 When @racket[x] is an exact integer, @racket[(gamma x)] is exact.
 
 @examples[#:eval untyped-eval
@@ -76,16 +76,17 @@ When @racket[x] is an exact integer, @racket[(gamma x)] is exact.
                  (gamma -1.0)
                  (gamma 0.0)
                  (gamma -0.0)
+                 (gamma 1+i)
                  (gamma 172.0)
                  (eval:alts
                   (bf (gamma 172))
                   (eval:result @racketresultfont{(bf "1.241018070217667823424840524103103992618e309")}))]
 
-Error is no more than 10 @tech{ulps} everywhere that has been tested, and is usually no more than 4
-ulps.
+On the real line the error is no more than 10 @tech{ulps} everywhere that has been tested, and is usually no more than 4
+ulps. In the rest of the complex plane the relative error is smaller than 1e-13 (@tech{ulps}=450).
 }
 
-@defproc[(log-gamma [x Real]) (U Zero Flonum)]{
+@defproc[(log-gamma [x Number]) Number]{
 Like @racket[(log (abs (gamma x)))], but more accurate and without unnecessary overflow.
 The only exact cases are @racket[(log-gamma 1) = 0] and @racket[(log-gamma 2) = 0].
 
@@ -99,25 +100,27 @@ The only exact cases are @racket[(log-gamma 1) = 0] and @racket[(log-gamma 2) =
                  (log-gamma -1)
                  (log-gamma -1.0)
                  (log-gamma 0.0)
+                 (log-gamma 1+i)
                  (log (abs (gamma 172.0)))
                  (log-gamma 172.0)]
 
-Error is no more than 11 @tech{ulps} everywhere that has been tested, and is usually no more than 2
+On the real line error is no more than 11 @tech{ulps} everywhere that has been tested, and is usually no more than 2
 ulps. Error reaches its maximum near negative roots.
 }
 
-@defproc[(psi0 [x Real]) Flonum]{
+@defproc[(psi0 [x Number]) Number]{
 Computes the @hyperlink["http://en.wikipedia.org/wiki/Digamma_function"]{digamma function},
 the logarithmic derivative of the gamma function.
 
 @examples[#:eval untyped-eval
                  (plot (function psi0 -2.5 4.5) #:y-min -5 #:y-max 5)
                  (psi0 0)
                  (psi0 1)
+                 (psi0 1+i)
                  (- gamma.0)]
 
-Except near negative roots, maximum observed error is 2 @tech{ulps}, but is usually no more
-than 1.
+On the real line, except near negative roots, maximum observed error is 2 @tech{ulps}, but is usually no more
+than 1. In the complex plane the relative error is around 1e-13.
 
 Near negative roots, which occur singly between each pair of negative integers, @racket[psi0]
 exhibits @tech{catastrophic cancellation} from using the reflection formula, meaning that
@@ -129,6 +132,8 @@ error.
 If you need low relative error near negative roots, use @racket[bfpsi0].
 }
 
+@section[#:tag "real-functions"]{Real Functions}
+
 @defproc[(psi [m Integer] [x Real]) Flonum]{
 Computes a @hyperlink["http://en.wikipedia.org/wiki/Polygamma_function"]{polygamma function},
 or the @racket[m]th logarithmic derivative of the gamma function. The order @racket[m] must be

math-lib/math/private/functions/gamma.rkt

@@ -272,13 +272,38 @@ Approximations:
                   [(and (x . fl> . -4.5) (x . fl< . 4.5))  (flgamma-taylor x)]
                   [else  (flgamma-lanczos x)])))]))
 
+(define √2π 2.5066282746310005024157652848110)
+(: complex-gamma (Number -> Number))
+(define (complex-gamma z+1)
+  (define neg? (< (real-part z+1)0))
+  (define z (- (if neg? (- z+1) z+1) 1))
+  (cond
+    [(or (= z 1)(= z 0))1]
+    [else
+     (define zh (+ z 0.5))
+     (define zgh (+ zh lanczos-complex-g))
+     (define zp (expt zgh (* 0.5 zh)))
+
+     (define ans
+       (* zp (exp (- zgh)) zp
+          √2π (+ (car lanczos-complex-c)
+                 (for/sum : Number ([a (in-list (cdr lanczos-complex-c))]
+                                    [i (in-naturals 1)])
+                   (/ a (+ z i))))))
+     (if neg?
+         (- (/ pi (* ans z+1 (sin (* pi z+1)))))
+         ans)]))
+
 (: gamma (case-> (One -> One)
                  (Integer -> Positive-Integer)
                  (Float -> Float)
-                 (Real -> (U Positive-Integer Flonum))))
-(define (gamma x)
+                 (Real -> (U Positive-Integer Flonum))
+                 (Number -> Number)))
+(define (gamma z)
+  (define x (if (= (imag-part z) 0) (real-part z) z))
   (cond [(double-flonum? x)  (flgamma x)]
         [(exact-integer? x)
          (cond [(x . > . 0)  (factorial (- x 1))]
                [else  (raise-argument-error 'gamma "Real, not Zero or Negative-Integer" x)])]
-        [else  (flgamma (fl x))]))
+        [(real? x)  (flgamma (fl x))]
+        [else (complex-gamma z)]))

math-lib/math/private/functions/lanczos.rkt

@@ -2,7 +2,7 @@
 
 (require "../../flonum.rkt")
 
-(provide lanczos-sum lanczos-g)
+(provide (all-defined-out))
 
 ;; Lanczos polynomial for N=13 G=6.024680040776729583740234375
 ;; Max experimental error (with arbitary precision arithmetic) is 1.196214e-17
@@ -36,3 +36,25 @@
     1.0)))
 
 (define lanczos-g 6.024680040776729583740234375)
+
+
+;; Lanczos series approximation for the complex plane
+;; ref Paul Godfrey's MATLAB implementations | https://my.fit.edu/~gabdo/paulbio.html
+;; ->relative error < 1e-13
+(define lanczos-complex-c
+  (list    0.99999999999999709182
+          57.156235665862923517
+         -59.597960355475491248
+          14.136097974741747174
+          -0.49191381609762019978
+           0.33994649984811888699e-4
+           0.46523628927048575665e-4
+          -0.98374475304879564677e-4
+           0.15808870322491248884e-3
+          -0.21026444172410488319e-3
+           0.21743961811521264320e-3
+          -0.16431810653676389022e-3
+           0.84418223983852743293e-4
+          -0.26190838401581408670e-4
+           0.36899182659531622704e-5))
+(define lanczos-complex-g 607/128)

math-lib/math/private/functions/log-gamma.rkt

@@ -4,7 +4,8 @@
          "../../flonum.rkt"
          "../../base.rkt"
          "log-gamma-zeros.rkt"
-         "gamma.rkt")
+         "gamma.rkt"
+         "lanczos.rkt")
 
 (provide fllog-gamma log-gamma)
 
@@ -398,15 +399,43 @@
                   [(x . fl< . 4.5)  (fllog-gamma-small-positive x)]
                   [else  (fllog-gamma-large-positive x)])))]))
 
+(define log-√2π 0.9189385332046727417803297)
+(define log-π 1.14472988584940017414342735)
+(: complex-log-gamma (Number -> Number))
+(define (complex-log-gamma z)
+  (define neg? (< (real-part z)0))
+  (define z* (if neg? (- z) z))
+  (cond
+    [(or (= z* 1)(= z* 2))0]
+    [else
+     (define z- (- z* 0.5))
+     (define zg (+ z- lanczos-complex-g))
+     
+     (define ans
+       (+ (- zg) (* z- (log zg))
+          log-√2π (log (+ (car lanczos-complex-c)
+                          (for/sum : Number ([a (in-list (cdr lanczos-complex-c))]
+                                             [i (in-naturals 0)])
+                            (/ a (+ z* i)))))))
+     (cond
+       [neg?
+        (define lpi (+ log-π (* (if (< (imag-part z) 0) +i -i) pi)))
+        (- lpi (log z) ans (log (sin (* pi z))))]
+       [else ans])]))
+
+
 (: log-gamma (case-> (One -> Zero)
                      (Flonum -> Flonum)
-                     (Real -> (U Zero Flonum))))
-(define (log-gamma x)
+                     (Real -> (U Zero Flonum))
+                     (Number -> Number)))
+(define (log-gamma z)
+  (define x (if (= (imag-part z) 0) (real-part z) z))
   (cond [(flonum? x)  (fllog-gamma x)]
         [(single-flonum? x)  (fllog-gamma (fl x))]
         [(integer? x)
          (cond [(x . <= . 0)
                 (raise-argument-error 'log-gamma "Real, not Zero or Negative-Integer" x)]
                [(eqv? x 1)  0]
                [else  (fllog-factorial (fl (- x 1)))])]
-        [else  (fllog-gamma (fl x))]))
+        [(real? x)  (fllog-gamma (fl x))]
+        [else (complex-log-gamma z)]))

math-lib/math/private/functions/psi.rkt

@@ -6,7 +6,8 @@
          "../number-theory/bernoulli.rkt"
          "gamma.rkt"
          "hurwitz-zeta.rkt"
-         "tan-diff.rkt")
+         "tan-diff.rkt"
+         "lanczos.rkt")
 
 (provide flpsi0 flpsi psi0 psi)
 
@@ -208,6 +209,29 @@
         [(x . = . +inf.0)  +inf.0]
         [else  +nan.0]))
 
+(: complex-psi0 (Number -> Number))
+(define (complex-psi0 z*)
+  (define neg? (< (real-part z*) 0.5))
+  (define z (if neg? (- 1 z*) z*))
+
+  (define-values (n d)
+    (for/fold : (Values Number Number)
+      ([n : Number 0][d : Number 0])
+      ([c (in-list (reverse (cdr lanczos-complex-c)))]
+       [k (in-range (length lanczos-complex-c) -1 -1)])
+      (define dz (/ 1 (+ z k -2)))
+      (define dd (* c dz))
+      (values (- n (* dd dz)) (+ d dd))))
+
+  (define gg (+ z lanczos-complex-g -0.5))
+
+  (define ans (+ (log gg) (- (/ n (+ d (car lanczos-complex-c)))
+                             (/ lanczos-complex-g gg))))
+  
+  (if neg?
+      (- ans (/ pi (tan (* pi z*))))
+      ans))
+
 (define pi.128 267257146016241686964920093290467695825/85070591730234615865843651857942052864)
 
 (: flexppi (Flonum -> Flonum))
@@ -232,13 +256,16 @@
          (fl- (fl* sgn (flpsi m (fl- 1.0 x))) t)]
         [else  +nan.0]))
 
-(: psi0 (Real -> Flonum))
-(define (psi0 x)
+(: psi0 (case-> (Real -> Flonum)
+                (Number -> Number)))
+(define (psi0 z)
+  (define x (if (= (imag-part z) 0) (real-part z) z))
   (cond [(flonum? x)  (flpsi0 x)]
         [(single-flonum? x)  (flpsi0 (fl x))]
         [(and (integer? x) (x . <= . 0))
          (raise-argument-error 'psi0 "Real, not Zero or Negative-Integer" x)]
-        [else  (flpsi0 (fl x))]))
+        [(real? x)  (flpsi0 (fl x))]
+        [else (complex-psi0 z)]))
 
 (: psi (Integer Real -> Flonum))
 (define (psi m x)

math-lib/math/special-functions.rkt

@@ -9,19 +9,23 @@
          (except-in "private/functions/lambert.rkt" lambert)
          (except-in "private/functions/zeta.rkt" eta zeta)
          (except-in "private/functions/hurwitz-zeta.rkt" hurwitz-zeta)
-         "private/functions/psi.rkt"
+         (except-in "private/functions/psi.rkt" psi0)
          "private/functions/incomplete-gamma.rkt"
          "private/functions/incomplete-beta.rkt"
          "private/functions/stirling-error.rkt"
          "private/functions/fresnel.rkt")
 
 (require/untyped-contract
  "private/functions/gamma.rkt"
- [gamma  (Real -> Real)])
+ [gamma  (Number -> Number)])
 
 (require/untyped-contract
  "private/functions/log-gamma.rkt"
- [log-gamma  (Real -> Real)])
+ [log-gamma  (Number -> Number)])
+
+(require/untyped-contract
+ "private/functions/psi.rkt"
+ [psi0  (Number -> Number)])
 
 (require/untyped-contract
  "private/functions/beta.rkt"
@@ -61,6 +65,7 @@
           "private/functions/fresnel.rkt")
          gamma
          log-gamma
+         psi0
          beta log-beta
          erf erfc
          lambert

math-test/math/tests/function-tests.rkt

@@ -46,7 +46,31 @@
 (check-exn exn:fail:contract? (λ () (gamma 0)))
 (check-equal? (gamma 4) 6)
 (check-equal? (gamma 4.0) (flgamma 4.0))
-
+;checks calculated at https://wolframalpha.com/input/?i=gamma(z) / log-gamma(z) / digamma(z)
+(let ([z 1.5+0.2i]
+      [ans 0.869862133805271779271655350309163858030873370542544143+0.00730170831118402614972641679459165642341675697339983422i])
+  (check-= (/ (gamma z) ans) 1 1e-13))
+(let ([z 1.5+4i]
+      [ans -0.018686297879039387689174933326088586248898218845351961+0.0026241318932335593901663964582546948097937382268327333i])
+  (check-= (/ (gamma z) ans) 1 1e-13))
+(let ([z 196+285i]
+      [ans 5.72038483529989178989490686567823916571020691258673e291-2.00024207656771597752271221062051820963168450311374e291i])
+  (check-= (/ (gamma z) ans) 1 1e-13))
+(let ([z 28.45-0.05i]
+      [ans 4.788665109060112126493799529515601713890011469557001e28-8.048791086601393564384851295197107978228121888100461e27i])
+  (check-= (/ (gamma z) ans) 1 1e-13))
+(let ([z 0.5-0.00333i]
+      [ans 1.772367471310459385811039514905244407455502656029765525+0.01158858570803772503302975596977297632272041615128277137i])
+  (check-= (/ (gamma z) ans) 1 1e-13))
+(let ([z -0.5-0.00333i]
+      [ans -3.544732073864574326606517758100934969124507449000061+0.0004307441958626149491398963294062742489283873077748653i])
+  (check-= (/ (gamma z) ans) 1 1e-13))
+(let ([z -1+i]
+      [ans -0.9974967895818289935938328922330359491032588243198502190-4.228631137454774745965835886896275145906361523162647513i])
+  (check-= (/ (log-gamma z) ans) 1 1e-13))
+(let ([z -1+i]
+      [ans 0.59465032062247697727187848272191072247626297176354162323+2.5766740474685811741340507947500004904456562664038166655i])
+  (check-= (/ (psi0 z) ans) 1 1e-13))
 ;; ---------------------------------------------------------------------------------------------------
 ;; hypot