Notebooks · posts/gauss-legendre-quadrature-and-where-tanh-sinh-wins

Gauss-Legendre quadrature against tanh-sinh on a smooth and a singular integral

Adam Zettel · September 25, 2026 · written with an agent

Five Gauss-Legendre points integrate every polynomial of degree 9 or less exactly, and miss degree 10 by exactly what the error formula says. On $e^{x}$ seven points reach rounding. On $1/ 1 - x^{2}$ , which is infinite at both ends of the interval, 128 points still miss by $1.4 \times 1 0^{- 2}$ , and the error only halves when the points double. The tanh-sinh rule gets that integral to rounding with 37 points.

\int_{- 1}^{1} p (x) d x = i = 1 \sum n w_{i} p (x_{i}) for every polynomial p of degree \leq 2 n - 1

A quadrature rule replaces an integral by a weighted sum of the integrand's values: $\int_{- 1}^{1} f (x) d x \approx \sum_{i} w_{i} f (x_{i})$ . With $n$ points there are $2 n$ numbers to choose, the nodes $x_{i}$ and the weights $w_{i}$ . The Gauss-Legendre rule chooses them so that the sum is exact for every polynomial of degree up to $2 n - 1$ , which is as high as $2 n$ free numbers can reach. The nodes turn out to be the roots of the Legendre polynomial $P_{n}$ .

The cells below compute the nodes and weights, check the degree claim, and then try the rule on two integrands: $e^{x}$ , which is smooth, and $1/ 1 - x^{2}$ , whose integral is $π$ although the integrand is infinite at $x = \pm 1$ . On the second one a different rule, tanh-sinh, does far better.

The reference is posts/refs/gauss-legendre-quadrature-and-where-tanh-sinh-wins.py (output beside it): numpy's leggauss, which gets the nodes from the eigenvalues of a tridiagonal matrix rather than by Newton's method as here, closed-form integrals, and the tanh-sinh rule written in numpy.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 30–92

PI = 4.0 * atan(1.0)
#: NEWTON_TOL  stop Newton's method once a step is this small
NEWTON_TOL = 1.0e-15
#: NEWTON_MAX  a cap on Newton steps per node; each node needs about five
NEWTON_MAX = 50

#| [Pₙ(x), Pₙ'(x)]. Pₙ from the three-term recurrence
#| k·Pₖ = (2k − 1)·x·Pₖ₋₁ − (k − 1)·Pₖ₋₂, starting at P₀ = 1,
#| P₁ = x; the slope from (x² − 1)·Pₙ' = n·(x·Pₙ − Pₙ₋₁).
def legendre_and_slope(n, x) {
    p0 = 1.0
    p1 = x
    for k from 2 to n + 1 {
        p2 = ((2.0 * k - 1.0) * x * p1 - (k - 1.0) * p0) / k
        p0 = p1
        p1 = p2
    }
    return [p1, n * (x * p1 - p0) / (x * x - 1.0)]
}

#| The n-point Gauss-Legendre rule, into `xs` (ascending) and `ws`. Each node
#| is a root of Pₙ, found by Newton's method from the estimate
#| cos(π·(i + 3/4)/(n + 1/2)); its weight is 2/((1 − x²)·Pₙ'(x)²).
def gauss_legendre(n, xs, ws) {
    for i to n {
        x = cos(PI * (i + 0.75) / (n + 0.5))
        dx = 1.0
        steps = 0
        while abs(dx) > NEWTON_TOL and steps < NEWTON_MAX {
            ps = legendre_and_slope(n, x)
            dx = ps[0] / ps[1]
            x = x - dx
            steps += 1
        }
        ps = legendre_and_slope(n, x)
        xs[n - 1 - i] = x
        ws[n - 1 - i] = 2.0 / ((1.0 - x * x) * ps[1] * ps[1])
    }
}

#: N5  the rule the title is about
N5 = 5
X5 = zeros(N5)
W5 = zeros(N5)
gauss_legendre(N5, X5, W5)
show X5
show W5

