``````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",
``````Gauss-Legendre  5-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.035577718385561