Mathematics
f(z)    Complex Numbers and Functions     Tech Note
Complex Gamma Function

 This page shows how to compute the gamma function for a complex argument.  This computation first involves finding the natural logarithm of the complex gamma function.

CLnGamma:  Complex Log Gamma function

The following asymptotic formula is applied over the ranges explained below:

This formula is from the Table of the Gamma Function for Complex Arguments (U.S. National Bureau of Standards Applied Mathematics Series No. 34), August 1954. p. VII.  The Method of Computation section (p. VIII-X) mentioned this formula was first computed to 15 decimal places in the rectangular region 9 £ x £ 10 and 0 £ y £ 10.  The terms 1/zm were retained for up to m = 13, inclusive, at which point the truncating error did not exceed a half-unit in the 15th place.

A truncated version of the above formula is seen as Formula 6.1.41 in Handbook of Mathematical Functions by Abramowitz and Stegun, p. 257.  Another source of information is Mathematical Tables of ln G(z) for Complex Argument by A. A. Abramov (translated from the Russian), Pergamon Press, 1960.

For values “below” the range, 9 £ x £ 10, the following recurrence relationship was used to “start” with a value in that asymptotic range (this includes mapping values from the 2nd quadrant to the 1st quadrant):

This equation is derived from the more common recurrence formula:

Two other “rules” include:

CAsymptoticLnGamma and CLnGamma functions

// This is the asymptotic approximation used in the National Bureau of
// Standards "Table of the Gamma Function for Complex Arguments,"
// Applied Mathematics Series 34, 1954.  The NBS table was created
// using this approximation over the area  9 <= Re(z) <= 10 and
// 0 <= Im(z) <= 10.  Other table values were computed using the
// relationship ln gamma(z+1) = ln z + ln gamma(z).
FUNCTION CAsymptoticLnGamma (CONST a:  TComplex):  TComplex;
  CONST
    c:  ARRAY[1..8] OF TReal
     =  (1/12, -1/360, 1/1260, -1/1680, 1/1188, -691/360360,
         1/156, -3617/122400);
  VAR
    aTemp :  TComplex;
    i     :  WORD;
    powers:  ARRAY[1..8] OF TComplex;
    temp1 :  TComplex;
    temp2 :  TComplex;
BEGIN
  aTemp := CConvert (a,cfRectangular);
  temp1 := CLn(aTemp);                           // ln(z)
  temp2 := CSet (aTemp.x-0.5, aTemp.y, cfRectangular);   // a - 0.5
  RESULT := CMult(temp1, temp2);                 // (a - 0.5)*ln(a)
  RESULT := CSub(RESULT, aTemp);                 // (a - 0.5)*ln(a) - a
  RESULT.x := RESULT.x + HalfLn2PI;
 
  temp1 := ComplexOne;
  powers[1] := CDiv(temp1, aTemp);              // a^-1
  temp2 := CMult(powers[1],powers[1]);          // a^-2
  FOR i := 2 TO 8 DO
    powers[i] := CMult(powers[i-1],temp2);
 
  FOR i := 8 DOWNTO 1 DO                        // add in reverse order
  BEGIN
    temp1 := CSet(c[i]*powers[i].x, c[i]*powers[i].y, cfRectangular);
    RESULT := CAdd(RESULT, temp1);
  END
END {CAsymptoticLnGamma};
 
 
FUNCTION CLnGamma (CONST a:  TComplex):  TComplex;
  VAR
    aTemp:  TComplex;
    lna  :  TComplex;
    temp :  TComplex;
BEGIN
  aTemp := CConvert (a, cfRectangular);
 
  IF   (aTemp.x <= 0.0) AND (DeFuzz(aTemp.y) = 0.0)
  THEN BEGIN
    IF   DeFuzz( Abs(Frac(aTemp.x)) ) = 0.0     // negative integer?
    THEN BEGIN
      RESULT := CSet(PositiveInfinity, PositiveInfinity);
      EXIT
    END
  END;
 
  IF   aTemp.y < 0.0                 // 3rd or 4th quadrant?
  THEN BEGIN
    temp := CConjugate(aTemp);
    RESULT := CLnGamma(temp);        // try again in 1st or 2nd quadrant
    RESULT := CConjugate(RESULT)     // left this out!  1/3/91
  END
  ELSE BEGIN
    IF   aTemp.x < 9.0               // "left" of NBS table range
    THEN BEGIN
      lna := CLn(aTemp);
      temp := CSet (aTemp.x+1.0, aTemp.y, cfRectangular);
      temp := CLnGamma(temp);        // recursive call
      RESULT := CSub(temp, lna)
    END
    ELSE RESULT := CAsymptoticLnGamma(aTemp) // NBS range:  9 <= Re(z) <= 10
  END
END {CLnGamma};

Selected computed values that match published values

Complex LnGamma function CLnGamma = Ln[Gamma(z)]
 
                          Ln Gamma(z)
     z                    rectangular
