Integrate over a triangle in the 2D normal distribution

Question

I'd like to evaluate following expression efficiently (numerically).

$$g_a(x) := \int_0^x e^{-t^2} \int_0^{at} e^{-s^2} dsdt$$

If I want a given fixed accuracy, and evaluate both integrals using e.g. a summed quadrature formula for each integral, the cost is about $O(ax^2)$. Is there a better way than this naive one?

Answer 1

As the commenter noted, implicit in the question is whether it is better to do 2D quadrature than evaluate an error function at each sample point of a 1D quadrature; i.e., is there structure in the definition of the error function that allows speedup relative to the naive quadrature implementation? So I wrote a little code to find out: Using boost::math::erf and std::exp, I determined that it is faster to evaluate erf than std::exp:

Run on (8 X 3700 MHz CPU s)
------------------------------------------------------------
Benchmark                     Time           CPU Iterations
------------------------------------------------------------
BM_Erf<float>                 5 ns          5 ns  110919662
BM_Erf<double>                9 ns          9 ns   73097986
BM_Erf<long double>           8 ns          8 ns   73967212
BM_Bell<float>               11 ns         11 ns   59276009
BM_Bell<double>              12 ns         12 ns   52942894
BM_Bell<long double>        197 ns        197 ns   13041274

So it's faster to evaluate erf than exp (though keep in mind that boost is not bound by the 1ULP accuracy goal of the standard library math functions). Next you just need a good quadrature rule, say, tanh-sinh for small $x$ and trapezoidal for large. These will converge exponentially in about 30 function evaluations (double precision), so your total compute time should be about $(12+9)*30 \mathrm{ns} \approx 0.5\mu s$.

That's as fast as I can think to get it.

Code to compute the speed:

#include <cmath>
#include <ostream>
#include <random>
#include <boost/math/special_functions/erf.hpp>
#include <benchmark/benchmark.h>

template<typename Real>
static void BM_Erf(benchmark::State& state) {
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<Real> dis(0, 1000);
    auto s = dis(gen);
    Real y;
    while (state.KeepRunning()) {
        benchmark::DoNotOptimize(y = boost::math::erf(s));
    }
    std::ostream cnull(nullptr);
    cnull << y;
}

template<typename Real>
static void BM_Bell(benchmark::State& state) {
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<Real> dis(0, 1000);
    auto s = dis(gen);
    Real y;
    while (state.KeepRunning()) {
        benchmark::DoNotOptimize(y = std::exp(-s*s));
    }
    std::ostream cnull(nullptr);
    cnull << y;
}


BENCHMARK_TEMPLATE(BM_Erf, float);
BENCHMARK_TEMPLATE(BM_Erf, double);
BENCHMARK_TEMPLATE(BM_Erf, long double);

BENCHMARK_TEMPLATE(BM_Bell, float);
BENCHMARK_TEMPLATE(BM_Bell, double);
BENCHMARK_TEMPLATE(BM_Bell, long double);


BENCHMARK_MAIN();