Error functions of a real ball

arb_hypgeom_erf and arb_hypgeom_erfc compute the error function $\operatorname{erf}(x) = \frac{2}{\sqrt\pi} \int_0^x e^{-t^2} dt$ and its complement $\operatorname{erfc}(x) = 1 - \operatorname{erf}(x)$ for a real ball $x$; arb_hypgeom_erfinv and arb_hypgeom_erfcinv compute the inverse functions. Unlike most real functions of this module they do not call the complex implementation: for real $x$ one of the two convergent series has only positive terms, the asymptotic expansion has a simple error bound, and at very high precision the incomplete gamma function $\Gamma(1/2, x^2)$ can be computed with the bit-burst algorithm. acb_hypgeom_erf calls arb_hypgeom_erf for real input. The complex algorithms and the incomplete gamma functions in general are described in Error functions, incomplete gamma functions and exponential integrals.

The methods#

arb_hypgeom_erf chooses among the following.

• Tiny or huge $|x|$. For $|x| < 2^{-\mathit{prec}/16}$ the series below needs few terms. If $x^2/\log 2$ exceeds the precision (or $|x| > 2^{60}$), the function returns $\operatorname{sgn}(x) \pm \varepsilon$ with $\varepsilon = e^{-x^2}/(|x|\sqrt\pi)$, which bounds $|\operatorname{erfc}(|x|)|$.

• Convergent series. The default method is

  $$\operatorname{erf}(x) = \frac{2x}{\sqrt\pi} e^{-x^2} \sum_{k=0}^\infty \frac{(x^2)^k}{(3/2)_k},$$

  in which all terms are positive, so that no precision is lost to cancellation. The number of terms $N$ is estimated with the Lambert $W$ function, the tail is bounded by the tail of the exponential series ($(x^2)^k/(3/2)_k \le (x^2)^k/k!$), and the sum is computed by rectangular splitting with _arb_hypgeom_gamma_lower_sum_rs_1, so that most multiplications are by small integers.

• Asymptotic expansion. For $|x| > 2$, the expansion

  $$\operatorname{erfc}(x) = \frac{e^{-x^2}}{x\sqrt\pi} \left(\sum_{k=0}^{N-1} \frac{(1/2)_k}{(-x^2)^k} + R_N\right), \qquad |R_N| \le \frac{(1/2)_N}{x^{2N}} \le \frac{N!}{x^{2N}},$$

  is used when some $N$ makes the bound smaller than $2^{-\mathit{prec}}$; the sum is computed at the reduced precision that the factor $e^{-x^2}$ allows.

• Bit-burst algorithm. Above 30000 bits and for $150 e^{-0.004\sqrt{\mathit{prec}}} < |x| < 0.6\sqrt{\mathit{prec}}$, the function computes $\operatorname{erfc}(|x|) = \Gamma(1/2, x^2)/\sqrt\pi$ with arb_hypgeom_erf_bb: $x^2$ is approached by numbers $z_0, z_1, \ldots$ with 64, 128, 256, ... bits, the value at $z_0$ comes from the asymptotic or the convergent series of the incomplete gamma function, and each step from $z_j$ to $z_{j+1}$ sums a Taylor series by binary splitting. The cost is $O(M(n) \log^2 n)$ at $n$ bits, asymptotically less than the $O(n)$ multiplications by small integers and $O(\sqrt n)$ full multiplications of rectangular splitting. See The bit-burst algorithm.

For inexact $x$, the series are evaluated at the midpoint and the propagated error is bounded by $\frac{2}{\sqrt\pi} e^{-m^2} r$, where $m$ is the lower bound of $|x|$ and $r$ the radius (and by 2); the precision is reduced to the accuracy of $x$ plus the number of bits that $e^{-x^2}$ gains.

arb_hypgeom_erfc computes $1 - \operatorname{erf}(x)$ for $x \le 1$. For larger $x$, where the subtraction would cancel, it uses the asymptotic expansion whenever it reaches the target accuracy (with a single term for $x > 2^{\mathit{prec}/2+10}$), the bit-burst algorithm above 30000 bits, and otherwise the convergent series for $\operatorname{erf}$ at a working precision increased by $(x^2 + \log x)/\log 2$ bits.

