Fórmulas de Exponenciais Duplas

integração numérica

fortran

python

julia

Integração numérica pela quadratura tanh-sinh.

Autor

Pedro H. N. Vieira

Data de Publicação

5 de dezembro de 2025

O desenvolvimento da Fórmula de Exponencial Dupla (DE), também conhecida como Quadratura tanh-sinh, surge apartir da regra trapezoidal primitiva \[ \int_{-\infty}^{\infty} f(\xi) \, d\xi \approx h \sum_{k=-\infty}^{\infty} f(k h) \] sendo \(h>0\). O erro de discretização é \(\mathcal{O}(\exp(-2\pi w/h))\) quando \(h \to 0\), sendo \(w\) a distância do eixo real para a singularidade mais próxima do integrando. Como a convergência é exponencial, dividindo \(h\) pela metade faz dobrar o número de dígitos corretos no resultado para valores suficientimente pequenos de \(h\) (Takahasi e Mori 1974; Bailey, Jeyabalan, e Li 2005; Vanherck, Sorée, e Magnus 2020; Engelen 2021).

Intervalos Finitos

Para aplicar tal método a um intervalo finito, faz-se um mapa de \((-\infty,\infty)\) para \([-1,1]\) usando uma função monótona \(\varphi(t)\), tal que

\[ \int_{-1}^{1} f(\xi) \, d\xi = \int_{-\infty}^{\infty} f[\varphi(t)] \, \dot{\varphi}(t) \, dt \approx h \sum_{k=-n}^{n} \dot{\varphi}(k h)\,f[\varphi(kh)] \]

onde \(\dot{\varphi} \equiv d\varphi/dt\) é a Jacobiana da transformada. Usa-se uma transformada que force a derivada do integrando a zerar nos pontos extremos (finais). Pode-se escolher

\[ \varphi(t) = \tanh g(t) \\ \dot{\varphi}(t) = \dot{g}(t) \, \text{sech}^2 g(t) \]

onde \(g(t)\) é uma função ímpar que leva a

\[ \int_{-1}^{1} f(\xi) \, d\xi \approx h \sum_{k=-n}^{n} w_k f(\xi_k) \]

sendo

\[ \xi_k = 1 - \delta_k \\ w_k = 2\, \dot{g}(kh) \, \delta_k \, (1 + q_k)^{-1} \\ \delta_k = 2\,q_k\,(1 + q_k)^{-1} \\ q_k = e^{-2\,g(kh)} \]

Usando o fato que \(\xi_{-k} = -\xi_k\) e \(w_{-k} = w_k\), podemos expressar a integração por

\[ \int_{-1}^{1} f(\xi) \, d\xi \approx h \left \{ \dot{g}(0) \, f(0) + \sum_{k=1}^{n} w_k \left [ f(-1 + \delta_k) + f(+1 - \delta_k) \right ] \right \}. \]

Um intervalo arbitrário \([a,b]\) pode ser mapeado em \([-1,1]\) pela transformada linear \[ \xi = \sigma\,x+\gamma \\ \sigma=(b-a)/2 \\ \gamma = (a+b)/2, \]

resultando em:

\[ \int_a^b f(\xi)\,d\xi = \sigma \int_{-1}^1 f(\sigma \, x + \gamma)\,dx \]

\[ \int_{a}^{b} f(\xi) \, d\xi \approx \sigma h \left \{ \dot{g}(0) \, f(\gamma) + \sum_{k=1}^{n} w_k \left [ f(a + \sigma\,\delta_k) + f(b - \sigma\,\delta_k) \right ] \right \}. \]

A função \(g(t)\) pode ser escolhida tal que

\[ g(t) = \eta \sinh t, \quad \dot{g}(t) = \eta \cosh t, \]

