acb_hypgeom.h – hypergeometric functions in the complex numbers

The generalized hypergeometric function is formally defined by

\[{}_pF_q(a_1,\ldots,a_p;b_1,\ldots,b_q;z) = \sum_{k=0}^\infty \frac{(a_1)_k\dots(a_p)_k}{(b_1)_k\dots(b_q)_k} \frac {z^k} {k!}.\]

It can be interpreted using analytic continuation or regularization when the sum does not converge. In a looser sense, we understand “hypergeometric functions” to be linear combinations of generalized hypergeometric functions with prefactors that are products of exponentials, powers, and gamma functions.

Convergent series

In this section, we define

\[T(k) = \frac{\prod_{i=0}^{p-1} (a_i)_k}{\prod_{i=0}^{q-1} (b_i)_k} z^k\]

and

\[{}_pH_{q}(a_0,\ldots,a_{p-1}; b_0 \ldots b_{q-1}; z) = {}_{p+1}F_{q}(a_0,\ldots,a_{p-1},1; b_0 \ldots b_{q-1}; z) = \sum_{k=0}^{\infty} T(k)\]

For the conventional generalized hypergeometric function \({}_pF_{q}\), compute \({}_pH_{q+1}\) with the explicit parameter \(b_q = 1\), or remove a 1 from the \(a_i\) parameters if there is one.

void acb_hypgeom_pfq_bound_factor(mag_t C, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, ulong n)

Computes a factor C such that

\[\left|\sum_{k=n}^{\infty} T(k)\right| \le C |T(n)|.\]

We check that \(\operatorname{Re}(b+n) > 0\) for all lower parameters b. If this does not hold, C is set to infinity. Otherwise, we cancel out pairs of parameters \(a\) and \(b\) against each other. We have

\[\left|\frac{a+k}{b+k}\right| = \left|1 + \frac{a-b}{b+k}\right| \le 1 + \frac{|a-b|}{|b+n|}\]

and

\[\left|\frac{1}{b+k}\right| \le \frac{1}{|b+n|}\]

for all \(k \ge n\). This gives us a constant D such that \(T(k+1) \le D T(k)\) for all \(k \ge n\). If \(D \ge 1\), we set C to infinity. Otherwise, we take \(C = \sum_{k=0}^{\infty} D^k = (1-D)^{-1}\).

As currently implemented, the bound becomes infinite when \(n\) is too small, even if the series converges.

long acb_hypgeom_pfq_choose_n(acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long prec)

Heuristically attempts to choose a number of terms n to sum of a hypergeometric series at a working precision of prec bits.

Uses double precision arithmetic internally. As currently implemented, it can fail to produce a good result if the parameters are extremely large or extremely close to nonpositive integers.

Numerical cancellation is assumed to be significant, so truncation is done when the current term is prec bits smaller than the largest encountered term.

This function will also attempt to pick a reasonable truncation point for divergent series.

void acb_hypgeom_pfq_sum_forward(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)
void acb_hypgeom_pfq_sum_rs(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)
void acb_hypgeom_pfq_sum_bs(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)
void acb_hypgeom_pfq_sum_fme(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)
void acb_hypgeom_pfq_sum(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)

Computes \(s = \sum_{k=0}^{n-1} T(k)\) and \(t = T(n)\). Does not allow aliasing between input and output variables. We require \(n \ge 0\).

The forward version computes the sum using forward recurrence.

The bs version computes the sum using binary splitting.

The rs version computes the sum in reverse order using rectangular splitting. It only computes a magnitude bound for the value of t.

The fme version uses fast multipoint evaluation.

The default version automatically chooses an algorithm depending on the inputs.

void acb_hypgeom_pfq_sum_bs_invz(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t w, long n, long prec)
void acb_hypgeom_pfq_sum_invz(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, const acb_t w, long n, long prec)

Like acb_hypgeom_pfq_sum(), but taking advantage of \(w = 1/z\) possibly having few bits.

void acb_hypgeom_pfq_direct(acb_t res, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)

Computes