Inverse functions#

arb_hypgeom_erfinv solves $\operatorname{erf}(y) = x$ by Newton's method with $1/\operatorname{erf}'(y) = \frac{\sqrt\pi}{2} e^{y^2}$. The starting value comes from a rational approximation in double precision [Bla1976], from the first terms of the asymptotic expansion of $\operatorname{erf}^{-1}$ near $\pm 1$ when $1 - |x| < 2^{-52}$, and from the same expansion in ball arithmetic when $1 - |x| < 2^{-1000}$. Each iteration doubles the precision; when an iterate is an interval Newton step contained in the previous enclosure, the existence of the root in it is proved, and the loop stops once a proved enclosure has the target accuracy. Near $|x| = 1$ the iteration works with $1 - |x|$ and $\operatorname{erfc}$, which avoids cancellation. arb_hypgeom_erfcinv passes $x$ itself in this form when $x \le 0.01$, so that arguments such as $10^{-1000}$ lose no accuracy.

Example#

At 64 bits, $\operatorname{erf}(30)$ is $1$ up to a tiny error, and only arb_hypgeom_erfc resolves the difference; arb_hypgeom_erfcinv inverts it for a tiny argument:

#include <flint/arb.h>
#include <flint/arb_hypgeom.h>

int main(void)
{
    arb_t x, y, z;

    arb_init(x);
    arb_init(y);
    arb_init(z);

    arb_set_ui(x, 30);
    arb_hypgeom_erf(y, x, 64);
    flint_printf("erf(30)          = "); arb_printn(y, 20, 0); flint_printf("\n");
    arb_hypgeom_erfc(y, x, 64);
    flint_printf("erfc(30)         = "); arb_printn(y, 20, 0); flint_printf("\n");
    /* 1 - erf(30) cannot resolve the tiny difference at 64 bits */
    arb_hypgeom_erf(y, x, 64);
    arb_sub_ui(y, y, 1, 64);
    arb_neg(y, y);
    flint_printf("1 - erf(30)      = "); arb_printn(y, 20, 0); flint_printf("\n");

    /* inverse functions; erfcinv works directly with tiny arguments */
    arb_set_str(x, "0.5", 64);
    arb_hypgeom_erfinv(y, x, 64);
    flint_printf("erfinv(0.5)      = "); arb_printn(y, 20, 0); flint_printf("\n");
    arb_set_str(x, "1e-1000", 64);
    arb_hypgeom_erfcinv(y, x, 64);
    flint_printf("erfcinv(1e-1000) = "); arb_printn(y, 20, 0); flint_printf("\n");
    arb_hypgeom_erfc(z, y, 64);
    flint_printf("erfc of that     = "); arb_printn(z, 20, 0); flint_printf("\n");

    arb_clear(x);
    arb_clear(y);
    arb_clear(z);
    flint_cleanup();
    return 0;
}

erf(30)          = [1.0000000000000000000 +/- 2.57e-393]
erfc(30)         = [2.5646562037561116000e-393 +/- 3.65e-413]
1 - erf(30)      = [+/- 2.57e-393]
erfinv(0.5)      = [0.4769362762044698734 +/- 6.83e-20]
erfcinv(1e-1000) = [47.93894640426451069 +/- 3.99e-18]
erfc of that     = [1.000000000000000e-1000 +/- 3.75e-1016]

The enclosure of $\operatorname{erf}(30)$ is the bound $1 \pm e^{-900}/(30\sqrt\pi)$, which is correct but says nothing about $\operatorname{erfc}(30)$ itself; computing $\operatorname{erfc}$ directly gives 64 significant bits. The string "1e-1000" is not exactly representable, so the input to arb_hypgeom_erfcinv is a ball of relative radius $2^{-64}$, and the last line recovers it within that radius.