a qual dá origem à regra tanh-sinh, onde \(\eta\) é um parâmetro positivo. A Jacobiana \(\dot{g}(t)\) desta transformada tem decaimento duplo-exponencial para \(t \to \pm \infty\). É recomendado escolher \(\eta=1\) por conveniência, pois a transformada é, praticamente, insensível a valores de \(\eta\) próximos da unidade. Para um \(n\) fixo, recomenda-se escolher \(h \approx \log(\pi^2 n) / n\) o qual, para \(\eta = 1\), leva a um erro total da ordem de \(\mathcal{O}(\exp(-\pi^2 n / \log(\pi^2 n)))\).

Algo que torna este método de integração estranho é que ele não é exato para nenhum polinômio, nem sequer para uma constante, contudo, ele é notavelmente acurado. Outra característica deste método é que uma amostragem uniforme no eixo \(t\) leva a um agrupamento dos nós \(\xi_k\) perto dos pontos extremos do intervalo de integração, o qual, devido ao decaimento rápido dos pesos \(w_k\), permite uma supressão de qualquer singularidade nos pontos extremos de \(f(\xi)\).

Para estimar o erro, arbitra-se um tamanho inicial de malha \(h\), dobra-se o número de pontos a cada iteração (fazendo \(h_{k+1} = h_k / 2\)), pois, assim, pode-se reutilizar os pontos calculados na iteração anterior, até que se atinja a convergência especificada \(\epsilon\), isto é, são feitas sucessivas iterações \(k\) até que a condição seguinte seja satisfeita. Tipicamente, o erro real será \(\epsilon^2\) após a convergência.

\[ |T_n| < \epsilon \left | T_0 + \sum_{k=1}^n T_k \right |, \quad n =1,2,3,\dots \]

Intervalos Semi-infinitos

Para um intervalo semi-infinito \([a,\infty)\), uma transformada DE mista apropriada é

\[ \xi \leftarrow \varphi(t) = a + e^{t - \exp(-t)} \\ \dot{\varphi}(t) = (1 + e^{-t})\, e^{t - \exp(-t)} \]

a qual leva a

\[ \int_a^\infty d(\xi) d\xi \approx h \sum_{k=-n_1}^{n_2} w_k \, f(a + \delta_k) \]

com abscissas e pesos dados por

\[ \delta_k = e^{kh - \exp(-kh)} \\ w_k = (1 + e^{-kh}) \, \delta_k \]

Note que, devido à falta de simetria da transformada, o somatório não pode ser “dobrado” em \(k=0\) como para os interfalos finitos. Assim sendo, as somas parciais devem ser calculadas para \(k\) positivo e negativo, separadamente.

Códigos

Para testar o método, faremos as seguintes integrações procurando usando uma tolerância de \(1 \times 10^{-6}\) como critério de convergência.

\[ y_1 = \int_{-1}^{1} \frac{1}{(x - 2)\,(1 - x)^{1/4}\,(1 + x)^{3/4}} \approx -1.9490 \]

\[ y_2 = \int_{-1}^{1} \frac{\cos(\pi \, x)}{\sqrt{1 - x}} \approx -0.69049 \]

\[ y_3 = \int_{0}^{\infty} \frac{\exp(-1 - x)}{1 + x} \approx 0.21938 \]

\[ y_4 = \int_{-\infty}^{\infty} \frac{1}{(1 + x^2)^{5 / 4}} \approx 2.3962 \]

\[ y_5 = \int_{-\infty}^{\infty} \frac{1}{1 + x^4} \approx 2.2214 \]

\[ y_6 = \int_{0.1}^{1} \frac{1}{x^2} = 9 \]

Existe o pacote DoubleExponentialFormulas que podemos usar.

using DoubleExponentialFormulas

# y1
function f1(x)
    return 1.0 / ((x - 2.0) * (1.0 - x)^0.25 * (1.0 + x)^0.75)
end

y1 = quadde(f1, -1.0, 1.0, atol=1e-6)
println("====================================\ny1")
println(y1)

# y2
function f2(x)
    return cos(pi * x) / (1.0 - x)^0.5
end

y2 = quadde(f2, -1.0, 1.0, atol=1e-6)
println("====================================\ny2")
println(y2)