\[{}_pH_{q}(z) = \sum_{k=0}^{\infty} T(k) = \sum_{k=0}^{n-1} T(k) + \varepsilon\]

directly from the defining series, including a rigorous bound for the truncation error \(\varepsilon\) in the output.

If \(n < 0\), this function chooses a number of terms automatically using acb_hypgeom_pfq_choose_n().

void acb_hypgeom_pfq_series_direct(acb_poly_t res, const acb_poly_struct * a, long p, const acb_poly_struct * b, long q, const acb_poly_t z, int regularized, long n, long len, long prec)

Computes \({}_pH_{q}(z)\) directly using the defining series, given parameters and argument that are power series. The result is a power series of length len.

An error bound is computed automatically as a function of the number of terms n. If \(n < 0\), the number of terms is chosen automatically.

If regularized is set, the regularized hypergeometric function is computed instead.

Asymptotic series

Let \(U(a,b,z)\) denote the confluent hypergeometric function of the second kind with the principal branch cut, and let \(U^{*} = z^a U(a,b,z)\). For all \(z \ne 0\) and \(b \notin \mathbb{Z}\) (but valid for all \(b\) as a limit), we have (DLMF 13.2.42)

\[U(a,b,z) = \frac{\Gamma(1-b)}{\Gamma(a-b+1)} M(a,b,z) + \frac{\Gamma(b-1)}{\Gamma(a)} z^{1-b} M(a-b+1,2-b,z).\]

Moreover, for all \(z \ne 0\) we have

\[\frac{{}_1F_1(a,b,z)}{\Gamma(b)} = \frac{(-z)^{-a}}{\Gamma(b-a)} U^{*}(a,b,z) + \frac{z^{a-b} e^z}{\Gamma(a)} U^{*}(b-a,b,-z)\]

which is equivalent to DLMF 13.2.41 (but simpler in form).

We have the asymptotic expansion

\[U^{*}(a,b,z) \sim {}_2F_0(a, a-b+1, -1/z)\]

where \({}_2F_0(a,b,z)\) denotes a formal hypergeometric series, i.e.

\[U^{*}(a,b,z) = \sum_{k=0}^{n-1} \frac{(a)_k (a-b+1)_k}{k! (-z)^k} + \varepsilon_n(z).\]

The error term \(\varepsilon_n(z)\) is bounded according to DLMF 13.7. A case distinction is made depending on whether \(z\) lies in one of three regions which we index by \(R\). Our formula for the error bound increases with the value of \(R\), so we can always choose the larger out of two indices if \(z\) lies in the union of two regions.

Let \(r = |b-2a|\). If \(\operatorname{Re}(z) \ge r\), set \(R = 1\). Otherwise, if \(\operatorname{Im}(z) \ge r\) or \(\operatorname{Re}(z) \ge 0 \land |z| \ge r\), set \(R = 2\). Otherwise, if \(|z| \ge 2r\), set \(R = 3\). Otherwise, the bound is infinite. If the bound is finite, we have

\[|\varepsilon_n(z)| \le 2 \alpha C_n \left|\frac{(a)_n (a-b+1)_n}{n! z^n} \right| \exp(2 \alpha \rho C_1 / |z|)\]

in terms of the following auxiliary quantities