# numpy.polynomial.legendre.leggauss(5): posts/refs/gauss-legendre-quadrature-and-where-tanh-sinh-wins.out
NUMPY_X5 = [-0.9061798459386641, -0.5384693101056831, 0.0, 0.5384693101056831, 0.9061798459386641]
NUMPY_W5 = [0.23692688505618908, 0.47862867049936625, 0.568888888888889, 0.47862867049936625, 0.23692688505618908]
node_diff = max(abs(X5 - NUMPY_X5))
weight_diff = max(abs(W5 - NUMPY_W5))
show node_diff
show weight_diff

test "the five nodes and weights are numpy's leggauss(5), to 3e-16" {
    assert_lt(node_diff, 3.0e-16)
    assert_lt(weight_diff, 3.0e-16)
}
test "the weights add up to 2, the length of the interval" {
    assert_near(sum(W5), 2.0, 3.0e-16)
}

PI = 4.0 * atan(1.0)
#: NEWTON_TOL  stop Newton's method once a step is this small
NEWTON_TOL = 1.0e-15
#: NEWTON_MAX  a cap on Newton steps per node; each node needs about five
NEWTON_MAX = 50

#| [Pₙ(x), Pₙ'(x)]. Pₙ from the three-term recurrence
#| k·Pₖ = (2k − 1)·x·Pₖ₋₁ − (k − 1)·Pₖ₋₂, starting at P₀ = 1,
#| P₁ = x; the slope from (x² − 1)·Pₙ' = n·(x·Pₙ − Pₙ₋₁).
def legendre_and_slope(n, x) {
    p0 = 1.0
    p1 = x
    for k from 2 to n + 1 {
        p2 = ((2.0 * k - 1.0) * x * p1 - (k - 1.0) * p0) / k
        p0 = p1
        p1 = p2
    }
    return [p1, n * (x * p1 - p0) / (x * x - 1.0)]
}

#| The n-point Gauss-Legendre rule, into `xs` (ascending) and `ws`. Each node
#| is a root of Pₙ, found by Newton's method from the estimate
#| cos(π·(i + 3/4)/(n + 1/2)); its weight is 2/((1 − x²)·Pₙ'(x)²).
def gauss_legendre(n, xs, ws) {
    for i to n {
        x = cos(PI * (i + 0.75) / (n + 0.5))
        dx = 1.0
        steps = 0
        while abs(dx) > NEWTON_TOL and steps < NEWTON_MAX {
            ps = legendre_and_slope(n, x)
            dx = ps[0] / ps[1]
            x = x - dx
            steps += 1
        }
        ps = legendre_and_slope(n, x)
        xs[n - 1 - i] = x
        ws[n - 1 - i] = 2.0 / ((1.0 - x * x) * ps[1] * ps[1])
    }
}

#: N5  the rule the title is about
N5 = 5
X5 = zeros(N5)
W5 = zeros(N5)
gauss_legendre(N5, X5, W5)
show X5
show W5

# numpy.polynomial.legendre.leggauss(5): posts/refs/gauss-legendre-quadrature-and-where-tanh-sinh-wins.out
NUMPY_X5 = [-0.9061798459386641, -0.5384693101056831, 0.0, 0.5384693101056831, 0.9061798459386641]
NUMPY_W5 = [0.23692688505618908, 0.47862867049936625, 0.568888888888889, 0.47862867049936625, 0.23692688505618908]
node_diff = max(abs(X5 - NUMPY_X5))
weight_diff = max(abs(W5 - NUMPY_W5))
show node_diff
show weight_diff

test "the five nodes and weights are numpy's leggauss(5), to 3e-16" {
    assert_lt(node_diff, 3.0e-16)
    assert_lt(weight_diff, 3.0e-16)
}
test "the weights add up to 2, the length of the interval" {
    assert_near(sum(W5), 2.0, 3.0e-16)
}

Exact through degree 9

The cell integrates $(1 + x)^{k}$ for $k = 0, 1, \dots, 13$ , whose exact integral is $2^{k + 1} / (k + 1)$ , with the 5-point rule and with the 3-point rule. (A power of $x$ alone would be a poor test: for odd $k$ the symmetric nodes give 0, which is right, for any rule.) The 5-point rule is exact to rounding through degree 9 and then misses. The 3-point rule misses from degree 6 on, as $2 n - 1 = 5$ says.