# y3
function f3(x)
    return exp(-1.0 - x) / (1.0 + x)
end

y3 = quadde(f3, 0.0, Inf, atol=1e-6)
println("====================================\ny3")
println(y3)

# y4
function f4(x)
    return 1.0 / (1.0 + x^2)^1.25
end

y4 = quadde(f4, -Inf, Inf, atol=1e-6)
println("====================================\ny4")
println(y4)

# y5
function f5(x)
    return 1.0 / (1.0 + x^4)
end

y5 = quadde(f5, -Inf, Inf, atol=1e-6)
println("====================================\ny5")
println(y5)

# y6
function f6(x)
    return 1.0 / x^2.0
end

y6 = quadde(f6, 0.1, 1.0, atol=1e-6)
println("====================================\ny6")
println(y6)
println("====================================")

====================================
y1 = -1.9489064907098714
erro1 = 4.1439840463189424e-7
====================================
y2 = -0.6904945733550009
erro2 = 2.0692066870700658e-11
====================================
y3 = 0.21938354361828977
erro3 = 4.1047020901104285e-8
====================================
y4 = 2.3962804694711894
erro4 = 1.9725917688128037e-11
====================================
y5 = 2.2214414690791875
erro5 = 6.796168248378037e-13
====================================
y6 = 9.000000000000012
erro6 = 3.0485303698109525e-10
====================================

O SciPy tem uma implementação deste esquema de integração numérica.

from numpy import exp, cos, pi, inf
from scipy.integrate import tanhsinh

# y1
def f1(x):
    return 1.0 / ((x - 2.0) * (1.0 - x)**0.25 * (1.0 + x)**0.75)

y1 = tanhsinh(f1, -1.0, 1.0, atol=1e-6)
print("====================================\ny1")
print(y1)

# y2
def f2(x):
    return cos(pi * x) / (1.0 - x)**0.5

y2 = tanhsinh(f2, -1.0, 1.0, atol=1e-6)
print("====================================\ny2")
print(y2)

# y3
def f3(x):
    return exp(-1.0 - x) / (1.0 + x)

y3 = tanhsinh(f3, 0.0, inf, atol=1e-6)
print("====================================\ny3")
print(y3)

# y4
def f4(x):
    return 1.0 / (1.0 + x**2)**1.25

y4 = tanhsinh(f4, -inf, inf, atol=1e-6)
print("====================================\ny4")
print(y4)

# y5
def f5(x):
    return 1.0 / (1.0 + x**4)

y5 = tanhsinh(f5, -inf, inf, atol=1e-6)
print("====================================\ny5")
print(y5)

# y6
def f6(x):
    return 1.0 / x**2.0

y6 = tanhsinh(f6, 0.1, 1.0, atol=1e-6)
print("====================================\ny6")
print(y6)
print("====================================")

====================================
y1
     success: False
      status: -2
    integral: -1.948955261697972
       error: 0.0010072696028446116
        nfev: 16387
    maxlevel: 10
====================================
y2
     success: True
      status: 0
    integral: -0.6904945815613805
       error: 6.187935331412465e-07
        nfev: 131
    maxlevel: 3
====================================
y3
     success: True
      status: 0
    integral: 0.2193844684452573
       error: 4.288750071284711e-07
        nfev: 67
    maxlevel: 2
====================================
y4
     success: True
      status: 0
    integral: 2.3962804694708084
       error: 2.2912128284329816e-12
        nfev: 131
    maxlevel: 3
====================================
y5
     success: True
      status: 0
    integral: 2.221441469083877
       error: 3.163386563835863e-11
        nfev: 259
    maxlevel: 4
====================================
y6
     success: True
      status: 0
    integral: 8.99999999345614
       error: 2.4112420575291934e-09
        nfev: 67
    maxlevel: 2
====================================

O valor de \(y_1\) falhou em convergir para a tolerância desejada.

O código para realizar a integração numérica foi organizado em um módulo, disponível no Github.

