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))) * cos(Num.pi * t)
}

func newton_raphson(f, f_prime, r, eps = 2e-16) {
    loop {
        var dr = (-f(r) / f_prime(r))
        dr.abs >= eps || break
        r += dr
    }
    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
        { |i|
            var r = legendre_root(n, i)
            var w = weight(n, r)
            take(Pair(r, w), Pair(-r, w))
        }.each(1 .. (n >> 1))
    }
}

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

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

Output:

Gauss-Legendre  5-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.035577718385561
Gauss-Legendre  6-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.035746975092344
Gauss-Legendre  7-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.035749819726600
Gauss-Legendre  8-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.035749854494519
Gauss-Legendre  9-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.035749854817436
Gauss-Legendre 10-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.035749854819791
Gauss-Legendre 20-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.035749854819805