At degree $2 n$ the miss is known exactly. The error of the $n$ -point rule is

\int_{- 1}^{1} f d x - i \sum w_{i} f (x_{i}) = \frac{2 ^{2 n + 1} ( n ! ) ^{4}}{( 2 n + 1 ) (( 2 n )! ) ^{3}} f^{(2 n)} (ξ)

for some $ξ$ in the interval, and for $f = (1 + x)^{10}$ the tenth derivative is the constant $10!$ , so with $n = 5$ the error is $2^{11} (5!)^{4} / (11 (10!)^{2}) = 2.9318 \times 1 0^{- 3}$ .

In the figure an error of exactly 0 is drawn at $1 0^{- 17}$ , since a log axis has no 0.

[2]lines 116–184

#: DEGREES  the powers k = 0 .. 13 of (1 + x) that are integrated
DEGREES = 14
#: DEGREE_TOP  the highest power integrated
DEGREE_TOP = DEGREES - 1.0
#: FLOOR  where an error of exactly zero is drawn on a log axis
FLOOR = 1.0e-17

#| The relative error of the rule (`xs`, `ws`) on ∫ (1 + x)^k dx over [−1, 1],
#| for k = 0 .. DEGREES - 1.
def power_errors(xs, ws) {
    errs = zeros(DEGREES)
    for k to DEGREES {
        exact = pow(2.0, k + 1.0) / (k + 1.0)
        errs[k] = abs(sum(ws * pow(1.0 + xs, k)) - exact) / exact
    }
    return errs
}

#: N3  a smaller rule, to show the cliff moving with n
N3 = 3
X3 = zeros(N3)
W3 = zeros(N3)
gauss_legendre(N3, X3, W3)
ERR5 = power_errors(X5, W5)
ERR3 = power_errors(X3, W3)
show ERR5
show ERR3

#| m!
def factorial(m) {
    p = 1.0
    for k from 2 to m + 1 {
        p = p * k
    }
    return p
}

#: DEG10  the first degree the 5-point rule gets wrong
DEG10 = 10
miss10 = pow(2.0, DEG10 + 1.0) / (DEG10 + 1.0) - sum(W5 * pow(1.0 + X5, DEG10))
formula10 = pow(2.0, DEG10 + 1.0) * pow(factorial(N5), 4.0) / ((DEG10 + 1.0) * pow(factorial(DEG10), 2.0))
show miss10
show formula10

degree_axis = linspace(0.0, DEGREES - 1.0, DEGREES)
figure("Relative error on the integral of $(1 + x)^k$", "degree k", "relative error")
yscale("log")
plot(degree_axis, max(ERR5, FLOOR), "5 points")
plot(degree_axis, max(ERR3, FLOOR), "3 points")
style("5 points", "s1 circle")
style("3 points", "s2 circle")

test "the 5-point rule is exact to rounding (under 1e-15) for every degree up to 9" {
    for k to 2 * N5 {
        assert_lt(ERR5[k], 1.0e-15)
    }
}
test "and degree 10 misses by 2^11 (5!)^4 / (11 (10!)^2) = 2.9318e-3, as numpy's rule does" {
    assert_near(miss10, 0.0029318124556219796, 1.0e-12)
    assert_near(formula10, 0.0029318124556219796, 1.0e-17)
    # numpy's leggauss(5) on the same integral
    assert_near(miss10, 0.002931812455585714, 1.0e-12)
}
test "the 3-point rule is exact through degree 5 and off by 2.50e-3 (relative) at degree 6, as numpy finds" {
    for k to 2 * N3 {
        assert_lt(ERR3[k], 1.0e-15)
    }
    assert_near(ERR3[2 * N3], 2.50e-3, 0.005e-3)
}

#: DEGREES  the powers k = 0 .. 13 of (1 + x) that are integrated
DEGREES = 14
#: DEGREE_TOP  the highest power integrated
DEGREE_TOP = DEGREES - 1.0
#: FLOOR  where an error of exactly zero is drawn on a log axis
FLOOR = 1.0e-17