------------  ----------------------------------
 0.0 -  1.0i   -0.650923199302 + 1.872436647262i 
 0.0 +  1.0i   -0.650923199302 - 1.872436647262i
 0.0 -  3.0i   -4.342756588258 + 0.517445555726i
 0.0 +  3.0i   -4.342756588258 - 0.517445555726i
 
 0.5 +  0.5i    0.112387242810 - 0.750729202121i
 
 1.0 +  0.0i    0.000000000000 + 0.000000000000i
 1.0 +  1.0i   -0.650923199302 - 0.301640320468i
 1.0 -  1.0i   -0.650923199302 + 0.301640320468i
 
 5.0 +  0.0i    3.178053830348 + 0.000000000000i
 5.0 +  3.0i    2.244246717020 + 4.714089538905i

Note: 

CGamma:  Complex Gamma function

The CGamma function involves both the CLnGamma (from above) and CExp functions:

CGamma function

FUNCTION CGamma (CONST a:  TComplex):  TComplex;
  VAR
    lna:  TComplex;
BEGIN
  lna := CLnGamma(a);
  IF   lna.x >= 75.0    // arbitrary cutoff for infinity
  THEN RESULT := CSet(PositiveInfinity, PositiveInfinity)
  ELSE
    IF   lna.x < -200.0 // arbitrary cutoff for underflow
    THEN RESULT := ComplexZero        // avoid underflow
    ELSE RESULT := CExp(lna);
END {CGamma};

This code calls the CGamma function for various complex values, a[k].

Memo.Lines.Add('Complex Gamma function CGamma = Gamma(z)');
Memo.Lines.Add('');
Memo.Lines.Add('                             Gamma(z)           ' +
               '             Gamma(z)');
Memo.Lines.Add('          z                rectangular          ' +
               '              polar');
Memo.Lines.Add('     ------------  ----------------------------   ' +
               '--------------------------------');
 
FOR k := Low(a) TO High(a) DO
BEGIN
  z := CGamma(a[k]);
 
  Memo.Lines.Add(   Format('%2d  %-s %-s %-s',
                            [k, CToRectStr(a[k],  5,1),
                                CToRectStr(z, 13,8),
                                CToPolarStr(z,  13,8)]))
END;

The code above puts its result in a memo box and looks like the following:

Sample computed values

Complex Gamma function CGamma = Gamma(z)
 
                             Gamma(z)                        Gamma(z)
          z                rectangular                        polar
     ------------  ----------------------------   --------------------------------
 1    0.0 +  0.0i           INF +          INFi           INF * CIS(   0.00000000)
 2    0.5 +  0.5i    0.81816400 -   0.76331383i    1.11894608 * CIS(  -0.75072920)
 3   -0.5 +  0.5i   -1.58147783 -   0.05485017i    1.58242872 * CIS(  -3.10692369)
 4   -0.5 -  0.5i   -1.58147783 +   0.05485017i    1.58242872 * CIS(   3.10692369)
 5    0.5 -  0.5i    0.81816400 +   0.76331383i    1.11894608 * CIS(   0.75072920)
 6    1.0 +  0.0i    1.00000000 +   0.00000000i    1.00000000 * CIS(   0.00000000)
 7    1.0 +  1.0i    0.49801567 -   0.15494983i    0.52156405 * CIS(  -0.30164032)
 8    0.0 +  1.0i   -0.15494983 -   0.49801567i    0.52156405 * CIS(  -1.87243665)
 9   -1.0 +  1.0i   -0.17153292 +   0.32648275i    0.36880147 * CIS(   2.05455417)
10   -1.0 +  0.0i           INF +          INFi           INF * CIS(   0.00000000)
11   -1.0 -  1.0i   -0.17153292 -   0.32648275i    0.36880147 * CIS(  -2.05455417)
12    0.0 -  1.0i   -0.15494983 +   0.49801567i    0.52156405 * CIS(   1.87243665)
13    1.0 -  1.0i    0.49801567 +   0.15494983i    0.52156405 * CIS(   0.30164032)
14    5.0 +  0.0i   24.00000000 +   0.00000000i   24.00000000 * CIS(   0.00000000)
15    5.0 +  3.0i    0.01604188 -   9.43329329i    9.43330693 * CIS(  -1.56909577)
16    0.0 +  3.0i    0.01129867 -   0.00643092i    0.01300064 * CIS(  -0.51744556)
17   -5.0 +  3.0i    0.00000790 +   0.00000476i    0.00000922 * CIS(   0.54212006)
18   -5.0 +  0.0i           INF +          INFi           INF * CIS(   0.00000000)
19   -5.0 -  3.0i    0.00000790 -   0.00000476i    0.00000922 * CIS(  -0.54212006)
20    0.0 -  3.0i    0.01129867 +   0.00643092i    0.01300064 * CIS(   0.51744556)
21   -5.0 -  3.0i    0.00000790 -   0.00000476i    0.00000922 * CIS(  -0.54212006)

Notes:

·        G(1 + 0i) = 0! = 1

·        G(5 + 0i) = 4! = 24

·        G(0.5) = =1.77245385…

·        See Tables of Functions with Formulae and Curves by Eugen Jahnke and Fritze Emde, Dover Publications, 1945, pp. 12-13 for contour diagram and 3D relief view of the gamma function.

 

Download ComplexMathLibrary

Updated 18 Feb 2002


since 24 June 2001