module quadde_module
    use, intrinsic :: ieee_arithmetic
    use, intrinsic :: iso_fortran_env, only: wp => real64  ! working precision
    implicit none

    private  ! All entities are now module-private by default
    public quadde, integrable_function

    abstract interface
        function integrable_function(x) result(y)
            import wp
            real(wp), intent(in) :: x
            complex(wp) :: y
        end function integrable_function
    end interface

contains

    !> Double exponential quadrature of f(x) in the interval [a, b] with at most
    !> n levels (more than 6 usually leads to no improvement on convergence)
    !> with requested error eps.
    function quadde(f, a, b, n, eps) result(v)
        procedure(integrable_function) :: f
        integer :: n, k, mode
        real(wp) :: a, b, c, d, e, h, eps, tol, t, r, x, w, u, eh, sign
        complex(wp) :: p, q, s, v, y, fp, fm
        logical :: is_finite

        if (eps < 0.0) stop "eps must be positive"
        if (n < 0) stop "n must be positive"

        tol = 10.0 * eps
        mode = 0  ! tanh-sinh
        c = 0.0
        d = 1.0
        sign = 1.0
        h = 2.0

        if (b < a) then  ! swap bounds
            x = b
            b = a
            a = x
            sign = -1.0
        end if

        if (ieee_is_finite(a) .and. ieee_is_finite(b)) then
            c = (a + b) * 0.5
            d = (b - a) * 0.5
            x = c
        else if (ieee_is_finite(a)) then
            mode = 1  ! exp-sinh
            c = a
            x = a + d
        else if (ieee_is_finite(b)) then
            mode = 1  ! exp-sinh
            d = -d
            sign = -sign
            c = b
            x = b + d
        else
            mode = 2  ! sinh-sinh
            x = 0.0
        end if

        s = f(x)
        do k = 0, n
            p = 0.0
            fp = 0.0
            fm = 0.0

            h = h * 0.5
            eh = exp(h)
            t = eh

            if (k > 0) eh = eh * eh

            if (mode == 0) then  ! tanh-sinh
                do
                    u = exp(1.0 / t - t)  ! = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2
                    r = 2.0 * u / (1.0 + u)  ! = 1 - tanh(sinh(j*h))
                    w = (t + 1.0 / t) * r / (1.0 + u)  ! = cosh(j*h) / cosh(sinh(j*h))^2
                    x = d * r

                    if (a + x > a) then  ! too close to a then reuse previous fp
                        y = f(a + x)
                        is_finite = ieee_is_finite(real(y)) .and. ieee_is_finite(aimag(y))
                        if (is_finite) fp = y
                    end if

                    if (b - x < b) then  ! too close to b then reuse previous fm
                        y = f(b - x)
                        is_finite = ieee_is_finite(real(y)) .and. ieee_is_finite(aimag(y))
                        if (is_finite) fm = y
                    end if

                    q = w * (fp + fm)
                    p = p + q
                    t = t * eh

                    if (abs(q) < eps * abs(p)) exit
                end do
            else
                t = t * 0.5
                do
                    r = exp(t - 0.25 / t)  ! exp(sinh(j*h))
                    w = r
                    q = 0.0
                    if (mode == 1) then  ! exp-sinh
                        x = c + d / r
                        if (abs(x - c) < epsilon(x)) exit  ! x hit the finite endpoint
                        y = f(x)
                        is_finite = ieee_is_finite(real(y)) .and. ieee_is_finite(aimag(y))
                        if (is_finite) q = q + y / w
                    else  ! sinh-sinh
                        r = (r - 1.0 / r) * 0.5  ! sinh(sinh(j*h))
                        w = (w + 1.0 / w) * 0.5  ! cosh(sinh(j*h))
                        x = c - d * r
                        y = f(x)
                        if (is_finite) q = q + y * w
                    end if
                    x = c + d * r
                    y = f(x)
                    if (is_finite) q = q + y * w
                    q = q * (t + 0.25 / t)  ! cosh(j*h)
                    p = p + q
                    t = t * eh
                    if (abs(q) <= eps * abs(p)) exit
                end do
            end if
            v = s - p
            s = s + p
            if (abs(v) <= tol * abs(s)) exit
        end do
        e = abs(v) / (abs(s) + eps)  ! TODO return or print e
        v = sign * d * s * h  ! result with estimated relative error e
    end function quadde