#| The relative error of the rule (`xs`, `ws`) on ∫ (1 + x)^k dx over [−1, 1],
#| for k = 0 .. DEGREES - 1.
def power_errors(xs, ws) {
    errs = zeros(DEGREES)
    for k to DEGREES {
        exact = pow(2.0, k + 1.0) / (k + 1.0)
        errs[k] = abs(sum(ws * pow(1.0 + xs, k)) - exact) / exact
    }
    return errs
}

#: N3  a smaller rule, to show the cliff moving with n
N3 = 3
X3 = zeros(N3)
W3 = zeros(N3)
gauss_legendre(N3, X3, W3)
ERR5 = power_errors(X5, W5)
ERR3 = power_errors(X3, W3)
show ERR5
show ERR3

#| m!
def factorial(m) {
    p = 1.0
    for k from 2 to m + 1 {
        p = p * k
    }
    return p
}

#: DEG10  the first degree the 5-point rule gets wrong
DEG10 = 10
miss10 = pow(2.0, DEG10 + 1.0) / (DEG10 + 1.0) - sum(W5 * pow(1.0 + X5, DEG10))
formula10 = pow(2.0, DEG10 + 1.0) * pow(factorial(N5), 4.0) / ((DEG10 + 1.0) * pow(factorial(DEG10), 2.0))
show miss10
show formula10

degree_axis = linspace(0.0, DEGREES - 1.0, DEGREES)
figure("Relative error on the integral of $(1 + x)^k$", "degree k", "relative error")
yscale("log")
plot(degree_axis, max(ERR5, FLOOR), "5 points")
plot(degree_axis, max(ERR3, FLOOR), "3 points")
style("5 points", "s1 circle")
style("3 points", "s2 circle")

test "the 5-point rule is exact to rounding (under 1e-15) for every degree up to 9" {
    for k to 2 * N5 {
        assert_lt(ERR5[k], 1.0e-15)
    }
}
test "and degree 10 misses by 2^11 (5!)^4 / (11 (10!)^2) = 2.9318e-3, as numpy's rule does" {
    assert_near(miss10, 0.0029318124556219796, 1.0e-12)
    assert_near(formula10, 0.0029318124556219796, 1.0e-17)
    # numpy's leggauss(5) on the same integral
    assert_near(miss10, 0.002931812455585714, 1.0e-12)
}
test "the 3-point rule is exact through degree 5 and off by 2.50e-3 (relative) at degree 6, as numpy finds" {
    for k to 2 * N3 {
        assert_lt(ERR3[k], 1.0e-15)
    }
    assert_near(ERR3[2 * N3], 2.50e-3, 0.005e-3)
}

A smooth integrand and a singular one

$e^{x}$ is not a polynomial, but it is smooth everywhere, including off the real line, and polynomials approximate it well. The Gauss-Legendre error falls faster than any power of $n$ : $8.2 \times 1 0^{- 10}$ at 5 points, rounding at 7.

$1/ 1 - x^{2}$ is a different case. Its integral is $π$ , finite, but the function goes to infinity at both ends, and no polynomial follows it there. The rule never evaluates at $\pm 1$ , which is good, since the value there is infinite; but its outermost node sits a distance of order $1/ n^{2}$ from the end, and the integral over that last stretch, $\int d x / 2 (1 - x)$ over a length of order $1/ n^{2}$ , is of order $1/ n$ . The error behaves like that: it halves each time $n$ doubles, and $n$ times the error settles near 1.73. At that rate each further digit costs ten times as many points.

[3]lines 203–247

E_EXACT = exp(1.0) - exp(-1.0)
#: N_SMOOTH  Gauss-Legendre rules from 1 to 12 points on e^x
N_SMOOTH = 12
gl_smooth_n = linspace(1.0, N_SMOOTH, N_SMOOTH)
gl_smooth_err = zeros(N_SMOOTH)
for i to N_SMOOTH {
    n = i + 1
    xs = zeros(n)
    ws = zeros(n)
    gauss_legendre(n, xs, ws)
    gl_smooth_err[i] = abs(sum(ws * exp(xs)) - E_EXACT)
}
show gl_smooth_err

