Trapezoidal integration with error bounds using Arb

trapezoid.c

/*
Naive implementation of the trapezoidal rule for integration.
(A non-naive implementation would provide adaptivity, etc.)

WARNING: I have not checked carefully that the code is correct.
Let me know if there is a bug.
*/

#include "arb.h"
#include "arb_poly.h"
#include "arb_hypgeom.h"

/* Sets res to f(g) where g is an element of R[[t]], truncating the result
   to O(t^len). Here f(x) = erf(1+erf(x^2)). */
void
integrand(arb_poly_t res, const arb_poly_t g, slong len, slong prec)
{
    arb_poly_mullow(res, g, g, len, prec);
    arb_hypgeom_erf_series(res, res, len, prec);
    arb_poly_add_si(res, res, 1, prec);
    arb_hypgeom_erf_series(res, res, len, prec);
}

/* Sets res to f^{(k)}(x), computed by evaluating the power series
   and reading off the last coefficient. This could be overloaded
   for speed when a single derivative can be computed faster.
   poly1 and poly2 are just used as temporaries here. */
void
fk(arb_t res, const arb_t x, slong k, slong prec,
    arb_poly_t poly1, arb_poly_t poly2)
{
    /* poly1 = x + t   (t is the dummy power series variable) */
    arb_poly_set_arb(poly1, x);
    arb_poly_set_coeff_si(poly1, 1, 1);

    integrand(poly2, poly1, k + 1, prec);

    /* read off the coefficient */
    arb_poly_get_coeff_arb(res, poly2, k);

    /* don't forget to multiply by k! */
    if (k >= 2)
    {
        arb_t t;
        arb_init(t);
        arb_fac_ui(t, k, prec);
        arb_mul(res, res, t, prec);
        arb_clear(t);
    }
}

/* Sets res to the integral of f on [a,b], evaluated using
   the (N+1)-point trapezoid rule, with Nbound interval evaluations of f''
   to bound the error. */
void
trapezoid(arb_t res, const arb_t a, const arb_t b, slong N, slong Nbound, slong prec)
{
    arb_t t, u, h, bound;
    arb_poly_t poly1, poly2;
    slong k;

    arb_init(t);
    arb_init(u);
    arb_init(h);
    arb_init(bound);
    arb_poly_init(poly1);
    arb_poly_init(poly2);

    /* First, we compute the error bound. */
    arb_zero(bound);

    /* h = (b - a) / Nbound */
    arb_sub(h, b, a, prec);
    arb_div_ui(h, h, Nbound, prec);

    /* Bound |f''| on subintervals. */
    for (k = 0; k < Nbound; k++)
    {
        arb_set(t, a);
        arb_addmul_ui(t, h, k, prec);
        arb_add(u, t, h, prec);
        arb_union(t, t, u, prec);

        fk(t, t, 2, prec, poly1, poly2);

        arb_abs(t, t);
        arb_max(bound, bound, t, prec);
    }

    /* Trapezoid rule error bound = (b-a)^3 / (12 N^2) * bound */
    arb_sub(t, b, a, prec);
    arb_pow_ui(t, t, 3, prec);
    /* Warning: 12*N*N could overflow in the input for huge N --
       should split this computation to avoid. */
    arb_div_ui(t, t, 12 * N * N, prec);
    arb_mul(bound, bound, t, prec);

    /* Now compute the trapezoid sum. */
    arb_zero(res);

    /* h = (b - a) / N */
    arb_sub(h, b, a, prec);
    arb_div_ui(h, h, N, prec);

    for (k = 0; k <= N; k++)
    {
        arb_set(t, a);
        arb_addmul_ui(t, h, k, prec);

        fk(t, t, 0, prec, poly1, poly2);

        if (k == 0 || k == N)  /* first and last points have half weight */
            arb_mul_2exp_si(t, t, -1);
        arb_add(res, res, t, prec);
    }

    arb_mul(res, res, h, prec);

    /* add bound to the radius of res */
    arb_add_error(res, bound);

    arb_clear(t);
    arb_clear(u);
    arb_clear(h);
    arb_clear(bound);
    arb_poly_clear(poly1);
    arb_poly_clear(poly2);
}

int main()
{
    arb_t a, b, res;
    slong N, Nbound, prec;

    arb_init(a);
    arb_init(b);
    arb_init(res);

    prec = 53;
    N = 1000;
    Nbound = 100;

    /* integrate f(x) on [sqrt(2), pi] */
    arb_sqrt_ui(a, 2, prec);
    arb_const_pi(b, prec);
    trapezoid(res, a, b, N, Nbound, prec);

    arb_printn(res, 30, 0);
    printf("\n");

    arb_clear(a);
    arb_clear(b);
    arb_clear(res);

    flint_cleanup();
}