\[\sigma = |(b-2a)/z|\]\[\begin{split}C_n = \begin{cases} 1 & \text{if } R = 1 \\ \chi(n) & \text{if } R = 2 \\ (\chi(n) + \rho \nu^2 n) \nu^n & \text{if } R = 3 \end{cases}\end{split}\]\[\nu = \left(\tfrac{1}{2} + \tfrac{1}{2}\sqrt{1-4\sigma^2}\right)^{-1/2} \le 1 + 2 \sigma^2\]\[\chi(n) = \sqrt{\pi} \Gamma(\tfrac{1}{2}n+1) / \Gamma(\tfrac{1}{2} n + \tfrac{1}{2})\]\[\begin{split}\sigma' = \begin{cases} \sigma & \text{if } R \ne 3 \\ \nu \sigma & \text{if } R = 3 \end{cases}\end{split}\]\[\alpha = (1 - \sigma')^{-1}\]\[\rho = \tfrac{1}{2} |2a^2-2ab+b| + \sigma' (1+ \tfrac{1}{4} \sigma') (1-\sigma')^{-2}\]
void acb_hypgeom_u_asymp(acb_t res, const acb_t a, const acb_t b, const acb_t z, long n, long prec)

Sets res to \(U^{*}(a,b,z)\) computed using n terms of the asymptotic series, with a rigorous bound for the error included in the output. We require \(n \ge 0\).

int acb_hypgeom_u_use_asymp(const acb_t z, long prec)

Heuristically determines whether the asymptotic series can be used to evaluate \(U(a,b,z)\) to prec accurate bits (assuming that a and b are small).

Confluent hypergeometric functions

void acb_hypgeom_u_1f1_series(acb_poly_t res, const acb_poly_t a, const acb_poly_t b, const acb_poly_t z, long len, long prec)

Computes \(U(a,b,z)\) as a power series truncated to length len, given \(a, b, z \in \mathbb{C}[[x]]\). If \(b[0] \in \mathbb{Z}\), it computes one extra derivative and removes the singularity (it is then assumed that \(b[1] \ne 0\)). As currently implemented, the output is indeterminate if \(b\) is nonexact and contains an integer.

void acb_hypgeom_u_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, long prec)

Computes \(U(a,b,z)\) as a sum of two convergent hypergeometric series. If \(b \in \mathbb{Z}\), it computes the limit value via acb_hypgeom_u_1f1_series(). As currently implemented, the output is indeterminate if \(b\) is nonexact and contains an integer.

void acb_hypgeom_u(acb_t res, const acb_t a, const acb_t b, const acb_t z, long prec)

Computes \(U(a,b,z)\) using an automatic algorithm choice. The function acb_hypgeom_u_asymp() is used if \(a\) or \(a-b+1\) is a nonpositive integer (in which case the asymptotic series terminates), or if z is sufficiently large. Otherwise acb_hypgeom_u_1f1() is used.

void acb_hypgeom_m_asymp(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, long prec)
void acb_hypgeom_m_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, long prec)
void acb_hypgeom_m(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, long prec)

Computes the confluent hypergeometric function \(M(a,b,z) = {}_1F_1(a,b,z)\), or \(\mathbf{M}(a,b,z) = \frac{1}{\Gamma(b)} {}_1F_1(a,b,z)\) if regularized is set.

The error function

void acb_hypgeom_erf_1f1a(acb_t res, const acb_t z, long prec)
void acb_hypgeom_erf_1f1b(acb_t res, const acb_t z, long prec)
void acb_hypgeom_erf_asymp(acb_t res, const acb_t z, long prec, long prec2)
void acb_hypgeom_erf(acb_t res, const acb_t z, long prec)

Computes the error function respectively using

\[\operatorname{erf}(z) = \frac{2z}{\sqrt{\pi}} {}_1F_1(\tfrac{1}{2}, \tfrac{3}{2}, -z^2)\]\[\operatorname{erf}(z) = \frac{2z e^{-z^2}}{\sqrt{\pi}} {}_1F_1(1, \tfrac{3}{2}, z^2)\]\[\operatorname{erf}(z) = \frac{z}{\sqrt{z^2}} \left(1 - \frac{e^{-z^2}}{\sqrt{\pi}} U(\tfrac{1}{2}, \tfrac{1}{2}, z^2)\right).\]

and an automatic algorithm choice. The asymp version takes a second precision to use for the U term.

void acb_hypgeom_erfc(acb_t res, const acb_t s, const acb_t z, long prec)

Computes the complementary error function \(\operatorname{erfc}(z) = 1 - \operatorname{erf}(z)\). This function avoids catastrophic cancellation for large positive z.

void acb_hypgeom_erfi(acb_t res, const acb_t s, const acb_t z, long prec)