#: SING_N  Gauss-Legendre rules on the singular integrand: 4, 8, ..., 128 points
SING_N = [4, 8, 16, 32, 64, 128]
NSING = len(SING_N)
gl_sing_n = zeros(NSING)
gl_sing_err = zeros(NSING)
for i to NSING {
    n = SING_N[i]
    xs = zeros(n)
    ws = zeros(n)
    gauss_legendre(n, xs, ws)
    gl_sing_n[i] = n
    gl_sing_err[i] = abs(sum(ws / sqrt(1.0 - xs * xs)) - PI)
}
show gl_sing_err
show gl_sing_err * gl_sing_n

#: GL_ROUNDING  the number of Gauss-Legendre points at which e^x reaches rounding
GL_ROUNDING = 7
test "Gauss-Legendre misses the integral of e^x by 8.248e-10 with 5 points and by under 3e-15 from 7 on, as with numpy's nodes" {
    assert_near(gl_smooth_err[4], 8.248e-10, 0.001e-10)
    for i from GL_ROUNDING - 1 to N_SMOOTH {
        assert_lt(gl_smooth_err[i], 3.0e-15)
    }
}
test "128 Gauss-Legendre points still miss the integral of 1/sqrt(1 - x^2) by 1.3552e-2, numpy's figure" {
    assert_near(gl_sing_err[NSING - 1], 1.355156e-2, 1.0e-7)
}
test "and n times the error is 1.7279 at 64 points and 1.7346 at 128, as numpy finds" {
    assert_near(gl_sing_err[NSING - 2] * gl_sing_n[NSING - 2], 1.7279, 1.0e-4)
    assert_near(gl_sing_err[NSING - 1] * gl_sing_n[NSING - 1], 1.7346, 1.0e-4)
}

E_EXACT = exp(1.0) - exp(-1.0)
#: N_SMOOTH  Gauss-Legendre rules from 1 to 12 points on e^x
N_SMOOTH = 12
gl_smooth_n = linspace(1.0, N_SMOOTH, N_SMOOTH)
gl_smooth_err = zeros(N_SMOOTH)
for i to N_SMOOTH {
    n = i + 1
    xs = zeros(n)
    ws = zeros(n)
    gauss_legendre(n, xs, ws)
    gl_smooth_err[i] = abs(sum(ws * exp(xs)) - E_EXACT)
}
show gl_smooth_err

#: SING_N  Gauss-Legendre rules on the singular integrand: 4, 8, ..., 128 points
SING_N = [4, 8, 16, 32, 64, 128]
NSING = len(SING_N)
gl_sing_n = zeros(NSING)
gl_sing_err = zeros(NSING)
for i to NSING {
    n = SING_N[i]
    xs = zeros(n)
    ws = zeros(n)
    gauss_legendre(n, xs, ws)
    gl_sing_n[i] = n
    gl_sing_err[i] = abs(sum(ws / sqrt(1.0 - xs * xs)) - PI)
}
show gl_sing_err
show gl_sing_err * gl_sing_n

#: GL_ROUNDING  the number of Gauss-Legendre points at which e^x reaches rounding
GL_ROUNDING = 7
test "Gauss-Legendre misses the integral of e^x by 8.248e-10 with 5 points and by under 3e-15 from 7 on, as with numpy's nodes" {
    assert_near(gl_smooth_err[4], 8.248e-10, 0.001e-10)
    for i from GL_ROUNDING - 1 to N_SMOOTH {
        assert_lt(gl_smooth_err[i], 3.0e-15)
    }
}
test "128 Gauss-Legendre points still miss the integral of 1/sqrt(1 - x^2) by 1.3552e-2, numpy's figure" {
    assert_near(gl_sing_err[NSING - 1], 1.355156e-2, 1.0e-7)
}
test "and n times the error is 1.7279 at 64 points and 1.7346 at 128, as numpy finds" {
    assert_near(gl_sing_err[NSING - 2] * gl_sing_n[NSING - 2], 1.7279, 1.0e-4)
    assert_near(gl_sing_err[NSING - 1] * gl_sing_n[NSING - 1], 1.7346, 1.0e-4)
}

