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:

• Map values from 3rd or 4th quadrant to 1st or 2nd quadrant:

• For z = 0, and any negative integer, .

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.