Computes the imaginary error function \(\operatorname{erfi}(z) = -i\operatorname{erf}(iz)\). This is a trivial wrapper of acb_hypgeom_erf().

Bessel functions

void acb_hypgeom_bessel_j_asymp(acb_t res, const acb_t nu, const acb_t z, long prec)

Computes the Bessel function of the first kind via acb_hypgeom_u_asymp(). For all complex \(\nu, z\), we have

\[J_{\nu}(z) = \frac{z^{\nu}}{2^{\nu} e^{iz} \Gamma(\nu+1)} {}_1F_1(\nu+\tfrac{1}{2}, 2\nu+1, 2iz) = A_{+} B_{+} + A_{-} B_{-}\]

where

\[A_{\pm} = z^{\nu} (z^2)^{-\tfrac{1}{2}-\nu} (\mp i z)^{\tfrac{1}{2}+\nu} (2 \pi)^{-1/2} = (\pm iz)^{-1/2-\nu} z^{\nu} (2 \pi)^{-1/2}\]\[B_{\pm} = e^{\pm i z} U^{*}(\nu+\tfrac{1}{2}, 2\nu+1, \mp 2iz).\]

Nicer representations of the factors \(A_{\pm}\) can be given depending conditionally on the parameters. If \(\nu + \tfrac{1}{2} = n \in \mathbb{Z}\), we have \(A_{\pm} = (\pm i)^{n} (2 \pi z)^{-1/2}\). And if \(\operatorname{Re}(z) > 0\), we have \(A_{\pm} = \exp(\mp i [(2\nu+1)/4] \pi) (2 \pi z)^{-1/2}\).

void acb_hypgeom_bessel_j_0f1(acb_t res, const acb_t nu, const acb_t z, long prec)

Computes the Bessel function of the first kind from

\[J_{\nu}(z) = \frac{1}{\Gamma(\nu+1)} \left(\frac{z}{2}\right)^{\nu} {}_0F_1\left(\nu+1, -\frac{z^2}{4}\right).\]
void acb_hypgeom_bessel_j(acb_t res, const acb_t nu, const acb_t z, long prec)

Computes the Bessel function of the first kind \(J_{\nu}(z)\) using an automatic algorithm choice.

void acb_hypgeom_bessel_y(acb_t res, const acb_t nu, const acb_t z, long prec)

Computes the Bessel function of the second kind \(Y_{\nu}(z)\) from the formula

\[Y_{\nu}(z) = \frac{\cos(\nu \pi) J_{\nu}(z) - J_{-\nu}(z)}{\sin(\nu \pi)}\]

unless \(\nu = n\) is an integer in which case the limit value

\[Y_n(z) = -\frac{2}{\pi} \left( i^n K_n(iz) + \left[\log(iz)-\log(z)\right] J_n(z) \right)\]

is computed. As currently implemented, the output is indeterminate if \(\nu\) is nonexact and contains an integer.

void acb_hypgeom_bessel_i_asymp(acb_t res, const acb_t nu, const acb_t z, long prec)
void acb_hypgeom_bessel_i_0f1(acb_t res, const acb_t nu, const acb_t z, long prec)
void acb_hypgeom_bessel_i(acb_t res, const acb_t nu, const acb_t z, long prec)

Computes the modified Bessel function of the first kind \(I_{\nu}(z) = z^{\nu} (iz)^{-\nu} J_{\nu}(iz)\) respectively using asymptotic series (see acb_hypgeom_bessel_j_asymp()), the convergent series

\[I_{\nu}(z) = \frac{1}{\Gamma(\nu+1)} \left(\frac{z}{2}\right)^{\nu} {}_0F_1\left(\nu+1, \frac{z^2}{4}\right),\]

or an automatic algorithm choice.

void acb_hypgeom_bessel_k_asymp(acb_t res, const acb_t nu, const acb_t z, long prec)

Computes the modified Bessel function of the second kind via via acb_hypgeom_u_asymp(). For all \(\nu\) and all \(z \ne 0\), we have

