#include <blitz/vector.h>

using namespace blitz;

int main()
{
    // Integrate the expression
    //            x
    //   2        /
    // -------    | exp(-t^2) dt
    // sqrt(Pi)   /
    //           0
    // to estimate the error function erf(x) on the interval [0,1].
    //
    // Three methods are compared:
    //  1. Naive summation
    //  2. Two-point trapezoid
    //  3. Three-point Simpson's

    const int numSamples = 1024;       // Number of samples
    double dt = 1.0 / numSamples;      // Distance between samples
    Range R(0, numSamples - 1);        // Index set for sample vectors
    Range displayRange(0, 900, 100);   // Indices to output to screen

    // For comparison purposes, use the built-in erf(x) function.
    Vector<double> z = erf(R * dt);
    cout << "  Built-in: " << z(displayRange+1) << endl;


    // Naive summation
    Vector<double> z1 = accumulate(M_2_SQRTPI * exp(-sqr(R * dt)) * dt);
    cout << "     Naive: " << z1(displayRange+1) << endl;

    // Calculate the rms error
    double error1 = sqrt(mean(sqr(z1 - z)));
    cout << "RMS Error: " << error1 << endl;


    // For the following rules, it's easier to work with a
    // sampled integrand.
    Vector<double> a = M_2_SQRTPI * exp(-sqr(R * dt));


    // Two-point trapezoid
    Range J(1, numSamples-1);
    Vector<double> z2 = accumulate(dt / 2.0 * (a(J) + a(J-1)));

    cout << " Trapezoid: " << z2(displayRange) << endl;
    cout << "RMS Error: " << sqrt(mean(sqr(z2-z(J)))) << endl;


    // Three-point Simpson's (parabolic)
    Range I(1, numSamples-2);
    Vector<double> z3 = accumulate(dt / 6.0 * (a(I-1) + 4 * a(I) + a(I+1)));

    cout << " Parabolic: " << z3(displayRange) << endl;
    cout << "RMS Error: " << sqrt(mean(sqr(z3 - z(I)))) << endl;


    return 0;
}