Tanh-sinh

The tanh-sinh rule, from Takahasi and Mori (1974), changes variables first:

x = tanh (\frac{π}{2} sinh t), d x = \frac{π}{2} cosh t sech^{2} (\frac{π}{2} sinh t) d t .

As $t$ runs over the whole line, $x$ runs over $(- 1, 1)$ , and the endpoints move off to $t = \pm \infty$ . The factor $d x / d t$ falls like $exp (- \frac{π}{2} e^{∣ t ∣})$ , a double exponential, which is faster than $1/ 1 - x^{2}$ grows. So the new integrand is smooth and dies off quickly in both directions, and for such a function the plain trapezoid rule, equal steps $h$ in $t$ , converges very fast. Here the nodes are $t = k h$ with $∣ t ∣ \leq 4.5$ . The last term kept on the singular integrand is $2.1 \times 1 0^{- 17}$ at $h = 1$ , whose last node is $t = 4$ , and under $1 0^{- 28}$ for every smaller step; the terms left out are smaller still.

There is one trap. Once $\frac{π}{2} sinh t$ passes about 19, $tanh$ rounds to exactly 1, and $1 - x^{2}$ computed from $x$ is 0: the singular integrand becomes $1/0$ . With $h = 1/4$ that happens at 12 of the 37 nodes, and the sum is infinite. The fix is to compute $1 - x^{2}$ directly as $sech^{2} (\frac{π}{2} sinh t)$ , which is tiny there but not 0, and to write the integrand in terms of it.

Halving $h$ doubles the points. On the singular integral the error goes $1.9 \times 1 0^{- 3}$ , $2.0 \times 1 0^{- 8}$ and about $1 0^{- 15}$ , which is rounding, at 9, 19 and 37 points: the number of correct digits more than doubles with each halving, until rounding stops it. Gauss-Legendre with 128 points was at $1.4 \times 1 0^{- 2}$ . On $e^{x}$ the roles reverse: tanh-sinh needs 73 points for rounding where Gauss-Legendre needed 7. The substitution spends most of its points near the ends, and a smooth integrand has nothing there that needs them.

[4]lines 283–390

#: TMAX  the trapezoid rule runs over t in [-TMAX, TMAX]
TMAX = 4.5
#: LEVELS  steps h = 1, 1/2, ..., 1/64
LEVELS = 7

#| The tanh-sinh rule with step `h` on both integrands:
#| [number of points, sum for e^x, sum for 1/sqrt(1 - x^2), the singular
#| integrand's term at the last node, nodes where 1 - x*x rounds to 0].
#| 1 − x² is taken as sech²(u), never computed from x.
def tanh_sinh(h) {
    kmax = floor(TMAX / h)
    s_exp = 0.0
    s_sing = 0.0
    end_term = 0.0
    zeros_1mx2 = 0.0
    for k from -kmax to kmax + 1 {
        t = k * h
        u = 0.5 * PI * sinh(t)
        x = tanh(u)
        sech2 = 1.0 / (cosh(u) * cosh(u))
        w = h * 0.5 * PI * cosh(t) * sech2
        s_exp += w * exp(x)
        s_sing += w / sqrt(sech2)
        if abs(k) == kmax {
            end_term = max(end_term, w / sqrt(sech2))
        }
        if 1.0 - x * x == 0.0 {
            zeros_1mx2 += 1.0
        }
    }
    return [2.0 * kmax + 1.0, s_exp, s_sing, end_term, zeros_1mx2]
}

ts_points = zeros(LEVELS)
ts_exp_err = zeros(LEVELS)
ts_sing_err = zeros(LEVELS)
ts_zeros = zeros(LEVELS)
ts_end = zeros(LEVELS)
ts_step = zeros(LEVELS)
for level to LEVELS {
    ts_step[level] = pow(0.5, level)
    r = tanh_sinh(ts_step[level])
    ts_points[level] = r[0]
    ts_exp_err[level] = abs(r[1] - E_EXACT)
    ts_sing_err[level] = abs(r[2] - PI)
    ts_end[level] = r[3]
    ts_zeros[level] = r[4]
}
show ts_step
show ts_points
show ts_exp_err
show ts_sing_err
show ts_zeros
show ts_end[0]