\[K_{\nu}(z) = \left(\frac{2z}{\pi}\right)^{-1/2} e^{-z} U^{*}(\nu+\tfrac{1}{2}, 2\nu+1, 2z).\]
void acb_hypgeom_bessel_k_0f1_series(acb_poly_t res, const acb_poly_t nu, const acb_poly_t z, long len, long prec)

Computes the modified Bessel function of the second kind \(K_{\nu}(z)\) as a power series truncated to length len, given \(\nu, z \in \mathbb{C}[[x]]\). Uses the formula

\[K_{\nu}(z) = \frac{1}{2} \frac{\pi}{\sin(\pi \nu)} \left[ \left(\frac{z}{2}\right)^{-\nu} {}_0{\widetilde F}_1\left(1-\nu, \frac{z^2}{4}\right) - \left(\frac{z}{2}\right)^{\nu} {}_0{\widetilde F}_1\left(1+\nu, \frac{z^2}{4}\right) \right].\]

If \(\nu[0] \in \mathbb{Z}\), it computes one extra derivative and removes the singularity (it is then assumed that \(\nu[1] \ne 0\)). As currently implemented, the output is indeterminate if \(\nu[0]\) is nonexact and contains an integer.

void acb_hypgeom_bessel_k_0f1(acb_t res, const acb_t nu, const acb_t z, long prec)

Computes the modified Bessel function of the second kind from

\[K_{\nu}(z) = \frac{1}{2} \left[ \left(\frac{z}{2}\right)^{-\nu} \Gamma(\nu) {}_0F_1\left(1-\nu, \frac{z^2}{4}\right) - \left(\frac{z}{2}\right)^{\nu} \frac{\pi}{\nu \sin(\pi \nu) \Gamma(\nu)} {}_0F_1\left(\nu+1, \frac{z^2}{4}\right) \right]\]

if \(\nu \notin \mathbb{Z}\). If \(\nu \in \mathbb{Z}\), it computes the limit value via acb_hypgeom_bessel_k_0f1_series(). As currently implemented, the output is indeterminate if \(\nu\) is nonexact and contains an integer.

void acb_hypgeom_bessel_k(acb_t res, const acb_t nu, const acb_t z, long prec)

Computes the modified Bessel function of the second kind \(K_{\nu}(z)\) using an automatic algorithm choice.

Incomplete gamma functions

void acb_hypgeom_gamma_upper_asymp(acb_t res, const acb_t s, const acb_t z, int modified, long prec)
void acb_hypgeom_gamma_upper_1f1a(acb_t res, const acb_t s, const acb_t z, int modified, long prec)
void acb_hypgeom_gamma_upper_1f1b(acb_t res, const acb_t s, const acb_t z, int modified, long prec)
void acb_hypgeom_gamma_upper_singular(acb_t res, long s, const acb_t z, int modified, long prec)
void acb_hypgeom_gamma_upper(acb_t res, const acb_t s, const acb_t z, int modified, long prec)

Computes the upper incomplete gamma function respectively using

\[\Gamma(s,z) = e^{-z} U(1-s,1-s,z)\]\[\Gamma(s,z) = \Gamma(s) - \frac{z^s}{s} {}_1F_1(s, s+1, -z)\]\[\Gamma(s,z) = \Gamma(s) - \frac{z^s e^{-z}}{s} {}_1F_1(1, s+1, z)\]\[\Gamma(s,z) = \frac{(-1)^n}{n!} (\psi(n+1) - \log(z)) + \frac{(-1)^n}{(n+1)!} z \, {}_2F_2(1,1,2,2+n,-z) - z^{-n} \sum_{k=0}^{n-1} \frac{(-z)^k}{(k-n) k!}, \quad n = -s \in \mathbb{Z}_{\ge 0}\]

and an automatic algorithm choice. The automatic version also handles other special input such as \(z = 0\) and \(s = 1, 2, 3\). The singular version evaluates the finite sum directly and therefore assumes that s is not too large. If modified is set, computes the exponential integral \(z^{-s} \Gamma(s,z) = E_{1-s}(z)\) instead.

Exponential and trigonometric integrals