end module

Para testar o código, fizemos o seguinte programa, também disponível no Github.

program quadde_test
    use quadde_module
    use, intrinsic :: ieee_arithmetic
    use, intrinsic :: iso_fortran_env, only: wp => real64  ! working precision
    implicit none

    real(wp), parameter :: tol = 1e-6
    real(wp), parameter :: PI = 3.141592653589793238462643383279502884197169399375105820974944592307816406286198_wp
    real(wp) :: POS_INF, NEG_INF
    complex(wp) :: val
    logical :: is_nan

    POS_INF = ieee_value(1.0_wp, ieee_positive_inf)
    NEG_INF = ieee_value(-1.0_wp, ieee_negative_inf)

    print *, "Testing quadde..."

    val = quadde(f1, -1.0_wp, 1.0_wp, 6, tol)
    print *, "Test 1:", val
    is_nan = ieee_is_nan(real(val)) .or. ieee_is_nan(aimag(val))
    if (abs(-1.9490 - val) > 1e-4 .or. is_nan) stop "Test 1 failed"

    val = quadde(f2, -1.0_wp, 1.0_wp, 6, tol)
    print *, "Test 2:", val
    is_nan = ieee_is_nan(real(val)) .or. ieee_is_nan(aimag(val))
    if (abs(-0.69049 - val) > 1e-4 .or. is_nan) stop "Test 2 failed"

    val = quadde(f3, 0.0_wp, POS_INF, 6, tol)
    print *, "Test 3:", val
    is_nan = ieee_is_nan(real(val)) .or. ieee_is_nan(aimag(val))
    if (abs(0.21938 - val) > 1e-4 .or. is_nan) stop "Test 3 failed"

    val = quadde(f4, NEG_INF, POS_INF, 6, tol)
    print *, "Test 4:", val
    is_nan = ieee_is_nan(real(val)) .or. ieee_is_nan(aimag(val))
    if (abs(2.3962 - val) > 1e-4 .or. is_nan) stop "Test 4 failed"

    val = quadde(f5, NEG_INF, POS_INF, 6, tol)
    print *, "Test 5:", val
    is_nan = ieee_is_nan(real(val)) .or. ieee_is_nan(aimag(val))
    if (abs(2.2214 - val) > 1e-4 .or. is_nan) stop "Test 5 failed"

    val = quadde(f6, 0.1_wp, 1.0_wp, 6, tol)
    print *, "Test 6:", val
    is_nan = ieee_is_nan(real(val)) .or. ieee_is_nan(aimag(val))
    if (abs(9.0_wp - val) > 1e-4 .or. is_nan) stop "Test 6 failed"

    print *, "All tests passed!"

contains
    ! After contains, you can only define subroutines or functions.
    ! You cannot put assignments, if statements, or any "straight-line" code there.

    function f1(x) result(y)
        real(wp), intent(in) :: x
        complex(wp) :: y
        y = 1.0 / ( (x - 2.0) * (1.0 - x)**0.25 * (1.0 + x)**0.75 )
    end function


    function f2(x) result(y)
        real(wp), intent(in) :: x
        complex(wp) :: y
        y = cos(PI * x) / sqrt(1.0 - x)
    end function


    function f3(x) result(y)
        real(wp), intent(in) :: x
        complex(wp) :: y
        y = exp(-1.0 - x) / (1.0 + x)
    end function


    function f4(x) result(y)
        real(wp), intent(in) :: x
        complex(wp) :: y
        y = 1.0 / (1.0 + x**2.0)**(5.0 / 4.0)
    end function


    function f5(x) result(y)
        real(wp), intent(in) :: x
        complex(wp) :: y
        y = 1.0 / (1.0 + x**4.0)
    end function


    function f6(x) result(y)
        real(wp), intent(in) :: x
        complex(wp) :: y
        y = x**(-2.0)
    end function

