diff --git a/GaussTuranExampleTemp/example.jl b/GaussTuranExampleTemp/example.jl
new file mode 100644
index 0000000..2098d52
--- /dev/null
+++ b/GaussTuranExampleTemp/example.jl
@@ -0,0 +1,61 @@
+using Optim
+using DoubleFloats
+using TaylorDiff
+using PreallocationTools
+using Integrals
+
+# This fix is needed to avoid crashing: https://github.com/JuliaLang/julia/pull/54201
+include("linear_algebra_fix.jl")
+
+FT = Double64
+n = 5
+s = 1
+a = FT(0.0)
+b = FT(1.0)
+
+"""
+    Example functions ϕⱼ from [1]
+"""
+function f(x::T, j::Int)::T where {T <: Number}
+    pow = if j % 2 == 0
+        j / 2 - 1
+    else
+        (j - 1) / 2 - 1 / 3
+    end
+    x^pow
+end
+
+I = Integrals.GaussTuran(
+    f,
+    a,
+    b,
+    n,
+    s;
+    optimization_options = Optim.Options(;
+        x_abstol = FT(1e-250),
+        g_tol = FT(1e-250),
+        show_trace = true,
+        show_every = 100,
+        iterations = 10_000
+    )
+)
+
+# Integration error |I(fⱼ) - ₐ∫ᵇfⱼ(x)dx| first N functions fⱼ
+max_int_error_upper = maximum(abs(I(x -> f(x, j)) - I.cache.rhs_upper[j])
+for j in 1:(I.cache.N))
+@show max_int_error_upper
+# max_int_error_upper = 2.465190328815662e-32
+
+# Integration error |I(fⱼ) - ₐ∫ᵇfⱼ(x)dx| last n functions fⱼ
+max_int_error_lower = maximum(
+    abs(I(x -> f(x, j)) - I.cache.rhs_lower[j - I.cache.N])
+for
+j in (I.cache.N + 1):(I.cache.N + I.cache.n)
+)
+@show max_int_error_lower
+# max_int_error_lower = 9.079315240327527e-30
+
+# Example with eˣ
+exp_error = abs(I(Base.exp) - (Base.exp(1) - 1))
+@show exp_error;
+# exp_error = 6.43662141695295e-18
diff --git a/GaussTuranExampleTemp/linear_algebra_fix.jl b/GaussTuranExampleTemp/linear_algebra_fix.jl
new file mode 100644
index 0000000..1361bb7
--- /dev/null
+++ b/GaussTuranExampleTemp/linear_algebra_fix.jl
@@ -0,0 +1,296 @@
+using LinearAlgebra
+using LinearAlgebra: AdjOrTrans, require_one_based_indexing, _ustrip
+
+# manually hoisting b[j] significantly improves performance as of Dec 2015
+# manually eliding bounds checking significantly improves performance as of Dec 2015
+# replacing repeated references to A.data[j,j] with [Ajj = A.data[j,j] and references to Ajj]
+# does not significantly impact performance as of Dec 2015
+# in the transpose and conjugate transpose naive substitution variants,
+# accumulating in z rather than b[j,k] significantly improves performance as of Dec 2015
+function LinearAlgebra.generic_trimatdiv!(
+        C::AbstractVecOrMat, uploc, isunitc, tfun::Function,
+        A::AbstractMatrix, B::AbstractVecOrMat)
+    require_one_based_indexing(C, A, B)
+    mA, nA = size(A)
+    m, n = size(B, 1), size(B, 2)
+    if nA != m
+        throw(DimensionMismatch(lazy"second dimension of left hand side A, $nA, and first dimension of right hand side B, $m, must be equal"))
+    end
+    if size(C) != size(B)
+        throw(DimensionMismatch(lazy"size of output, $(size(C)), does not match size of right hand side, $(size(B))"))
+    end
+    iszero(mA) && return C
+    oA = oneunit(eltype(A))
+    @inbounds if uploc == 'U'
+        if isunitc == 'N'
+            if tfun === identity
+                for k in 1:n
+                    amm = A[m, m]
+                    iszero(amm) && throw(SingularException(m))
+                    Cm = C[m, k] = amm \ B[m, k]
+                    # fill C-column
+                    for i in (m - 1):-1:1
+                        C[i, k] = oA \ B[i, k] - _ustrip(A[i, m]) * Cm
+                    end
+                    for j in (m - 1):-1:1
+                        ajj = A[j, j]
+                        iszero(ajj) && throw(SingularException(j))
+                        Cj = C[j, k] = _ustrip(ajj) \ C[j, k]
+                        for i in (j - 1):-1:1
+                            C[i, k] -= _ustrip(A[i, j]) * Cj
+                        end
+                    end
+                end
+            else # tfun in (adjoint, transpose)
+                for k in 1:n
+                    for j in 1:m
+                        ajj = A[j, j]
+                        iszero(ajj) && throw(SingularException(j))
+                        Bj = B[j, k]
+                        for i in 1:(j - 1)
+                            Bj -= tfun(A[i, j]) * C[i, k]
+                        end
+                        C[j, k] = tfun(ajj) \ Bj
+                    end
+                end
+            end
+        else # isunitc == 'U'
+            if tfun === identity
+                for k in 1:n
+                    Cm = C[m, k] = oA \ B[m, k]
+                    # fill C-column
+                    for i in (m - 1):-1:1
+                        C[i, k] = oA \ B[i, k] - _ustrip(A[i, m]) * Cm
+                    end
+                    for j in (m - 1):-1:1
+                        Cj = C[j, k]
+                        for i in 1:(j - 1)
+                            C[i, k] -= _ustrip(A[i, j]) * Cj
+                        end
+                    end
+                end
+            else # tfun in (adjoint, transpose)
+                for k in 1:n
+                    for j in 1:m
+                        Bj = B[j, k]
+                        for i in 1:(j - 1)
+                            Bj -= tfun(A[i, j]) * C[i, k]
+                        end
+                        C[j, k] = oA \ Bj
+                    end
+                end
+            end
+        end
+    else # uploc == 'L'
+        if isunitc == 'N'
+            if tfun === identity
+                for k in 1:n
+                    a11 = A[1, 1]
+                    iszero(a11) && throw(SingularException(1))
+                    C1 = C[1, k] = a11 \ B[1, k]
+                    # fill C-column
+                    for i in 2:m
+                        C[i, k] = oA \ B[i, k] - _ustrip(A[i, 1]) * C1
+                    end
+                    for j in 2:m
+                        ajj = A[j, j]
+                        iszero(ajj) && throw(SingularException(j))
+                        Cj = C[j, k] = _ustrip(ajj) \ C[j, k]
+                        for i in (j + 1):m
+                            C[i, k] -= _ustrip(A[i, j]) * Cj
+                        end
+                    end
+                end
+            else # tfun in (adjoint, transpose)
+                for k in 1:n
+                    for j in m:-1:1
+                        ajj = A[j, j]
+                        iszero(ajj) && throw(SingularException(j))
+                        Bj = B[j, k]
+                        for i in (j + 1):m
+                            Bj -= tfun(A[i, j]) * C[i, k]
+                        end
+                        C[j, k] = tfun(ajj) \ Bj
+                    end
+                end
+            end
+        else # isunitc == 'U'
+            if tfun === identity
+                for k in 1:n
+                    C1 = C[1, k] = oA \ B[1, k]
+                    # fill C-column
+                    for i in 2:m
+                        C[i, k] = oA \ B[i, k] - _ustrip(A[i, 1]) * C1
+                    end
+                    for j in 2:m
+                        Cj = C[j, k]
+                        for i in (j + 1):m
+                            C[i, k] -= _ustrip(A[i, j]) * Cj
+                        end
+                    end
+                end
+            else # tfun in (adjoint, transpose)
+                for k in 1:n
+                    for j in m:-1:1
+                        Bj = B[j, k]
+                        for i in (j + 1):m
+                            Bj -= tfun(A[i, j]) * C[i, k]
+                        end
+                        C[j, k] = oA \ Bj
+                    end
+                end
+            end
+        end
+    end
+    return C
+end
+# conjugate cases
+function LinearAlgebra.generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, ::Function,
+        xA::AdjOrTrans, B::AbstractVecOrMat)
+    A = parent(xA)
+    require_one_based_indexing(C, A, B)
+    mA, nA = size(A)
+    m, n = size(B, 1), size(B, 2)
+    if nA != m
+        throw(DimensionMismatch(lazy"second dimension of left hand side A, $nA, and first dimension of right hand side B, $m, must be equal"))
+    end
+    if size(C) != size(B)
+        throw(DimensionMismatch(lazy"size of output, $(size(C)), does not match size of right hand side, $(size(B))"))
+    end
+    iszero(mA) && return C
+    oA = oneunit(eltype(A))
+    @inbounds if uploc == 'U'
+        if isunitc == 'N'
+            for k in 1:n
+                amm = conj(A[m, m])
+                iszero(amm) && throw(SingularException(m))
+                Cm = C[m, k] = amm \ B[m, k]
+                # fill C-column
+                for i in (m - 1):-1:1
+                    C[i, k] = oA \ B[i, k] - _ustrip(conj(A[i, m])) * Cm
+                end
+                for j in (m - 1):-1:1
+                    ajj = conj(A[j, j])
+                    iszero(ajj) && throw(SingularException(j))
+                    Cj = C[j, k] = _ustrip(ajj) \ C[j, k]
+                    for i in (j - 1):-1:1
+                        C[i, k] -= _ustrip(conj(A[i, j])) * Cj
+                    end
+                end
+            end
+        else # isunitc == 'U'
+            for k in 1:n
+                Cm = C[m, k] = oA \ B[m, k]
+                # fill C-column
+                for i in (m - 1):-1:1
+                    C[i, k] = oA \ B[i, k] - _ustrip(conj(A[i, m])) * Cm
+                end
+                for j in (m - 1):-1:1
+                    Cj = C[j, k]
+                    for i in 1:(j - 1)
+                        C[i, k] -= _ustrip(conj(A[i, j])) * Cj
+                    end
+                end
+            end
+        end
+    else # uploc == 'L'
+        if isunitc == 'N'
+            for k in 1:n
+                a11 = conj(A[1, 1])
+                iszero(a11) && throw(SingularException(1))
+                C1 = C[1, k] = a11 \ B[1, k]
+                # fill C-column
+                for i in 2:m
+                    C[i, k] = oA \ B[i, k] - _ustrip(conj(A[i, 1])) * C1
+                end
+                for j in 2:m
+                    ajj = conj(A[j, j])
+                    iszero(ajj) && throw(SingularException(j))
+                    Cj = C[j, k] = _ustrip(ajj) \ C[j, k]
+                    for i in (j + 1):m
+                        C[i, k] -= _ustrip(conj(A[i, j])) * Cj
+                    end
+                end
+            end
+        else # isunitc == 'U'
+            for k in 1:n
+                C1 = C[1, k] = oA \ B[1, k]
+                # fill C-column
+                for i in 2:m
+                    C[i, k] = oA \ B[i, k] - _ustrip(conj(A[i, 1])) * C1
+                end
+                for j in 1:m
+                    Cj = C[j, k]
+                    for i in (j + 1):m
+                        C[i, k] -= _ustrip(conj(A[i, j])) * Cj
+                    end
+                end
+            end
+        end
+    end
+    return C
+end
+
+function generic_mattridiv!(C::AbstractMatrix, uploc, isunitc, tfun::Function,
+        A::AbstractMatrix, B::AbstractMatrix)
+    require_one_based_indexing(C, A, B)
+    m, n = size(A)
+    if size(B, 1) != n
+        throw(DimensionMismatch(lazy"right hand side B needs first dimension of size $n, has size $(size(B,1))"))
+    end
+    if size(C) != size(A)
+        throw(DimensionMismatch(lazy"size of output, $(size(C)), does not match size of left hand side, $(size(A))"))
+    end
+    oB = oneunit(eltype(B))
+    unit = isunitc == 'U'
+    @inbounds if uploc == 'U'
+        if tfun === identity
+            for i in 1:m
+                for j in 1:n
+                    Aij = A[i, j]
+                    for k in 1:(j - 1)
+                        Aij -= C[i, k] * B[k, j]
+                    end
+                    unit || (iszero(B[j, j]) && throw(SingularException(j)))
+                    C[i, j] = Aij / (unit ? oB : B[j, j])
+                end
+            end
+        else # tfun in (adjoint, transpose)
+            for i in 1:m
+                for j in n:-1:1
+                    Aij = A[i, j]
+                    for k in (j + 1):n
+                        Aij -= C[i, k] * tfun(B[j, k])
+                    end
+                    unit || (iszero(B[j, j]) && throw(SingularException(j)))
+                    C[i, j] = Aij / (unit ? oB : tfun(B[j, j]))
+                end
+            end
+        end
+    else # uploc == 'L'
+        if tfun === identity
+            for i in 1:m
+                for j in n:-1:1
+                    Aij = A[i, j]
+                    for k in (j + 1):n
+                        Aij -= C[i, k] * B[k, j]
+                    end
+                    unit || (iszero(B[j, j]) && throw(SingularException(j)))
+                    C[i, j] = Aij / (unit ? oB : B[j, j])
+                end
+            end
+        else # tfun in (adjoint, transpose)
+            for i in 1:m
+                for j in 1:n
+                    Aij = A[i, j]
+                    for k in 1:(j - 1)
+                        Aij -= C[i, k] * tfun(B[j, k])
+                    end
+                    unit || (iszero(B[j, j]) && throw(SingularException(j)))
+                    C[i, j] = Aij / (unit ? oB : tfun(B[j, j]))
+                end
+            end
+        end
+    end
+    return C
+end
diff --git a/Project.toml b/Project.toml
index d281ac3..03e0ca0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -21,6 +21,9 @@ Cubature = "667455a9-e2ce-5579-9412-b964f529a492"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 MCIntegration = "ea1e2de9-7db7-4b42-91ee-0cd1bf6df167"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
+TaylorDiff = "b36ab563-344f-407b-a36a-4f200bebf99c"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [extensions]
@@ -29,6 +32,7 @@ IntegralsCubaExt = "Cuba"
 IntegralsCubatureExt = "Cubature"
 IntegralsFastGaussQuadratureExt = "FastGaussQuadrature"
 IntegralsForwardDiffExt = "ForwardDiff"
+IntegralsGaussTuranExt = ["Optim", "TaylorDiff", "PreallocationTools"]
 IntegralsMCIntegrationExt = "MCIntegration"
 IntegralsZygoteExt = ["Zygote", "ChainRulesCore"]
 
diff --git a/ext/IntegralsGaussTuranExt.jl b/ext/IntegralsGaussTuranExt.jl
new file mode 100644
index 0000000..8b54e84
--- /dev/null
+++ b/ext/IntegralsGaussTuranExt.jl
@@ -0,0 +1,278 @@
+
+module IntegralsGaussTuranExt
+using Base.Threads
+using Integrals
+if isdefined(Base, :get_extension)
+    using Optim
+    using TaylorDiff
+    using PreallocationTools
+else
+    using ..Optim
+    using ..TaylorDiff
+    using ..PreallocationTools
+end
+
+########################################################
+## Computing Gauss-Turán quadrature rules
+
+function DEFAULT_w(x::T)::T where {T}
+    one(T)
+end
+
+"""
+Cached data for the `GaussTuranLoss!` call.
+"""
+struct GaussTuranCache{T}
+    n::Int
+    s::Int
+    N::Int
+    a::T
+    b::T
+    ε::T
+    rhs_upper::Vector{T}
+    rhs_lower::Vector{T}
+    M_upper_buffer::LazyBufferCache{typeof(identity)}
+    M_lower_buffer::LazyBufferCache{typeof(identity)}
+    X_buffer::LazyBufferCache{typeof(identity)}
+    A_buffer::LazyBufferCache{typeof(identity)}
+    function GaussTuranCache(
+            n,
+            s,
+            N,
+            a::T,
+            b::T,
+            ε::T,
+            rhs_upper::Vector{T},
+            rhs_lower::Vector{T}
+    ) where {T}
+        new{T}(
+            n,
+            s,
+            N,
+            a,
+            b,
+            ε,
+            rhs_upper,
+            rhs_lower,
+            LazyBufferCache(),
+            LazyBufferCache(),
+            LazyBufferCache(),
+            LazyBufferCache()
+        )
+    end
+end
+
+"""
+Function whose root defines the quadrature rule.
+"""
+function GaussTuranLoss!(f, ΔX::AbstractVector{T}, cache) where {T}
+    (;
+    n,
+    s,
+    N,
+    a,
+    rhs_upper,
+    rhs_lower,
+    M_upper_buffer,
+    M_lower_buffer,
+    A_buffer,
+    X_buffer
+) = cache
+    M_upper = M_upper_buffer[ΔX, (N, N)]
+    M_lower = M_lower_buffer[ΔX, (n, N)]
+    A = A_buffer[ΔX, N]
+    X = X_buffer[ΔX, n]
+
+    # Compute X from ΔX
+    cumsum!(X, ΔX)
+    X .+= a
+
+    # Evaluating f
+    for (i, x) in enumerate(X)
+        # Threads.@threads for j = 1:N
+        for j in 1:N
+            M_upper[j, i:n:N] .= derivatives(x -> f(x, j), x, 1, Val(2s + 1)).value
+        end
+        Threads.@threads for j in (N + 1):(N + n)
+            M_lower[j - N, i:n:N] .= derivatives(x -> f(x, j), x, 1, Val(2s + 1)).value
+        end
+    end
+
+    # Solving for A
+    A .= M_upper \ rhs_upper
+
+    # Computing output
+    out = zero(eltype(ΔX))
+    for i in eachindex(X)
+        out_term = -rhs_lower[i]
+        for j in eachindex(A)
+            out_term += A[j] * M_lower[i, j]
+        end
+        out += out_term^2
+    end
+    sqrt(out)
+end
+
+"""
+    Callable result object of the Gauss-Turán quadrature rule
+    computation algorithm.
+"""
+struct GaussTuranResult{T, RType, dfType}
+    X::Vector{T}
+    A::Matrix{T}
+    res::RType
+    cache::GaussTuranCache
+    df::dfType
+
+    function GaussTuranResult(res, cache::GaussTuranCache{T}, df) where {T}
+        (; A_buffer, s, n, N, a) = cache
+        X = cumsum(res.minimizer) .+ a
+        df.f(res.minimizer)
+        A = reshape(A_buffer[T[], N], (n, 2s + 1))
+        new{T, typeof(res), typeof(df)}(X, A, res, cache, df)
+    end
+end
+
+"""
+    Input: function f(x, d) which gives the dth derivative of f
+"""
+function (I::GaussTuranResult{T} where {T})(integrand)
+    (; X, A, cache) = I
+    (; s) = cache
+    out = zero(eltype(X))
+    for (i, x) in enumerate(X)
+        derivs = derivatives(integrand, x, 1.0, Val(2s + 1)).value
+        for (m, deriv) in enumerate(derivs)
+            out += A[i, m] * deriv
+        end
+    end
+    out
+end
+
+"""
+    GaussTuran(f, a, b, n, s; w = DEFAULT_w, ε = nothing, X₀ = nothing)
+
+Determine a quadrature rule
+
+I(g) = ∑ₘ∑ᵢ Aₘᵢ * ∂ᵐ⁻¹g(xᵢ)         (m = 1, … 2s + 1, i = 1, …, n)
+
+that gives the precise integral ₐ∫ᵇf(x)dx for given linearly independent functions f₁, f₂, …, f₂₍ₛ₊₁₎ₙ.
+
+Method:
+
+The equations
+
+∑ₘ∑ᵢ Aₘᵢ * ∂ᵐ⁻¹fⱼ(xᵢ) = ₐ∫ᵇw(x)fⱼ(x)dx        j = 1, …, 2(s+1)n
+
+define an overdetermined linear system M(X)A = b in the weights Aₘᵢ for a given X = (x₁, x₂, …, xₙ).
+We split the matrix M into a square upper part M_upper of size (2s+1)n x (2s+1)n and a lower part M_lower of size n x (2s+1)n,
+and the right hand size b analogously. From this we obtain A = M_upper⁻¹ * b_upper. Then we can asses the correctness of X by comparing
+M_lower * A to b_lower, i.e. how well the last n equations holds. This yields the loss function
+
+loss(X) = ||M_lower * A - b_lower||₂ = ||M_lower * M_upper⁻¹ * b_upper - b_lower||₂.
+
+We have the constraints that we want X to be ordered and in the interval (a, b). To achieve this, we formulate the loss
+in terms of ΔX = (Δx₁, Δx₂, …, Δxₙ) = (x₁ - a, x₂ - x₁, …, xₙ - xₙ₋₁) on which we set the constraints
+
+ε ≤ Δxᵢ ≤ b - a - 2 * ε     i = 1, …, n
+nε ≤ a + ∑ΔX ≤ b - a - ε
+
+where ε is an enforced minimum distance between the nodes. This prevents that consecutive nodes converge towards eachother making
+M_upper singular.
+
+## Inputs
+
+  - `f`: Function with signature `f(x::T, j)::T` that returns fⱼ at x
+  - `a`: Integration lower bound
+  - `b`: Integration upper bound
+  - `n`: The number of nodes in the quadrature rule
+  - `s`: Determines the highest order derivative required from the functions fⱼ, currently 2(s + 1)
+
+## Keyword Arguments
+
+  - `w`: the integrand weighting function, must have signature w(x::Number)::Number. Defaults to `w(x) = 1`.
+  - `ε`: the minimum distance between nodes. Defaults to 1e-3 * (b - a) / (n + 1).
+  - `X₀`: The initial guess for the nodes. Defaults to uniformly distributed over (a, b).
+  - `integration_kwargs`: The key word arguments passed to `solve` for integrating w * fⱼ
+  - `optimization_kwargs`: The key word arguments passed to `Optim.Options` for the minization problem
+    for finding X.
+"""
+function Integrals.GaussTuran(
+        f,
+        a::T,
+        b::T,
+        n,
+        s;
+        w = DEFAULT_w,
+        ε = nothing,
+        X₀ = nothing,
+        integration_kwargs::NamedTuple = (; reltol = 1e-120),
+        optimization_options::Optim.Options = Optim.Options()
+) where {T <: AbstractFloat}
+    # Initial guess
+    if isnothing(X₀)
+        X₀ = collect(range(a, b, length = n + 2)[2:(end - 1)])
+    else
+        @assert length(X₀) == n
+    end
+    ΔX₀ = diff(X₀)
+    pushfirst!(ΔX₀, X₀[1] - a)
+
+    # Minimum distance between nodes
+    if isnothing(ε)
+        ε = 1e-3 * (b - a) / (n + 1)
+    else
+        @assert 0 < ε ≤ (b - a) / (n + 1)
+    end
+    ε = T(ε)
+
+    # Integrate w * f
+    integrand = (out, x, j) -> out[] = w(x) * f(x, j)
+    function integrate(j)
+        prob = IntegralProblem{true}(integrand, (a, b), j)
+        res = solve(prob, QuadGKJL(); integration_kwargs...)
+        res.u[]
+    end
+    N = (2s + 1) * n
+    rhs_upper = [integrate(j) for j in 1:N]
+    rhs_lower = [integrate(j) for j in (N + 1):(N + n)]
+
+    # Solving constrained non linear problem for ΔX, see
+    # https://julianlsolvers.github.io/Optim.jl/stable/examples/generated/ipnewton_basics/
+
+    # The cache for evaluating GaussTuranLoss
+    cache = GaussTuranCache(n, s, N, a, b, ε, rhs_upper, rhs_lower)
+
+    # The function whose root defines the quadrature rule
+    # Note: the optimization method requires a Hessian, 
+    # which brings the highest order derivative required to 2s + 2
+    func(ΔX) = GaussTuranLoss!(f, ΔX, cache)
+    df = TwiceDifferentiable(func, ΔX₀; autodiff = :forward)
+
+    # The constraints on ΔX
+    ΔX_lb = fill(ε, length(ΔX₀))
+    ΔX_ub = fill(b - a - 2 * ε, length(ΔX₀))
+
+    # Defining the variable and constraints nε ≤ a + ∑ΔX ≤ b - a - ε
+    sum_variable!(c, ΔX) = (c[1] = a + sum(ΔX); c)
+    sum_jacobian!(J, ΔX) = (J[1, :] .= one(eltype(ΔX)); J)
+    sum_hessian!(H, ΔX, λ) = nothing
+    sum_lb = [n * ε]
+    sum_ub = [b - a - ε]
+    constraints = TwiceDifferentiableConstraints(
+        sum_variable!,
+        sum_jacobian!,
+        sum_hessian!,
+        ΔX_lb,
+        ΔX_ub,
+        sum_lb,
+        sum_ub
+    )
+
+    # Solve for the quadrature rule by minimizing the loss function
+    res = Optim.optimize(df, constraints, T.(ΔX₀), IPNewton(), optimization_options)
+
+    GaussTuranResult(res, cache, df)
+end
+
+end # module IntegralsGaussTuranExt
diff --git a/src/Integrals.jl b/src/Integrals.jl
index d11fab9..715b0f1 100644
--- a/src/Integrals.jl
+++ b/src/Integrals.jl
@@ -20,6 +20,8 @@ include("sampled.jl")
 include("trapezoidal.jl")
 include("simpsons.jl")
 
+function GaussTuran end
+
 abstract type QuadSensitivityAlg end
 struct ReCallVJP{V}
     vjp::V
diff --git a/src/algorithms_sampled.jl b/src/algorithms_sampled.jl
index bc1085f..38ad0c3 100644
--- a/src/algorithms_sampled.jl
+++ b/src/algorithms_sampled.jl
@@ -10,7 +10,7 @@ Example with sampled data:
 ```julia
 using Integrals
 f = x -> x^2
-x = range(0, 1, length=20)
+x = range(0, 1, length = 20)
 y = f.(x)
 problem = SampledIntegralProblem(y, x)
 method = TrapezoidalRule()
@@ -31,7 +31,7 @@ Example with equidistant data:
 ```julia
 using Integrals
 f = x -> x^2
-x = range(0, 1, length=20)
+x = range(0, 1, length = 20)
 y = f.(x)
 problem = SampledIntegralProblem(y, x)
 method = SimpsonsRule()