#!/usr/bin/ruby

#
## http://rosettacode.org/wiki/Numerical_integration/Gauss-Legendre_Quadrature
#

func legendre_pair((1), x) { (x, 1) }
func legendre_pair( n,  x) {
    var (m1, m2) = legendre_pair(n - 1, x);
    var u = (1 - 1/n);
    ((1 + u)*x*m1 - u*m2, m1);
}

func legendre((0), _) { 1 }
func legendre( n,  x) { [legendre_pair(n, x)][0] }

func legendre_prime({ .is_zero }, _) { 0 }
func legendre_prime({ .is_one  }, _) { 1 }

func legendre_prime(n, x) {
    var (m0, m1) = legendre_pair(n, x);
    (m1 - x*m0) * n / (1 - x**2);
}

func approximate_legendre_root(n, k) {
    # Approximation due to Francesco Tricomi
    var t = ((4*k - 1) / (4*n + 2));
    (1 - ((n - 1)/(8 * n**3))) * (Num.pi * t -> cos);
}

func newton_raphson(f, f_prime, r, eps = 2e-16) {
    while (var dr = (-f(r) / f_prime(r)) -> abs >= eps) {
        r += dr.float!;
    }
    return r;
}

func legendre_root(n, k) {
    newton_raphson(legendre.method(:call, n), legendre_prime.method(:call, n),
                   approximate_legendre_root(n, k));
}

func weight(n, r) { 2 / ((1 - r**2) * legendre_prime(n, r)**2) }

func nodes(n) {
    gather {
        take(Pair(0, weight(n, 0))) if n.is_odd;
        for i in (1 .. n>>1) {
            var r = legendre_root(n, i);
            var w = weight(n, r);
            take(Pair(r, w), Pair(-r, w));
        }
    }
}

func quadrature(n, f, a, b, nds = nodes(n)) {
    func scale(x) { (x*(b - a) + a + b) / 2 }
    (b - a) / 2 * nds.map{ .second * f(scale(.first)) }.sum
}

[@|(5..10), 20].each { |i|
    printf("Gauss-Legendre %2d-point quadrature ∫₋₃⁺³ exp(x) dx ≈ %.15f\n",
        i, quadrature(i, {.exp}, -3, +3))
}