end program quadde_test

Compilar estes códigos e rodar pode ser feito com os seguintes comandos no terminal:

gfortran -fPIC -c quadde_module.f90
gfortran -fPIC -o test quadde_test.f90 quadde_module.o
./test

Cuja saída é:

 Testing quadde...
 Test 1:              (-1.9489702845825310,0.0000000000000000)
 Test 2:             (-0.69049458018520249,0.0000000000000000)
 Test 3:              (0.21938393436009282,0.0000000000000000)
 Test 4:               (2.3962804694707533,0.0000000000000000)
 Test 5:               (2.2214414406307439,0.0000000000000000)
 Test 6:               (8.9999999646063316,0.0000000000000000)
 All tests passed!

O código C, também no Github foi adaptado do apresentado no anexo de Engelen (2021).

#include "stdio.h"
#include "float.h"
#include "math.h"

// integrate function f, range a..b, max levels n, error tolerance eps
double quadde(double (*f)(double), double a, double b, int n, double eps) {
    const double tol = 10*eps;
    double c = 0, d = 1, s, sign = 1, e, v, h = 2;
    int k = 0, mode = 0; // Tanh-Sinh = 0, Exp-Sinh = 1, Sinh-Sinh = 2
    if (b < a) { // swap bounds
        v = b;
        b = a;
        a = v;
        sign = -1;
    }
    if (isfinite(a) && isfinite(b)) {
        c = (a+b)/2;
        d = (b-a)/2;
        v = c;
    }
    else if (isfinite(a)) {
        mode = 1; // Exp-Sinh
        // alternatively d = exp_sinh_opt_d(f, a, eps, d);
        c = a;
        v = a+d;
    }
    else if (isfinite(b)) {
        mode = 1; // Exp-Sinh
        d = -d; // alternatively d = exp_sinh_opt_d(f, b, eps, -d);
        sign = -sign;
        c = b;
        v = b+d;
    }
    else {
        mode = 2; // Sinh-Sinh
        v = 0;
    }
    s = f(v);
    do {
        double p = 0, q, fp = 0, fm = 0, t, eh;
        h /= 2;
        t = eh = exp(h);
        if (k > 0)
            eh *= eh;
        if (mode == 0) { // Tanh-Sinh
            do {
                double u = exp(1/t-t); // = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2
                double r = 2*u/(1+u); // = 1 - tanh(sinh(j*h))
                double w = (t+1/t)*r/(1+u); // = cosh(j*h)/cosh(sinh(j*h))^2
                double x = d*r;
                if (a+x > a) { // if too close to a then reuse previous fp
                    double y = f(a+x);
                    if (isfinite(y))
                        fp = y; // if f(x) is finite, add to local sum
                }
                if (b-x < b) { // if too close to a then reuse previous fp
                    double y = f(b-x);
                    if (isfinite(y))
                        fm = y; // if f(x) is finite, add to local sum
                }
                q = w*(fp+fm);
                p += q;
                t *= eh;
            } while (fabs(q) > eps*fabs(p));
        }
        else {
            t /= 2;
            do {
                double r = exp(t-.25/t); // = exp(sinh(j*h))
                double x, y, w = r;
                q = 0;
                if (mode == 1) { // Exp-Sinh
                    x = c + d/r;
                    if (fabs(x - c) < DBL_EPSILON) // if x hit the finite endpoint then break
                        break;
                    y = f(x);
                    if (isfinite(y)) // if f(x) is finite, add to local sum
                        q += y/w;
                }
                else { // Sinh-Sinh
                    r = (r-1/r)/2; // = sinh(sinh(j*h))
                    w = (w+1/w)/2; // = cosh(sinh(j*h))
                    x = c - d*r;
                    y = f(x);
                    if (isfinite(y)) // if f(x) is finite, add to local sum
                        q += y*w;
                }
                x = c + d*r;
                y = f(x);
                if (isfinite(y)) // if f(x) is finite, add to local sum
                    q += y*w;
                q *= t+.25/t; // q *= cosh(j*h)
                p += q;
                t *= eh;
            } while (fabs(q) > eps*fabs(p));
        }
        v = s-p;
        s += p;
        ++k;
    } while (fabs(v) > tol*fabs(s) && k <= n);
    e = fabs(v)/(fabs(s)+eps);
    return sign*d*s*h; // result with estimated relative error e
}