figure("$1/\\sqrt{1 - x^2}$ on $[-1, 1]$: error against points", "points", "|error|")
xscale("log")
yscale("log")
plot(gl_sing_n, gl_sing_err, "Gauss-Legendre")
plot(ts_points, max(ts_sing_err, FLOOR), "tanh-sinh")
style("Gauss-Legendre", "s1 circle")
style("tanh-sinh", "s2 circle")

figure("$e^x$ on $[-1, 1]$: error against points", "points", "|error|")
xscale("log")
yscale("log")
plot(gl_smooth_n, max(gl_smooth_err, FLOOR), "Gauss-Legendre")
plot(ts_points, max(ts_exp_err, FLOOR), "tanh-sinh")
style("Gauss-Legendre", "s1 circle")
style("tanh-sinh", "s2 circle")

test "h = 1, 1/2, 1/4 use 9, 19 and 37 points; the last term is 2.073e-17 at h = 1, under 1e-28 after" {
    assert_near(ts_step[1], 0.5, 0.0)
    assert_near(ts_step[2], 0.25, 0.0)
    assert_near(ts_points[0], 9, 0.0)
    assert_near(ts_points[1], 19, 0.0)
    assert_near(ts_points[2], 37, 0.0)
    # numpy's rule: posts/refs/gauss-legendre-quadrature-and-where-tanh-sinh-wins.out
    assert_near(ts_end[0], 2.073e-17, 0.001e-17)
    for level from 1 to LEVELS {
        assert_lt(ts_end[level], 1.0e-28)
    }
}
test "with h = 1/4, 1 - x*x rounds to 0 at 12 of the 37 nodes, as in numpy" {
    assert_near(ts_step[2], 0.25, 0.0)
    assert_near(ts_zeros[2], 12, 0.0)
}
test "tanh-sinh misses the integral of 1/sqrt(1 - x^2) by 1.915e-3 and 1.972e-8 at 9 and 19 points, as numpy's rule does" {
    assert_near(ts_sing_err[0], 1.915e-3, 0.001e-3)
    assert_near(ts_sing_err[1], 1.972e-8, 0.001e-8)
}
test "and under 2e-15 from 37 points on, where Gauss-Legendre at 128 points is off by 1.4e-2" {
    for level from 2 to LEVELS {
        assert_lt(ts_sing_err[level], 2.0e-15)
    }
    assert_lt(0.01, gl_sing_err[NSING - 1])
}
#: TS_ROUNDING  the level (h = 1/8, 73 points) at which tanh-sinh reaches rounding on e^x
TS_ROUNDING = 3
test "tanh-sinh misses the integral of e^x by 8.14e-11 at 37 points and by under 4e-15 from 73 on" {
    assert_near(ts_exp_err[2], 8.136e-11, 0.01e-11)
    for level from TS_ROUNDING to LEVELS {
        assert_lt(ts_exp_err[level], 4.0e-15)
    }
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

#: TMAX  the trapezoid rule runs over t in [-TMAX, TMAX]
TMAX = 4.5
#: LEVELS  steps h = 1, 1/2, ..., 1/64
LEVELS = 7

#| The tanh-sinh rule with step `h` on both integrands:
#| [number of points, sum for e^x, sum for 1/sqrt(1 - x^2), the singular
#| integrand's term at the last node, nodes where 1 - x*x rounds to 0].
#| 1 − x² is taken as sech²(u), never computed from x.
def tanh_sinh(h) {
    kmax = floor(TMAX / h)
    s_exp = 0.0
    s_sing = 0.0
    end_term = 0.0
    zeros_1mx2 = 0.0
    for k from -kmax to kmax + 1 {
        t = k * h
        u = 0.5 * PI * sinh(t)
        x = tanh(u)
        sech2 = 1.0 / (cosh(u) * cosh(u))
        w = h * 0.5 * PI * cosh(t) * sech2
        s_exp += w * exp(x)
        s_sing += w / sqrt(sech2)
        if abs(k) == kmax {
            end_term = max(end_term, w / sqrt(sech2))
        }
        if 1.0 - x * x == 0.0 {
            zeros_1mx2 += 1.0
        }
    }
    return [2.0 * kmax + 1.0, s_exp, s_sing, end_term, zeros_1mx2]
}

ts_points = zeros(LEVELS)
ts_exp_err = zeros(LEVELS)
ts_sing_err = zeros(LEVELS)
ts_zeros = zeros(LEVELS)
ts_end = zeros(LEVELS)
ts_step = zeros(LEVELS)
for level to LEVELS {
    ts_step[level] = pow(0.5, level)
    r = tanh_sinh(ts_step[level])
    ts_points[level] = r[0]
    ts_exp_err[level] = abs(r[1] - E_EXACT)
    ts_sing_err[level] = abs(r[2] - PI)
    ts_end[level] = r[3]
    ts_zeros[level] = r[4]
}
show ts_step
show ts_points
show ts_exp_err
show ts_sing_err
show ts_zeros
show ts_end[0]

figure("$1/\\sqrt{1 - x^2}$ on $[-1, 1]$: error against points", "points", "|error|")
xscale("log")
yscale("log")
plot(gl_sing_n, gl_sing_err, "Gauss-Legendre")
plot(ts_points, max(ts_sing_err, FLOOR), "tanh-sinh")
style("Gauss-Legendre", "s1 circle")
style("tanh-sinh", "s2 circle")

figure("$e^x$ on $[-1, 1]$: error against points", "points", "|error|")
xscale("log")
yscale("log")
plot(gl_smooth_n, max(gl_smooth_err, FLOOR), "Gauss-Legendre")
plot(ts_points, max(ts_exp_err, FLOOR), "tanh-sinh")
style("Gauss-Legendre", "s1 circle")
style("tanh-sinh", "s2 circle")

test "h = 1, 1/2, 1/4 use 9, 19 and 37 points; the last term is 2.073e-17 at h = 1, under 1e-28 after" {
    assert_near(ts_step[1], 0.5, 0.0)
    assert_near(ts_step[2], 0.25, 0.0)
    assert_near(ts_points[0], 9, 0.0)
    assert_near(ts_points[1], 19, 0.0)
    assert_near(ts_points[2], 37, 0.0)
    # numpy's rule: posts/refs/gauss-legendre-quadrature-and-where-tanh-sinh-wins.out
    assert_near(ts_end[0], 2.073e-17, 0.001e-17)
    for level from 1 to LEVELS {
        assert_lt(ts_end[level], 1.0e-28)
    }
}
test "with h = 1/4, 1 - x*x rounds to 0 at 12 of the 37 nodes, as in numpy" {
    assert_near(ts_step[2], 0.25, 0.0)
    assert_near(ts_zeros[2], 12, 0.0)
}
test "tanh-sinh misses the integral of 1/sqrt(1 - x^2) by 1.915e-3 and 1.972e-8 at 9 and 19 points, as numpy's rule does" {
    assert_near(ts_sing_err[0], 1.915e-3, 0.001e-3)
    assert_near(ts_sing_err[1], 1.972e-8, 0.001e-8)
}
test "and under 2e-15 from 37 points on, where Gauss-Legendre at 128 points is off by 1.4e-2" {
    for level from 2 to LEVELS {
        assert_lt(ts_sing_err[level], 2.0e-15)
    }
    assert_lt(0.01, gl_sing_err[NSING - 1])
}
#: TS_ROUNDING  the level (h = 1/8, 73 points) at which tanh-sinh reaches rounding on e^x
TS_ROUNDING = 3
test "tanh-sinh misses the integral of e^x by 8.14e-11 at 37 points and by under 4e-15 from 73 on" {
    assert_near(ts_exp_err[2], 8.136e-11, 0.01e-11)
    for level from TS_ROUNDING to LEVELS {
        assert_lt(ts_exp_err[level], 4.0e-15)
    }
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors