The gamma function of a complex ball

acb_hypgeom_gamma, acb_hypgeom_rgamma and acb_hypgeom_lgamma compute $\Gamma(z)$, $1/\Gamma(z)$ and the principal branch of $\log \Gamma(z)$ for a complex ball $z$. They are the implementations behind acb_gamma, acb_rgamma and acb_lgamma. Real arguments are passed to the real functions (see The gamma function of a real ball); otherwise the functions use a tabulated Taylor series near the real axis and the Stirling series elsewhere, as described in The gamma function.

The principal branch of log-gamma#

The function computed by acb_hypgeom_lgamma is not $\log(\Gamma(z))$ with the principal logarithm. It is the analytic continuation of $\log\Gamma(x)$ from the positive real axis to $\mathbb{C} \setminus (-\infty, 0]$, which is the convention of most computer algebra systems and of the DLMF [NIST2012, §5.4]: it satisfies $\log\Gamma(z+1) = \log\Gamma(z) + \log z$ exactly (not only modulo $2\pi i$), and its imaginary part grows like $\operatorname{Im}(z) \log |z|$ instead of staying in $(-\pi, \pi]$. The two functions differ by $2\pi i k(z)$, where the integer $k(z)$ changes along curves where $\operatorname{Im}\log\Gamma(z)$ is an odd multiple of $\pi$. On the branch cut, the value is the limit from the upper half-plane.

Computing this branch requires more than a logarithm of $\Gamma(z)$:

• The Stirling series gives $\log\Gamma(z)$ directly for $z$ of large modulus in the right half-plane. Its argument reduction $\log\Gamma(z) = \log\Gamma(z+r) - \log\,(z)_r$ needs the log-rising factorial $\log\,(z)_r = \sum_{k<r} \log(z+k)$ with the correct multiple of $2\pi i$, which acb_hypgeom_log_rising_ui computes with Hare's algorithm [Hare1997].
• The reflection formula for the left half-plane becomes $\log\Gamma(z) = \log\pi - \log\sin(\pi z) - \log\Gamma(1-z)$, where $\log\sin(\pi z)$ must be the analytic branch computed by acb_log_sin_pi.
• Near the real axis at moderate precision, the Taylor series of $1/\Gamma(1+u)$ is used for $u$ with real part in about $[-1/2, 1/2]$; the branch of $\log$ applied to it is fixed with a table of the lines where $\operatorname{Im}\log\Gamma$ crosses $(k + \tfrac12)\pi$.

Example#

At $z = -10.25 + 0.5i$ the principal branch of $\log\Gamma$ and the principal logarithm of $\Gamma$ differ by $-10\pi i$:

#include <flint/acb.h>
#include <flint/acb_hypgeom.h>

int main(void)
{
    acb_t z, y, t;

    acb_init(z);
    acb_init(y);
    acb_init(t);

    acb_set_d_d(z, -10.25, 0.5);

    acb_hypgeom_gamma(y, z, 64);
    flint_printf("gamma(z)      = "); acb_printn(y, 15, 0); flint_printf("\n");

    acb_hypgeom_lgamma(y, z, 64);
    flint_printf("lgamma(z)     = "); acb_printn(y, 15, 0); flint_printf("\n");

    acb_hypgeom_gamma(t, z, 64);
    acb_log(t, t, 64);
    flint_printf("log(gamma(z)) = "); acb_printn(t, 15, 0); flint_printf("\n");

    /* the difference is 2 pi i k */
    acb_sub(t, y, t, 64);
    acb_const_pi(y, 64);
    acb_mul_2exp_si(y, y, 1);
    acb_div(t, t, y, 64);
    flint_printf("difference / (2 pi) = "); acb_printn(t, 15, 0); flint_printf("\n");

    acb_hypgeom_rgamma(y, z, 64);
    flint_printf("rgamma(z)     = "); acb_printn(y, 15, 0); flint_printf("\n");

    acb_clear(z);
    acb_clear(y);
    acb_clear(t);
    flint_cleanup();
    return 0;
}

gamma(z)      = [7.08331784482241e-8 +/- 3.66e-23] + [-1.88626361829666e-7 +/- 2.32e-22]*I
lgamma(z)     = [-15.4175381117796 +/- 4.86e-15] + [-32.6274954179361 +/- 2.24e-14]*I
log(gamma(z)) = [-15.4175381117796 +/- 4.86e-15] + [-1.21156888203815 +/- 4.69e-15]*I
difference / (2 pi) = [+/- 2.62e-19] + [-5.00000000000000 +/- 1.09e-18]*I
rgamma(z)     = [1744778.02581992 +/- 4.92e-9] + [4646284.95319217 +/- 4.09e-9]*I

$\Gamma(z)$ is small here because $z$ is close to the negative real axis between two poles; $1/\Gamma(z)$, which is entire, is large.