Numerical integration in arbitrary-precision ball arithmetic

We present an implementation of arbitrary-precision numerical integration with rigorous error bounds in the Arb library. Rapid convergence is ensured for piecewise complex analytic integrals by use of the Petras algorithm, which combines adaptive bisection with adaptive Gaussian quadrature where error bounds are determined via complex magnitudes without evaluating derivatives. The code is general, easy to use, and efficient, often outperforming existing non-rigorous software.


Introduction
Many users can attest that there is a non-negligible chance of getting an incorrect answer when asking a numerical package or computer algebra system for an approximation of a definite integral b a f (x)dx, as rapid variation, narrow peaks, non-smooth points, cancellation or ill-conditioned numerical evaluation of f are prone to break widely used heuristic numerical integration methods.
One remedy is to compute rigorous error bounds using interval arithmetic. However, little work has been done to date on efficient arbitrary-precision implementations. Here, we present a new implementation of rigorous numerical integration in Arb, a C library for ball arithmetic 1 [4]. The integration code is easy to use directly in C, or can be wrapped from high-level languages. For example, an interface in Sage [11] exists (thanks to Marc Mezzarobba and Vincent Delecroix), which we demonstrate by computing We obtain nearly 100 digits with a rigorous error bound in 0.04 seconds (0.02 s when using C directly). This relatively difficult test integral (f (x) changes sign 950 times) was introduced by Rump [9] who observed that the quad function in Matlab took over a second only to return the erroneous 0.2511 (Rump's interval package Intlab computes 7 digits in about one second; see also [7]).

Algorithm and implementation
We consider integration of a function f : C → C on a segment [a, b], a, b ∈ C. We represent real numbers as mid-rad intervals (balls) [m ± r] and complex numbers as rectangles [m 1 ±r 1 ]+[m 2 ±r 2 ]i (which we also refer to as balls with slight abuse of terminology). True complex balls B(m 1 + m 2 i, r) would sometimes provide slightly better bounds, but rectangles are usually more convenient.
The user supplies the integrand f as a pointer to a C function func implementing its evaluation (we refer to the documentation for the detailed API). In effect, func gets called with the argument z and an extra flag d. If d = 0, func is to evaluate f (z) without any assumptions about regularity. If d = 1, func is to evaluate f (z) and also check that f is analytic on z, returning a non-finite ball (e.g. NaN) otherwise. For meromorphic f , the user can ignore d since f (z) automatically blows up at poles, but d needs to be handled for functions with branch cuts like √ z and log(z) (here by checking whether z overlaps (−∞, 0]). We use the Petras algorithm [8], which combines bisection with Gaussian quadrature of variable degree n. Error bounds for Gaussian quadrature use complex magnitudes. If f is analytic with |f | ≤ M on an ellipse E with foci ±1 and semiaxes X, Y, then | hits a singularity of f , but also improves convergence as n → ∞. Of course, the computed bound for M will not just depend on the function f but also on the stability of its evaluation in ball arithmetic if E is large.
Degree adaptivity ensures near-optimal complexity (O(p) evaluations of f ) for analytic f at high precision p, while space adaptivity (bisection) helps if there are singularities near [a, b] or if the ball enclosures are not optimal. For piecewise analytic f with discontinuities on [a, b] the complexity is typically O(p 2 ), i.e. a bit worse but still polynomial in p. Degree or space adaptivity used alone would give 2 O(p) complexity or fail to converge for common types of integrals.
Our version of the integration algorithm can be described as follows: . If I meets the tolerance goal, if α, β overlap, or if evaluation limits have been exceeded, set S ← S + I and go to 1. 3. Try to find an ellipse E with foci (α, β) and an n ≤ n max such that f is analytic on E (evaluating f (E) with d = 1) and the error bound for n-point Gaussian quadrature determined via |f (E)| meets the tolerance goal. If successful, compute this integral J, set S ← S + J and go to 1. 4. Interval bisection: let m = α+β 2 and extend Q with (α, m), (m, β). Compared to Petras [8], there are minor differences. Our ρ is not fixed; we try several sizes of E in step 3 to reduce n. The handling of tolerances is slightly different. We also compute quadrature nodes (w k , x k ) at runtime, without using pre-made tables. A key point is that generating nodes for high-precision Gaussian quadrature used to be considered too costly [1], but the recent work [6] solves this problem. 2 With default settings, computing nodes takes a few milliseconds for 100-digit precision and a few seconds for 1000 digits. 3 Nodes are automatically cached, so this cost is amortized for repeated integrations at the same or lower precision (possible n are restricted to a sparse sequence ≈ 2 k/2 to avoid computing nodes for many nearby n). As an optional tuning parameter, the user can change the allowed range of n which defaults to n max = 0.5p + 60.

Tolerances and evaluation limits
Besides the working precision p, the user specifies absolute and relative tolerances ε abs and ε rel . In effect, the algorithm attempts to achieve an error of max(ε abs , V ε rel ) where V is the magnitude of the integral. Reasonable values (used as defaults by the Sage wrapper) are ε abs = ε rel = 2 −p . Other values can be useful, e.g. if low accuracy is sufficient but a higher p must be used for numerical reasons. One might also set ε abs = 0 to use relative tolerance only, though for efficiency, it is better to supply ε abs ≈ V ε rel if an estimate for V is known when V ≈ 1. This Sage code shows computation of In reality, ε abs and ε rel are only guidelines and the algorithm does not strictly achieve the goal max(ε abs , V ε rel ). Indeed, due to the fixed working precision and possibly inexact parameters, the goal cannot generally be achieved. It is implied that the user will work with some guard bits and if needed adjust (p, ε abs , ε rel ) based on the reliable a posteriori information in the output ball radius.
Use of ε rel further depends circularly on V (V is essentially what we are trying to compute!), so the algorithm must guess V. A too large guess means loss of accuracy and a too small guess means unnecessary work. Our approach is to start with the tolerance ε abs and continuously update ε abs ← max(ε abs , I a ε rel ), ε abs ← max(ε abs , J a ε rel ) where |I| = [I a , I b ] and |J| = [J a , J b ] are intervals computed in steps 2 and 3; I a and J a will then be lower bounds for V (we err on the side of preserving accuracy), modulo global cancellation in the integral. As noted above, the user should exploit knowledge about V if possible since I a and J a may be pessimistic. More clever globally adaptive strategies are possible, but we settled for this simple approach in the present version.
To abort gracefully when convergence is too slow, evaluation limits include a bound on the number of calls to f (default 1000p + p 2 ) and a bound on the size N of the work queue Q (default 2p). By default, Q acts as a stack and step 4 puts the new subinterval with the larger error at the top; optionally, Q can be switched to a global priority queue, which may improve results if convergence is so slow that evaluation limits are exceeded. This has the downside of sometimes requiring N nearly as large as the number of calls to f (e.g. N ∼ p 2 ), whereas we always have N p with the stack. A more clever algorithm might use a top-level priority queue down to some depth before switching to a stack locally.

Benchmarks
We test various integrals with precision between about 10 and 1000 digits. Timings were obtained on an Intel Core i5-4300U CPU. We compare Arb to the heuristic arbitrary-precision integration routines intnum in Pari/GP [10] and quad in mpmath [5]. Both use double exponential quadrature without adaptive subdivision, although quad is degree-adaptive. Further comparisons with other numerical and interval packages (as well as alternative methods in Pari/GP and mpmath 4 ) would be useful, but out of scope for this brief overview.
We do not show the outputs, but note that in all cases, Arb computes correct balls with radius a small multiple of 2 −p . On some test cases, mpmath with default settings silently returns an inaccurate answer due to exceeding its limit on the quadrature degree, but it provides an optional mechanism to catch this. We increased the degree limit to let mpmath run to full accuracy in all cases, and have written (!) after a timing where the default is insufficient. Pari/GP is not adaptive and silently returns inaccurate answers without providing a catch mechanism or a way to increase the degree. It does provide an option to split the interval non-adaptively into 2 t parts, but it is up to the user to find a correct t. We have done so where necessary, which is also marked in the timings.  3.1 Integrals without singularities on the path Table 1 shows examples with smooth f on [a, b]. For meromorphic f , the number of subintervals largely depends on the location of the poles and does not change with p. The "spike integral" I 1 (Figure 1) is a well known pathological example [2,3]; all ordinary numerical integrators we have tested (Mathematica, GSL, SciPy, etc.) give inaccurate results with default settings. This integrand has poles near the real axis, forcing many local bisections. It is a piece of cake for the Petras algorithm, but Pari/GP and mpmath converge slowly unless the user manually splits the path at the peaks. I 2 could be sped up 40% in Arb by using cos 2 (x) = 1 2 (1+cos(2x)) for wide x to bound the denominator more tightly. For entire functions (I 4 , I 5 , I 6 ), the efficiency improves with larger p since arbitrarily large bounding ellipses can be used. I 5 is Rump's example again, and I 6 (whose graph has two sharp "bends") was provided by Silviu-Ioan Filip.
The code is seen to work well with special functions. In I 3 we integrate the Lambert W function, where we need to check for the branch cut on (−∞, −1/e] in the evaluation. I 7 also illustrates integration on a complex path. Overall, Arb is faster than Pari/GP and mpmath, despite the fact that rigorous error bounds create extra work. The speedup is in part explained by faster arithmetic and transcendental functions in Arb and lower overhead due to using C, as well as the advantage of Gaussian quadrature over the double exponential method for smooth integrands. However, if these differences are accounted for, we can still conclude that the Petras algorithm in ball arithmetic holds up extremely well for high-precision integration, on top of giving rigorous bounds. Table 2. Improper integrals and integrals with endpoint singularities. For integration with Arb, all improper integrals (i.e. excluding E0) have been truncated manually at a lower bound ε or upper bound N , chosen so that the omitted part is smaller than 2 −p .

Endpoint singularities and infinite intervals
The methods in Pari/GP and mpmath are designed to support typical integrals with infinite intervals or endpoint singularities, which often arise in applications. Arb requires finite a, b and a bounded f to return a finite result, but the user may provide a manual truncation (say to work around this restriction. Tail bounds must then be added based on symbolic knowledge about f . This is not ideal in terms of usability or efficiency, but since the Petras algorithm works well even with an endpoint very close to a singularity (or ∞), evaluating improper integrals to high precision in this way is at least feasible. 5 In Table 2, E 0 , E 1 and E 2 have algebraic or logarithmic singularities or decay, with E 1 requiring N ≈ 2 p and E 2 requiring ε ≈ 2 −p (no truncation is needed for E 0 , as f is bounded at the algebraic branch point singularity x = 1). Here Arb needs O(p) subintervals and O(p 2 ) evaluations, while the double exponential algorithm in Pari/GP and mpmath only needs roughly O(p) evaluations and therefore scales better. 6 For integrals with exponential decay (E 3 , E 4 and E 5 ), a cutoff of N ∼ p is sufficient, and here Arb retains excellent performance.
In a future extension of this work, some reasonable class of improper integrals could be supported more efficiently and conveniently (e.g. with the user providing a symbolic bound like |f (x)| < Cx α exp(−βx γ )). 5 An exception is when f has an essential singularity inducing oscillation combined with slow decay. Oscillation with exponential decay is not a problem (as in E4, E5), but integrals like work, so we can only hope for 5-10 digits without specialized oscillatory algorithms. 6 As a means to improve performance, we note the standard trick of manually changing variables to turn algebraic growth or decay into exponential decay. Indeed, x → sinh(x) gives E1 = E3. Similarly x → tanh(x) and x → e −x can be used in E0, E2. Table 3. Integrals with point discontinuities in f or f . Here p(x) = x 4 + 10x 3 + 19x 2 − 6x − 6 in D0, and u(x) = (x − x − 1 2 ), v = max(sin(x), cos(x)) in D3. For D3, the function evaluation limit had to be increased for convergence at p = 3333.

Piecewise and discontinuous functions
Piecewise real analytic functions can be integrated efficiently using piecewise complex analytic extensions. For example, |x| on R extends to the function √ z 2 of z = x + yi, which equals z in the right plane and −z in the left plane with a branch cut on Re(z) = 0. 7 We provide as library methods such extensions of sgn(x), |x|, x , x , max(x, y), min(x, y), with builtin branch cut detection. Table 3 shows integrals with mid-interval jumps or kinks, including one complex integral crossing a branch cut discontinuity (D 2 ). The example D 0 , where p(x) changes sign once on [0, 1], is due to Helfgott (see comments in [7]).
We see that a mid-interval singularity leads to use of O(p) subintervals and O(p 2 ) evaluations to isolate the problematic point by bisection. With k such points (D 1 and D 3 ), the cost simply increases by another factor k, and the user may have to raise the evaluation limits accordingly to let the algorithm complete (which we did for D 3 ). In contrast, Pari/GP and mpmath cope poorly with midinterval singularities and cannot achieve high accuracy on these examples unless the user manually splits the interval precisely at the problematic points.

Complex analysis
We conclude by illustrating integration as a tool for complex analysis. First, we consider computing derivatives via the Cauchy integral formula. Denote by ℘(z; τ ) = ∞ n=−2 a n (τ )z n the Weierstrass elliptic function for the lattice (1, τ ). We fix τ = i (placing the poles of ℘ at the Gaussian integers) and compute the Laurent coefficients a n = 1 2πi γ z −n−1 ℘(z)dz by integrating along the square connecting ±0.5 ± 0.5i. We ignore symmetry and compute all four segments. With p = 333, some results are (note that a −1 = a 100 = 0 and all a n are real): We lose about n bits of precision to cancellation due to the integrand magnitude growing with n. Apart from this, the difficulty increases quite slowly with n: a −2 takes 0.67 s while a 98 and a 100 take 0.85 s at this precision.
As We obtain balls that provably determine N (T ), and the method scales reasonably well. Unfortunately, the evaluation of ζ(s) in Arb is currently not well tuned for all s, which makes large T slower than necessary and can make this computation extremely slow with slightly different settings. In general, for complicated integrals, the user may need to customize the integrand evaluation to handle wide balls or large parameters optimally for a given path and precision.