double f1(double x) {
    return 1.0 / ( (x - 2.0) * pow(1.0 - x, 0.25) * pow(1.0 + x, 0.75) );
}

double f2(double x) {
    return cos(M_PI * x) / sqrt(1.0 - x);
}

double f3(double x) {
    return exp(-1.0 - x) / (1.0 + x);
}

double f4(double x) {
    return 1.0 / pow(1.0 + pow(x, 2.0), (5.0 / 4.0));
}

double f5(double x) {
    return 1.0 / (1.0 + pow(x, 4.0));
}

double f6(double x) {
    return 1.0 / pow(x, 2.0);
}

int main() {
    double tol = 1e-6;
    double val;

    val = quadde(f1, -1.0, 1.0, 6, tol);
    printf("Test 1: %.17g\n", val);
    if (fabs(-1.9490 - val) > 1e-4) {
        printf("Test 1 failed");
        return 1;
    }

    val = quadde(f2, -1.0, 1.0, 6, tol);
    printf("Test 2: %.17g\n", val);
    if (fabs(-0.69049 - val) > 1e-4) {
        return 2;
    }

    val = quadde(f3, 0.0, INFINITY, 6, tol);
    printf("Test 3: %.17g\n", val);
    if (fabs(0.21938 - val) > 1e-4) {
        return 3;
    }

    val = quadde(f4, -INFINITY, INFINITY, 6, tol);
    printf("Test 4: %.17g\n", val);
    if (fabs(2.3962 - val) > 1e-4) {
        return 4;
    }

    val = quadde(f5, -INFINITY, INFINITY, 6, tol);
    printf("Test 5: %.17g\n", val);
    if (fabs(2.2214 - val) > 1e-4) {
        return 5;
    }

    val = quadde(f6, 0.1, 1.0, 6, tol);
    printf("Test 6: %.17g\n", val);
    if (fabs(9.0 - val) > 1e-4) {
        return 6;
    }

    printf("All Tests passed!\n");

    return 0;
}

Pode-se compilar e rodar com os seguintes comandos no terminal:

gcc -o a.out quadde.c -lm && ./a.out

Com saída:

Test 1: -1.948970284582531
Test 2: -0.69049458018520249
Test 3: 0.21938393436009282
Test 4: 2.3962804694707533
Test 5: 2.2214414406307439
Test 6: 8.9999999646063316
All Tests passed!

Referências

Bailey, David H., Karthik Jeyabalan, e Xiaoye S. Li. 2005. «A Comparison of Three High-Precision Quadrature Schemes». Experimental Mathematics 14 (3): 317–29.

Engelen, Robert A. van. 2021. «Improving the Double Exponential Quadrature Tanh-Sinh, Sinh-Sinh and Exp-Sinh Formulas». Genivia, Inc. https://www.genivia.com/files/qthsh.pdf.

Takahasi, Hidetosi, e Masatake Mori. 1974. «Double Exponential Formulas for Numerical Integration». Publications of the Research Institute for Mathematical Sciences 9 (3): 721–41.

Vanherck, Joren, Bart Sorée, e Wim Magnus. 2020. «Tanh-sinh Quadrature for Single and Multiple Integration Using Floating-Point Arithmetic». arXiv. https://arxiv.org/abs/2007.15057.