The branch cut conventions of the following functions match Mathematica.

void acb_hypgeom_expint(acb_t res, const acb_t s, const acb_t z, long prec)

Computes the generalized exponential integral \(E_s(z)\). This is a trivial wrapper of acb_hypgeom_gamma_upper().

void acb_hypgeom_ei_asymp(acb_t res, const acb_t z, long prec)
void acb_hypgeom_ei_2f2(acb_t res, const acb_t z, long prec)
void acb_hypgeom_ei(acb_t res, const acb_t z, long prec)

Computes the exponential integral \(\operatorname{Ei}(z)\), respectively using

\[\operatorname{Ei}(z) = -e^z U(1,1,-z) - \log(-z) + \frac{1}{2} \left(\log(z) - \log\left(\frac{1}{z}\right) \right)\]\[\operatorname{Ei}(z) = z {}_2F_2(1, 1; 2, 2; z) + \gamma + \frac{1}{2} \left(\log(z) - \log\left(\frac{1}{z}\right) \right)\]

and an automatic algorithm choice.

void acb_hypgeom_si_asymp(acb_t res, const acb_t z, long prec)
void acb_hypgeom_si_1f2(acb_t res, const acb_t z, long prec)
void acb_hypgeom_si(acb_t res, const acb_t z, long prec)

Computes the sine integral \(\operatorname{Si}(z)\), respectively using

\[\operatorname{Si}(z) = \frac{i}{2} \left[ e^{iz} U(1,1,-iz) - e^{-iz} U(1,1,iz) + \log(-iz) - \log(iz) \right]\]\[\operatorname{Si}(z) = z {}_1F_2(\tfrac{1}{2}; \tfrac{3}{2}, \tfrac{3}{2}; -\tfrac{z^2}{4})\]

and an automatic algorithm choice.

void acb_hypgeom_ci_asymp(acb_t res, const acb_t z, long prec)
void acb_hypgeom_ci_2f3(acb_t res, const acb_t z, long prec)
void acb_hypgeom_ci(acb_t res, const acb_t z, long prec)

Computes the cosine integral \(\operatorname{Ci}(z)\), respectively using

\[\operatorname{Ci}(z) = \log(z) - \frac{1}{2} \left[ e^{iz} U(1,1,-iz) + e^{-iz} U(1,1,iz) + \log(-iz) + \log(iz) \right]\]\[\operatorname{Ci}(z) = -\tfrac{z^2}{4} {}_2F_3(1, 1; 2, 2, \tfrac{3}{2}; -\tfrac{z^2}{4}) + \log(z) + \gamma\]

and an automatic algorithm choice.

void acb_hypgeom_shi(acb_t res, const acb_t z, long prec)

Computes the hyperbolic sine integral \(\operatorname{Shi}(z) = -i \operatorname{Si}(iz)\). This is a trivial wrapper of acb_hypgeom_si().

void acb_hypgeom_chi_asymp(acb_t res, const acb_t z, long prec)
void acb_hypgeom_chi_2f3(acb_t res, const acb_t z, long prec)
void acb_hypgeom_chi(acb_t res, const acb_t z, long prec)

Computes the hyperbolic cosine integral \(\operatorname{Chi}(z)\), respectively using

\[\operatorname{Chi}(z) = -\frac{1}{2} \left[ e^{z} U(1,1,-z) + e^{-z} U(1,1,z) + \log(-z) - \log(z) \right]\]\[\operatorname{Chi}(z) = \tfrac{z^2}{4} {}_2F_3(1, 1; 2, 2, \tfrac{3}{2}; \tfrac{z^2}{4}) + \log(z) + \gamma\]

and an automatic algorithm choice.

void acb_hypgeom_li(acb_t res, const acb_t z, int offset, long prec)

If offset is zero, computes the logarithmic integral \(\operatorname{li}(z) = \operatorname{Ei}(\log(z))\).

If offset is nonzero, computes the offset logarithmic integral \(\operatorname{Li}(z) = \operatorname{li}(z) - \operatorname{li}(2)\).