text:
# [Lith, Boonkkamp, and IJzerman](https://doi.org/10.1016/j.amc.2017.09.003)
# A family of different methods that includes the secant method and Newton's method

# return f^(i-1)(x); not the same as default eval call
function evalf(F::Callable_Function{S,T,𝑭,P}, x, i) where {N,S<:Val{N},T<:Val{true},𝑭,P}
    F.f[i](x)
end
function evalf(F::Callable_Function{S,T,𝑭,P}, x, i) where {S<:Val{1},T<:Val{false},𝑭,P}
    F.f(x)
end
evalf(::Callable_Function, x, i) = throw(
    ArgumentError(
        "LithBoonkkampIJzerman expects functions to be defined singly or as tuples",
    ),
)

"""
    LithBoonkkampIJzerman{S,D} <: AbstractNewtonLikeMethod
    LithBoonkkampIJzerman(S,D)

Specify a linear multistep solver with `S` steps and `D` derivatives following [Lith, Boonkkamp, and
IJzerman](https://doi.org/10.1016/j.amc.2017.09.003).

## Examples

```jldoctest lith
julia> using Roots

julia> find_zero(sin, 3, Roots.LithBoonkkampIJzerman(2,0)) ≈ π # the secant method
true

julia> find_zero((sin,cos), 3, Roots.LithBoonkkampIJzerman(1,1)) ≈ π # Newton's method
true

julia> find_zero((sin,cos), 3, Roots.LithBoonkkampIJzerman(3,1)) ≈ π # Faster convergence rate
true

julia> find_zero((sin,cos, x->-sin(x)), 3, Roots.LithBoonkkampIJzerman(1,2)) ≈ π # Halley-like method
true
```

The method can be more robust to the intial condition. This example is from the paper (p13). Newton's method (the `S=1`, `D=1` case) fails if `|x₀| ≥ 1.089` but methods with more memory succeed.

```jldoctest lith
julia> fx =  ZeroProblem((tanh,x->sech(x)^2), 1.239); # zero at 0.0

julia> solve(fx, Roots.LithBoonkkampIJzerman(1,1)) |> isnan# Newton, NaN
true

julia> solve(fx, Roots.LithBoonkkampIJzerman(2,1)) |> abs |> <(eps())
true

julia> solve(fx, Roots.LithBoonkkampIJzerman(3,1)) |> abs |> <(eps())
true
```

Multiple derivatives can be constructed automatically using automatic differentiation. For example,

```jldoctest lith
julia> using ForwardDiff

julia> function δ(f, n::Int=1)
           n <= 0 && return f
           n == 1 && return x -> ForwardDiff.derivative(f,float(x))
           δ(δ(f,1),n-1)
       end;

julia> fs(f,n) = ntuple(i -> δ(f,i-1), Val(n+1));

julia> f(x) = cbrt(x) * exp(-x^2); # cf. Table 6 in paper, α = 0

julia> fx = ZeroProblem(fs(f,1), 0.1147);

julia> opts = (xatol=2eps(), xrtol=0.0, atol=0.0, rtol=0.0); # converge if |xₙ - xₙ₋₁| <= 2ϵ

julia> solve(fx, Roots.LithBoonkkampIJzerman(1, 1); opts...) |> isnan # NaN -- no convergence
true

julia> solve(fx, Roots.LithBoonkkampIJzerman(2, 1); opts...) |> abs |> <(eps()) # converges
true

julia> fx = ZeroProblem(fs(f,2), 0.06);                       # need better starting point

julia> solve(fx, Roots.LithBoonkkampIJzerman(2, 2); opts...) |> abs |> <(eps()) # converges
true
```

For the case `D=1`, a bracketing method based on this approach is implemented in [`LithBoonkkampIJzermanBracket`](@ref)

## Reference

In [Lith, Boonkkamp, and
IJzerman](https://doi.org/10.1016/j.amc.2017.09.003) an analysis is
given of the convergence rates when using linear multistep methods to
solve `0=f(x)` using `f⁻¹(0) = x` when `f` is a sufficiently smooth
linear function. The reformulation, attributed to Grau-Sanchez, finds
a differential equation for `f⁻¹`: `dx/dy = [f⁻¹]′(y) = 1/f′(x) = F` as
`x(0) = x₀ + ∫⁰_y₀ F(x(y)) dy`.

A linear multi-step method is used to solve this equation
numerically.  Let S be the number of memory steps (S= 1,2,...) and D be
the number of derivatives employed, then, with `F(x) = dx/dy`
`x_{n+S} = ∑_{k=0}^{S-1} aₖ x_{n+k} +∑d=1^D ∑_{k=1}^{S-1} aᵈ_{n+k}F⁽ᵈ⁾(x_{n+k})`.
 The `aₖ`s and `aᵈₖ`s are computed each step.

This table is from Tables 1 and 3 of the paper and gives the
convergence rate for simple roots identified therein:

```
s: number of steps remembered
d: number of derivatives uses
s/d  0    1    2    3    4
1    .    2    3    4    5
2    1.62 2.73 3.79 4.82 5.85
3    1.84 2.91 3.95 4.97 5.98
4    1.92 2.97 3.99 4.99 5.996
5    1.97 .    .    .    .
```

That is, more memory leads to a higher convergence rate; more
derivatives leads to a higher convergence rate. However, the
interval about `α`, the zero, where the convergence rate is guaranteed
may get smaller.

!!! note
    For the larger values of `S`, the expressions to compute the next value get quite involved.
    The higher convergence rate is likely only to be of help for finding solutions to high precision.

"""
struct LithBoonkkampIJzerman{S,D} <: AbstractNewtonLikeMethod end
LithBoonkkampIJzerman(s::Int, d::Int) = LithBoonkkampIJzerman{s,d}()
fn_argout(::LithBoonkkampIJzerman{S,D}) where {S,D} = 1 + D

struct LithBoonkkampIJzermanState{S′,D⁺,T,S} <: AbstractUnivariateZeroState{T,S}
    xn1::T
    xn0::T
    m::NTuple{S′,T}
    fxn1::S
    fxn0::S
    fm::NTuple{D⁺,NTuple{S′,S}}
end

function init_state(L::LithBoonkkampIJzerman{S,0}, F, x) where {S}
    x₀, x₁ = x₀x₁(x)
    fx₀, fx₁ = evalf(F, x₀, 1), evalf(F, x₁, 1)
    state = init_state(L, F, x₀, x₁, fx₀, fx₁)
end

function init_state(L::LithBoonkkampIJzerman{S,D}, F, x) where {S,D}
    x₀ = float(first(x))
    fx₀ = evalf(F, x₀, 1)
    state = init_state(L, F, nan(x₀), x₀, nan(fx₀), fx₀)
end

function init_state(
    L::LithBoonkkampIJzerman{S,D},
    F,
    x₀,
    x₁::R,
    fx₀,
    fx₁::T,
) where {S,D,R,T}
    xs, ys = init_lith(L, F, x₁, fx₁, x₀, fx₀) # [x₀,x₁,…,xₛ₋₁], ...
    # skip unit consideration here, as won't fit within storage of ys
    state = LithBoonkkampIJzermanState{S,D + 1,R,T}(
        xs[end],    # xₙ
        S > 1 ? xs[end - 1] : nan(xs[end]), # xₙ₋₁
        xs,         # all xs
        ys[1][end], # fₙ
        S > 1 ? ys[1][end - 1] : nan(ys[1]), # fₙ₋₁
        ys,        #
    )

    state
end
initial_fncalls(::LithBoonkkampIJzerman{S,D}) where {S,D} = S * (D + 1)

function update_state(
    L::LithBoonkkampIJzerman{S,D},
    F::Callable_Function,
    o::LithBoonkkampIJzermanState{S⁺,D′,R,T},
    options,
    l=NullTracks(),
) where {S,D,S⁺,D′,R,T}
    xs, ys = o.m, o.fm

    xᵢ::R = lmm(L, xs, ys...)
    isissue(o.xn1 - xᵢ) && return (o, true)

    for i in 1:(S - 1)
        @set! xs[i] = xs[i + 1]
    end
    @set! xs[end] = xᵢ

    for i in 0:D
        i′ = i + 1
        for j in 1:(S - 1)
            @set! ys[i′][j] = ys[i′][j + 1]
        end
        yij::T = evalf(F, xᵢ, i′)
        @set! ys[i′][end] = yij
    end
    incfn(l, 1 + D)

    @set! o.xn0 = o.xn1
    @set! o.xn1 = xᵢ
    @set! o.fxn0 = o.fxn1
    @set! o.fxn1 = ys[1][end]
    @set! o.m = xs
    @set! o.fm = ys

    return (o, false)
end

# manufacture initial xs, ys
# use lower memory terms to boot strap up. Secant uses initial default step
#D=0, geneate [x0].x1,...,xs
function init_lith(
    L::LithBoonkkampIJzerman{S,0},
    F::Callable_Function{Si,Tup,𝑭,P},
    x₁::R,
    fx₁::T,
    x₀::R,
    fx₀::T,
) where {S,Si,Tup,𝑭,P,R,T}
    xs = NTuple{S,R}(ntuple(_ -> one(R), Val(S)))
    yᵢ = NTuple{S,T}(ntuple(_ -> one(T), Val(S)))
    ys = NTuple{1,NTuple{S,T}}((yᵢ,))

    # build up to get S of them
    x0::R = zero(R)
    if isnan(x₀)
        x0 = _default_secant_step(x₁)
        fx0::T = evalf(F, x0, 1)
    else
        x0, fx0 = x₀, fx₀
    end

    @set! xs[1] = x0
    @set! xs[2] = x₁
    @set! ys[1][1] = fx0
    @set! ys[1][2] = fx₁

    for i in 3:S
        xᵢ::R = lmm(Val(i - 1), Val(0), xs, ys) # XXX allocates due to runtime i-1
        y1i::T = evalf(F, xᵢ, 1)
        @set! xs[i] = xᵢ
        @set! ys[1][i] = y1i
    end

    xs, ys
end

#D≥1. ignore x₀
function init_lith(
    L::LithBoonkkampIJzerman{S,D},
    F::Callable_Function{Si,Tup,𝑭,P},
    x₁::R,
    fx₁::T,
    x₀::R,
    fx₀::T,
) where {S,D,Si,Tup,𝑭,P,R,T}
    xs = NTuple{S,R}(ntuple(_ -> one(R), Val(S)))
    yᵢ = NTuple{S,T}(ntuple(_ -> one(T), Val(S)))
    ys = NTuple{D + 1,NTuple{S,T}}(ntuple(_ -> yᵢ, Val(D + 1)))

    @set! xs[1] = x₁
    @set! ys[1][1] = fx₁
    for i in 1:D
        yi1::T = evalf(F, x₁, i + 1)
        @set! ys[i + 1][1] = yi1
    end

    # build up to get S of them
    for i in 2:S
        xᵢ::R = lmm(Val(i - 1), Val(D), xs, ys) # XXX allocates! clean up
        @set! xs[i] = xᵢ
        for j in 0:D
            yji::T = evalf(F, xᵢ, j + 1)
            @set! ys[j + 1][i] = yji
        end
    end

    xs, ys
end

"""
    LithBoonkkampIJzermanBracket()

A bracketing method which is a modification of Brent's method due to
[Lith, Boonkkamp, and
IJzerman](https://doi.org/10.1016/j.amc.2017.09.003). The best
possible convergence rate is 2.91.

A function, its derivative, and a bracketing interval need to be specified.

The state includes the 3 points -- a bracket `[a,b]` (`b=xₙ` has
`f(b)` closest to `0`) and `c=xₙ₋₁` -- and the corresponding values
for the function and its derivative at these three points.

The next proposed step is either a `S=2` or `S=3` selection for the
[`LithBoonkkampIJzerman`](@ref) methods with derivative information
included only if it would be of help. The proposed is modified if it
is dithering. The proposed is compared against a bisection step; the
one in the bracket and with the smaller function value is chosen as
the next step.

"""
struct LithBoonkkampIJzermanBracket <: AbstractBracketing end
struct LithBoonkkampIJzermanBracketState{T,S,R} <: AbstractUnivariateZeroState{T,S}
    xn1::T
    xn0::T
    c::T
    fxn1::S
    fxn0::S
    fc::S
    fp1::R
    fp0::R
    fpc::R
end

function init_state(L::LithBoonkkampIJzermanBracket, F, x₀, x₁, fx₀, fx₁)
    a, b, fa, fb = x₀, x₁, fx₀, fx₁
    if abs(fa) < abs(fb)
        a, b, fa, fb = b, a, fb, fa
    end

    assert_bracket(fa, fb)

    f′a, f′b = evalf(F, a, 2), evalf(F, b, 2)
    c, fc, f′c = a, fa, f′a

    # skip unit consideration here, as won't fit within storage of ys
    state = LithBoonkkampIJzermanBracketState(
        b,    # xₙ
        a, # xₙ₋₁
        c,
        fb, # fₙ
        fa, # fₙ₋₁
        fc,
        f′b,
        f′a,
        f′c,
    )

    state
end

function update_state(
    M::LithBoonkkampIJzermanBracket,
    F,
    state::LithBoonkkampIJzermanBracketState{T,S,R},
    options::UnivariateZeroOptions,
    l=NullTracks(),
) where {T,S,R}
    b::T, c::T, a::T = state.xn1, state.c, state.xn0
    fb::S, fc::S, fa::S = state.fxn1, state.fc, state.fxn0
    f′a::R, f′c::R, f′b::R = state.fp0, state.fpc, state.fp1

    # Get next interpolating step
    # decide on S and D;
    # S is 3 if a,b,c are distinct; D=1 unless all derivative info will be of the wrong sign.
    s::Int = ((a == c) || (b == c)) ? 2 : 3

    # which derivatives do we include
    sₘ = sign((fb - fa) / (b - a))
    mc, mb = sign(f′c) == sₘ, sign(f′b) == sₘ

    d₀::S = zero(S)
    if s == 2
        if mc || mb
            # D = 1
            a2s, b2s = lmm_coefficients(LithBoonkkampIJzerman{2,1}(), (c, b), (fc, fb))
            h = -fb

            d₀ = -sum(a2s .* (c, b))
            mb && (d₀ += h * b2s[2] / f′b)
            mc && (d₀ += h * b2s[1] / f′c)
        else
            d₀ = lmm(LithBoonkkampIJzerman{2,0}(), (c, b), (fc, fb))
        end
    else
        ma = sign(f′a) == sₘ
        if mc || mb || ma
            # D = 1
            as, bs = lmm_coefficients(LithBoonkkampIJzerman{3,1}(), (a, c, b), (fa, fc, fb))
            h = -fb

            d₀ = -sum(as .* (a, c, b))
            mb && (d₀ += h * bs[end] / f′b) # only when helpful
            mc && (d₀ += h * bs[end - 1] / f′c)
            ma && (d₀ += h * bs[end - 2] / f′a)
        else
            d₀ = lmm(LithBoonkkampIJzerman{3,0}(), (a, c, b), (fa, fc, fb))
        end
    end

    # If the step is smaller than the tolerance, use the tolerance as step size.
    xatol, xrtol = options.xabstol, options.xreltol
    δ = xatol + abs(b) * xrtol
    Δ₀ = b - d₀
    if abs(Δ₀) <= δ
        d₀ = b - sign(Δ₀) * δ
    end

    # compare to bisection step; extra function evalution
    d₁ = a + (b - a) * (0.5) #_middle(a, b)
    f₀, f₁ = evalf(F, d₀, 1), evalf(F, d₁, 1)

    # interpolation outside a,b or bisection better use that
    d::T, fd::S, f′d::S = zero(T), zero(S), zero(S)
    if (abs(f₀) < abs(f₁)) && (min(a, b) < d₀ < max(a, b))
        d, fd, f′d = d₀, f₀, evalf(F, d₀, 2)# interp
    else
        d, fd, f′d = d₁, f₁, evalf(F, d₁, 2)#  bisection
    end

    # either [a,d] a bracket or [d,b]
    # [a < d] < b ...c -- b -> d, c-> b (update?)
    # a < [d < b] ...c -- a -> d (update?)
    if sign(fa) * sign(fd) <= 0
        c, fc, f′c = b, fb, f′b
        b, fb, f′b = d, fd, f′d
    else
        a, fa, f′a = d, fd, f′d
    end

    # a,b bracket; keep |fb| ≤ |fa|
    if abs(fa) < abs(fb)
        c, fc, f′c = b, fb, f′b
        a, b, fa, fb, f′a, f′b = b, a, fb, fa, f′b, f′a
    end

    incfn(l, 3)
    @set! state.xn1 = b
    @set! state.xn0 = a
    @set! state.c = c
    @set! state.fxn1 = fb
    @set! state.fxn0 = fa
    @set! state.fc = fc
    @set! state.fp0 = f′a
    @set! state.fpc = f′c
    @set! state.fp1 = f′b

    return (state, false)
end

function default_tolerances(
    ::M,
    ::Type{T},
    ::Type{S},
) where {M<:LithBoonkkampIJzermanBracket,T,S}
    xatol = 2eps(T)
    xrtol = zero(one(T))
    atol = zero(float(one(S))) * oneunit(S)
    rtol = 2eps(float(one(S))) * one(S)
    maxevals = typemax(Int)
    maxfnevals = typemax(Int)
    strict = true
    (xatol, xrtol, atol, rtol, maxevals, maxfnevals, strict)
end

### ------
# Script used to generate expressions
#
# Finds expressions by assuming an interpolating polynomial
# goes through the points. From (20) in the paper
# some cases take a **long** time to run.
# At some point, early on, this gets to be more of an academic exercise than
# a performant solution

#= -------
using SymPy

# see Ansatz (20) on p10
function inverse_polynomial_interpretation(s=2,d=2)

    @vars y
    hs = [Sym("h$i") for i ∈ 0:(1+d)*s-1]
    xs = [Sym("x$i") for i ∈ 0:s-1]
    fs = [Sym("f$i") for i ∈ 0:s-1]
    f′s = [Sym("f′$i") for i ∈ 0:s-1]
    f′′s = [Sym("f′′$i") for i ∈ 0:s-1]
    f′′′s = [Sym("f′′′$i") for i ∈ 0:s-1]
    f′′′′s = [Sym("f′′′′$i") for i ∈ 0:s-1]
    f′′′′′s = [Sym("f′′′′′$i") for i ∈ 0:s-1]
    f′′′′′′s = [Sym("f′′′′′′$i") for i ∈ 0:s-1]
    h0 = first(hs)

    H(y) = sum(hs[i+1]*(y-fs[s])^i for i ∈ 0:(1+d)*s-1)
    Hⁱ = H⁰ = H(y)
    Hⁱs = eltype(H⁰)[]
    for i ∈ 1:d
        Hⁱ = diff(Hⁱ,y)
        push!(Hⁱs, Hⁱ)
    end

    eqs = Sym[subs(H(fs[i]), Dict(h0=>xs[s])) - xs[i] for i ∈ 1:s-1]
    for i ∈ 1:s
        # cf. Liptag
        f1,f2,f3,f4,f5,f6 = f′s[i],f′′s[i],f′′′s[i],f′′′′s[i],f′′′′′s[i],f′′′′′′s[i]
        g′ = 1/f1
        g′′ = -f2/f1^3
        g′′′ = (3*f2^2 - f1*f3)/(f1^5)
        g′′′′ = -(15*f2^2 + 10*f1*f2*f3+f1^2*f4)/f1^7
        g′′′′′ = (105*f2^4 -105*f1*f2^2*f3 + 10*f1^2*f3^2 + 15*f1^2*f2*f4 -f1^3*f5)/f1^9
        g′′′′′′ = (-f1^4*f6 + 21*f1^3*f2*f5 + 35*f1^3*f3*f4 - 210*f1^2*f2^2*f4 - 280*f1^2*f2*f3^2 + 1260*f1*f2^3*f3 - 945*f2^5)/f1^11
        gⁱs = [g′,g′′,g′′′, g′′′′,g′′′′′,g′′′′′′]

        for j ∈ 1:d
            push!(eqs, subs(Hⁱs[j], Dict(y=>fs[i], h0=>xs[s])) - gⁱs[j])
        end

    end

    ϕ = sympy.linsolve(eqs, hs[2:end]...)
    ϕ = first(elements(ϕ))
    ϕ = Sym.(convert(Tuple, ϕ.__pyobject__))
    D = Dict{Any,Any}(h0=>xs[s])
    for i in 1:(d+1)*s-1
        D[hs[i+1]] = ϕ[i]
    end
    subs(H(0), D) |> simplify
end

# For g = f⁻¹ return [g', g'', g''',..., g⁽ⁿ⁾]
# (cf [Liptaj](https://vixra.org/pdf/1703.0295v1.pdf)
function liptag(N)

    @vars x₀ Δₓ
    fs = [Sym("f$i") for i ∈ 1:N]
    gs = [Sym("g$i") for i ∈ 1:N]
    a(i) = fs[i]/factorial(i)
    b(i) = gs[i]/factorial(i)

    gᵏs = [1/fs[1]]
    for n ∈ 2:N
        Δy = sum(a(j) * Δₓ^j for j ∈ 1:n)
        l = x₀ + Δₓ
        r = x₀ + sum(b(i)*Δy^i for i ∈ 1:n)
        ϕ = solve(l-r, gs[n])[1]

        for j ∈ 1:n-1
            ϕ = subs(ϕ, gs[j] => gᵏs[j])
        end
        L = limit(ϕ, Δₓ => 0)
        push!(gᵏs, L)
    end
    gᵏs
end
=#

# have computed these
# S/D 0 1 2 3 4  5 6
# 1   x ✓ ✓ ✓ ✓ ✓ ✓
# 2   ✓ ✓ ✓ ✓ ✓ ✓ ✓
# 3   ✓ ✓ - - x x x
# 4   ✓ - x x x x x
# 5   ✓ x x x x x x
# 6   - x x x x x x

# - can be found with script, but answers are too long for
#   inclusion here

## We have two means to do this:
## Using coefficients as,bs, ... returned by lmm_coefficients
## x = ∑ aᵢxᵢ + ∑ⱼ₊₁ⁿ ∑ᵢ bʲᵢFʲᵢ, where Fʲ is the jth derivative of g⁻¹ (F¹ = 1/f'...)
## Using a polynomial interpolant, H(y), going through (xᵢ,fʲ(xᵢ)), j ∈ 0:N)

function lmm(::Val{S}, ::Val{D}, xs, ys) where {S,D}
    xi = ntuple(ii -> xs[ii], Val(S))
    yi = ntuple(ii -> ntuple(j -> ys[ii][j], Val(S)), Val(D + 1))
    lmm(LithBoonkkampIJzerman{S,D}(), xi, yi...)
end

# secant
function lmm(::LithBoonkkampIJzerman{2,0}, xs, fs)
    x0, x1 = xs
    f0, f1 = fs

    (f0 * x1 - f1 * x0) / (f0 - f1)
end

function lmm(::LithBoonkkampIJzerman{3,0}, xs, fs)
    x0, x1, x2 = xs
    f0, f1, f2 = fs

    (
        f0^2 * f1 * x2 - f0^2 * f2 * x1 - f0 * f1^2 * x2 + f0 * f2^2 * x1 + f1^2 * f2 * x0 -
        f1 * f2^2 * x0
    ) / (f0^2 * f1 - f0^2 * f2 - f0 * f1^2 + f0 * f2^2 + f1^2 * f2 - f1 * f2^2)
end

function lmm(::LithBoonkkampIJzerman{4,0}, xs, fs)
    x0, x1, x2, x3 = xs
    f0, f1, f2, f3 = fs

    (
        f0^3 * f1^2 * f2 * x3 - f0^3 * f1^2 * f3 * x2 - f0^3 * f1 * f2^2 * x3 +
        f0^3 * f1 * f3^2 * x2 +
        f0^3 * f2^2 * f3 * x1 - f0^3 * f2 * f3^2 * x1 - f0^2 * f1^3 * f2 * x3 +
        f0^2 * f1^3 * f3 * x2 +
        f0^2 * f1 * f2^3 * x3 - f0^2 * f1 * f3^3 * x2 - f0^2 * f2^3 * f3 * x1 +
        f0^2 * f2 * f3^3 * x1 +
        f0 * f1^3 * f2^2 * x3 - f0 * f1^3 * f3^2 * x2 - f0 * f1^2 * f2^3 * x3 +
        f0 * f1^2 * f3^3 * x2 +
        f0 * f2^3 * f3^2 * x1 - f0 * f2^2 * f3^3 * x1 - f1^3 * f2^2 * f3 * x0 +
        f1^3 * f2 * f3^2 * x0 +
        f1^2 * f2^3 * f3 * x0 - f1^2 * f2 * f3^3 * x0 - f1 * f2^3 * f3^2 * x0 +
        f1 * f2^2 * f3^3 * x0
    ) / (
        f0^3 * f1^2 * f2 - f0^3 * f1^2 * f3 - f0^3 * f1 * f2^2 +
        f0^3 * f1 * f3^2 +
        f0^3 * f2^2 * f3 - f0^3 * f2 * f3^2 - f0^2 * f1^3 * f2 +
        f0^2 * f1^3 * f3 +
        f0^2 * f1 * f2^3 - f0^2 * f1 * f3^3 - f0^2 * f2^3 * f3 +
        f0^2 * f2 * f3^3 +
        f0 * f1^3 * f2^2 - f0 * f1^3 * f3^2 - f0 * f1^2 * f2^3 +
        f0 * f1^2 * f3^3 +
        f0 * f2^3 * f3^2 - f0 * f2^2 * f3^3 - f1^3 * f2^2 * f3 +
        f1^3 * f2 * f3^2 +
        f1^2 * f2^3 * f3 - f1^2 * f2 * f3^3 - f1 * f2^3 * f3^2 + f1 * f2^2 * f3^3
    )
end

function lmm(::LithBoonkkampIJzerman{5,0}, xs, fs)
    x0, x1, x2, x3, x4 = xs
    f0, f1, f2, f3, f4 = fs

    (
        f0^4 * f1^3 * f2^2 * f3 * x4 - f0^4 * f1^3 * f2^2 * f4 * x3 -
        f0^4 * f1^3 * f2 * f3^2 * x4 +
        f0^4 * f1^3 * f2 * f4^2 * x3 +
        f0^4 * f1^3 * f3^2 * f4 * x2 - f0^4 * f1^3 * f3 * f4^2 * x2 -
        f0^4 * f1^2 * f2^3 * f3 * x4 +
        f0^4 * f1^2 * f2^3 * f4 * x3 +
        f0^4 * f1^2 * f2 * f3^3 * x4 - f0^4 * f1^2 * f2 * f4^3 * x3 -
        f0^4 * f1^2 * f3^3 * f4 * x2 +
        f0^4 * f1^2 * f3 * f4^3 * x2 +
        f0^4 * f1 * f2^3 * f3^2 * x4 - f0^4 * f1 * f2^3 * f4^2 * x3 -
        f0^4 * f1 * f2^2 * f3^3 * x4 +
        f0^4 * f1 * f2^2 * f4^3 * x3 +
        f0^4 * f1 * f3^3 * f4^2 * x2 - f0^4 * f1 * f3^2 * f4^3 * x2 -
        f0^4 * f2^3 * f3^2 * f4 * x1 +
        f0^4 * f2^3 * f3 * f4^2 * x1 +
        f0^4 * f2^2 * f3^3 * f4 * x1 - f0^4 * f2^2 * f3 * f4^3 * x1 -
        f0^4 * f2 * f3^3 * f4^2 * x1 + f0^4 * f2 * f3^2 * f4^3 * x1 -
        f0^3 * f1^4 * f2^2 * f3 * x4 +
        f0^3 * f1^4 * f2^2 * f4 * x3 +
        f0^3 * f1^4 * f2 * f3^2 * x4 - f0^3 * f1^4 * f2 * f4^2 * x3 -
        f0^3 * f1^4 * f3^2 * f4 * x2 +
        f0^3 * f1^4 * f3 * f4^2 * x2 +
        f0^3 * f1^2 * f2^4 * f3 * x4 - f0^3 * f1^2 * f2^4 * f4 * x3 -
        f0^3 * f1^2 * f2 * f3^4 * x4 +
        f0^3 * f1^2 * f2 * f4^4 * x3 +
        f0^3 * f1^2 * f3^4 * f4 * x2 - f0^3 * f1^2 * f3 * f4^4 * x2 -
        f0^3 * f1 * f2^4 * f3^2 * x4 +
        f0^3 * f1 * f2^4 * f4^2 * x3 +
        f0^3 * f1 * f2^2 * f3^4 * x4 - f0^3 * f1 * f2^2 * f4^4 * x3 -
        f0^3 * f1 * f3^4 * f4^2 * x2 +
        f0^3 * f1 * f3^2 * f4^4 * x2 +
        f0^3 * f2^4 * f3^2 * f4 * x1 - f0^3 * f2^4 * f3 * f4^2 * x1 -
        f0^3 * f2^2 * f3^4 * f4 * x1 +
        f0^3 * f2^2 * f3 * f4^4 * x1 +
        f0^3 * f2 * f3^4 * f4^2 * x1 - f0^3 * f2 * f3^2 * f4^4 * x1 +
        f0^2 * f1^4 * f2^3 * f3 * x4 - f0^2 * f1^4 * f2^3 * f4 * x3 -
        f0^2 * f1^4 * f2 * f3^3 * x4 +
        f0^2 * f1^4 * f2 * f4^3 * x3 +
        f0^2 * f1^4 * f3^3 * f4 * x2 - f0^2 * f1^4 * f3 * f4^3 * x2 -
        f0^2 * f1^3 * f2^4 * f3 * x4 +
        f0^2 * f1^3 * f2^4 * f4 * x3 +
        f0^2 * f1^3 * f2 * f3^4 * x4 - f0^2 * f1^3 * f2 * f4^4 * x3 -
        f0^2 * f1^3 * f3^4 * f4 * x2 +
        f0^2 * f1^3 * f3 * f4^4 * x2 +
        f0^2 * f1 * f2^4 * f3^3 * x4 - f0^2 * f1 * f2^4 * f4^3 * x3 -
        f0^2 * f1 * f2^3 * f3^4 * x4 +
        f0^2 * f1 * f2^3 * f4^4 * x3 +
        f0^2 * f1 * f3^4 * f4^3 * x2 - f0^2 * f1 * f3^3 * f4^4 * x2 -
        f0^2 * f2^4 * f3^3 * f4 * x1 +
        f0^2 * f2^4 * f3 * f4^3 * x1 +
        f0^2 * f2^3 * f3^4 * f4 * x1 - f0^2 * f2^3 * f3 * f4^4 * x1 -
        f0^2 * f2 * f3^4 * f4^3 * x1 + f0^2 * f2 * f3^3 * f4^4 * x1 -
        f0 * f1^4 * f2^3 * f3^2 * x4 +
        f0 * f1^4 * f2^3 * f4^2 * x3 +
        f0 * f1^4 * f2^2 * f3^3 * x4 - f0 * f1^4 * f2^2 * f4^3 * x3 -
        f0 * f1^4 * f3^3 * f4^2 * x2 +
        f0 * f1^4 * f3^2 * f4^3 * x2 +
        f0 * f1^3 * f2^4 * f3^2 * x4 - f0 * f1^3 * f2^4 * f4^2 * x3 -
        f0 * f1^3 * f2^2 * f3^4 * x4 +
        f0 * f1^3 * f2^2 * f4^4 * x3 +
        f0 * f1^3 * f3^4 * f4^2 * x2 - f0 * f1^3 * f3^2 * f4^4 * x2 -
        f0 * f1^2 * f2^4 * f3^3 * x4 +
        f0 * f1^2 * f2^4 * f4^3 * x3 +
        f0 * f1^2 * f2^3 * f3^4 * x4 - f0 * f1^2 * f2^3 * f4^4 * x3 -
        f0 * f1^2 * f3^4 * f4^3 * x2 +
        f0 * f1^2 * f3^3 * f4^4 * x2 +
        f0 * f2^4 * f3^3 * f4^2 * x1 - f0 * f2^4 * f3^2 * f4^3 * x1 -
        f0 * f2^3 * f3^4 * f4^2 * x1 +
        f0 * f2^3 * f3^2 * f4^4 * x1 +
        f0 * f2^2 * f3^4 * f4^3 * x1 - f0 * f2^2 * f3^3 * f4^4 * x1 +
        f1^4 * f2^3 * f3^2 * f4 * x0 - f1^4 * f2^3 * f3 * f4^2 * x0 -
        f1^4 * f2^2 * f3^3 * f4 * x0 +
        f1^4 * f2^2 * f3 * f4^3 * x0 +
        f1^4 * f2 * f3^3 * f4^2 * x0 - f1^4 * f2 * f3^2 * f4^3 * x0 -
        f1^3 * f2^4 * f3^2 * f4 * x0 +
        f1^3 * f2^4 * f3 * f4^2 * x0 +
        f1^3 * f2^2 * f3^4 * f4 * x0 - f1^3 * f2^2 * f3 * f4^4 * x0 -
        f1^3 * f2 * f3^4 * f4^2 * x0 +
        f1^3 * f2 * f3^2 * f4^4 * x0 +
        f1^2 * f2^4 * f3^3 * f4 * x0 - f1^2 * f2^4 * f3 * f4^3 * x0 -
        f1^2 * f2^3 * f3^4 * f4 * x0 +
        f1^2 * f2^3 * f3 * f4^4 * x0 +
        f1^2 * f2 * f3^4 * f4^3 * x0 - f1^2 * f2 * f3^3 * f4^4 * x0 -
        f1 * f2^4 * f3^3 * f4^2 * x0 +
        f1 * f2^4 * f3^2 * f4^3 * x0 +
        f1 * f2^3 * f3^4 * f4^2 * x0 - f1 * f2^3 * f3^2 * f4^4 * x0 -
        f1 * f2^2 * f3^4 * f4^3 * x0 + f1 * f2^2 * f3^3 * f4^4 * x0
    ) / (
        f0^4 * f1^3 * f2^2 * f3 - f0^4 * f1^3 * f2^2 * f4 - f0^4 * f1^3 * f2 * f3^2 +
        f0^4 * f1^3 * f2 * f4^2 +
        f0^4 * f1^3 * f3^2 * f4 - f0^4 * f1^3 * f3 * f4^2 - f0^4 * f1^2 * f2^3 * f3 +
        f0^4 * f1^2 * f2^3 * f4 +
        f0^4 * f1^2 * f2 * f3^3 - f0^4 * f1^2 * f2 * f4^3 - f0^4 * f1^2 * f3^3 * f4 +
        f0^4 * f1^2 * f3 * f4^3 +
        f0^4 * f1 * f2^3 * f3^2 - f0^4 * f1 * f2^3 * f4^2 - f0^4 * f1 * f2^2 * f3^3 +
        f0^4 * f1 * f2^2 * f4^3 +
        f0^4 * f1 * f3^3 * f4^2 - f0^4 * f1 * f3^2 * f4^3 - f0^4 * f2^3 * f3^2 * f4 +
        f0^4 * f2^3 * f3 * f4^2 +
        f0^4 * f2^2 * f3^3 * f4 - f0^4 * f2^2 * f3 * f4^3 - f0^4 * f2 * f3^3 * f4^2 +
        f0^4 * f2 * f3^2 * f4^3 - f0^3 * f1^4 * f2^2 * f3 +
        f0^3 * f1^4 * f2^2 * f4 +
        f0^3 * f1^4 * f2 * f3^2 - f0^3 * f1^4 * f2 * f4^2 - f0^3 * f1^4 * f3^2 * f4 +
        f0^3 * f1^4 * f3 * f4^2 +
        f0^3 * f1^2 * f2^4 * f3 - f0^3 * f1^2 * f2^4 * f4 - f0^3 * f1^2 * f2 * f3^4 +
        f0^3 * f1^2 * f2 * f4^4 +
        f0^3 * f1^2 * f3^4 * f4 - f0^3 * f1^2 * f3 * f4^4 - f0^3 * f1 * f2^4 * f3^2 +
        f0^3 * f1 * f2^4 * f4^2 +
        f0^3 * f1 * f2^2 * f3^4 - f0^3 * f1 * f2^2 * f4^4 - f0^3 * f1 * f3^4 * f4^2 +
        f0^3 * f1 * f3^2 * f4^4 +
        f0^3 * f2^4 * f3^2 * f4 - f0^3 * f2^4 * f3 * f4^2 - f0^3 * f2^2 * f3^4 * f4 +
        f0^3 * f2^2 * f3 * f4^4 +
        f0^3 * f2 * f3^4 * f4^2 - f0^3 * f2 * f3^2 * f4^4 + f0^2 * f1^4 * f2^3 * f3 -
        f0^2 * f1^4 * f2^3 * f4 - f0^2 * f1^4 * f2 * f3^3 +
        f0^2 * f1^4 * f2 * f4^3 +
        f0^2 * f1^4 * f3^3 * f4 - f0^2 * f1^4 * f3 * f4^3 - f0^2 * f1^3 * f2^4 * f3 +
        f0^2 * f1^3 * f2^4 * f4 +
        f0^2 * f1^3 * f2 * f3^4 - f0^2 * f1^3 * f2 * f4^4 - f0^2 * f1^3 * f3^4 * f4 +
        f0^2 * f1^3 * f3 * f4^4 +
        f0^2 * f1 * f2^4 * f3^3 - f0^2 * f1 * f2^4 * f4^3 - f0^2 * f1 * f2^3 * f3^4 +
        f0^2 * f1 * f2^3 * f4^4 +
        f0^2 * f1 * f3^4 * f4^3 - f0^2 * f1 * f3^3 * f4^4 - f0^2 * f2^4 * f3^3 * f4 +
        f0^2 * f2^4 * f3 * f4^3 +
        f0^2 * f2^3 * f3^4 * f4 - f0^2 * f2^3 * f3 * f4^4 - f0^2 * f2 * f3^4 * f4^3 +
        f0^2 * f2 * f3^3 * f4^4 - f0 * f1^4 * f2^3 * f3^2 +
        f0 * f1^4 * f2^3 * f4^2 +
        f0 * f1^4 * f2^2 * f3^3 - f0 * f1^4 * f2^2 * f4^3 - f0 * f1^4 * f3^3 * f4^2 +
        f0 * f1^4 * f3^2 * f4^3 +
        f0 * f1^3 * f2^4 * f3^2 - f0 * f1^3 * f2^4 * f4^2 - f0 * f1^3 * f2^2 * f3^4 +
        f0 * f1^3 * f2^2 * f4^4 +
        f0 * f1^3 * f3^4 * f4^2 - f0 * f1^3 * f3^2 * f4^4 - f0 * f1^2 * f2^4 * f3^3 +
        f0 * f1^2 * f2^4 * f4^3 +
        f0 * f1^2 * f2^3 * f3^4 - f0 * f1^2 * f2^3 * f4^4 - f0 * f1^2 * f3^4 * f4^3 +
        f0 * f1^2 * f3^3 * f4^4 +
        f0 * f2^4 * f3^3 * f4^2 - f0 * f2^4 * f3^2 * f4^3 - f0 * f2^3 * f3^4 * f4^2 +
        f0 * f2^3 * f3^2 * f4^4 +
        f0 * f2^2 * f3^4 * f4^3 - f0 * f2^2 * f3^3 * f4^4 + f1^4 * f2^3 * f3^2 * f4 -
        f1^4 * f2^3 * f3 * f4^2 - f1^4 * f2^2 * f3^3 * f4 +
        f1^4 * f2^2 * f3 * f4^3 +
        f1^4 * f2 * f3^3 * f4^2 - f1^4 * f2 * f3^2 * f4^3 - f1^3 * f2^4 * f3^2 * f4 +
        f1^3 * f2^4 * f3 * f4^2 +
        f1^3 * f2^2 * f3^4 * f4 - f1^3 * f2^2 * f3 * f4^4 - f1^3 * f2 * f3^4 * f4^2 +
        f1^3 * f2 * f3^2 * f4^4 +
        f1^2 * f2^4 * f3^3 * f4 - f1^2 * f2^4 * f3 * f4^3 - f1^2 * f2^3 * f3^4 * f4 +
        f1^2 * f2^3 * f3 * f4^4 +
        f1^2 * f2 * f3^4 * f4^3 - f1^2 * f2 * f3^3 * f4^4 - f1 * f2^4 * f3^3 * f4^2 +
        f1 * f2^4 * f3^2 * f4^3 +
        f1 * f2^3 * f3^4 * f4^2 - f1 * f2^3 * f3^2 * f4^4 - f1 * f2^2 * f3^4 * f4^3 +
        f1 * f2^2 * f3^3 * f4^4
    )
end

function lmm(::LithBoonkkampIJzerman{6,0}, xs, fs)
    x0, x1, x2, x3, x4, x5 = xs
    f0, f1, f2, f3, f4, f5 = fs

    error("not implemented")
end

## d = 1; Newton-like

# return (as, bs⁰,[bs¹,...,bsⁿ⁻¹])
# where -∑ aᵢ xᵢ + h ⋅ ∑ₙ (∑ bsʲᵢ Fʲᵢ)l
function lmm_coefficients(::LithBoonkkampIJzerman{1,1}, xs, fs)
    a0 = -one(xs[1])
    b0 = one(fs[1])

    return (a0,), (b0,)
end

function lmm_coefficients(::LithBoonkkampIJzerman{2,1}, xs, fs)
    q = fs[1] / fs[2]

    # from the paper
    # x2 + a1 x1 + a0x0 =  h3 * (b1 * 1/fp1 + b0 * 1/fp0)
    a0 = (1 - 3q) / (q - 1)^3
    a1 = -1 - a0
    b0 = q / (q - 1)^2
    b1 = q * b0

    return (a0, a1), (b0, b1)
end

function lmm_coefficients(::LithBoonkkampIJzerman{3,1}, xs, fs)

    # from the paper
    q0 = fs[3 - 2] / fs[3]
    q1 = fs[3 - 1] / fs[3]

    a0 = (q1^2 * (q0 * (3 + 3q1 - 5q0) - q1)) / ((q0 - 1)^3 * (q0 - q1)^3)
    a1 = (q0^2 * (q1 * (5q1 - 3q0 - 3) + q0)) / ((q1 - 1)^3 * (q0 - q1)^3)
    a2 = (q0^2 * q1^2 * (3q1 - q0 * (q1 - 3) - 5)) / ((q0 - 1)^3 * (q1 - 1)^3) # minor typo in (27c)

    b0 = (q0 * q1^2) / ((q0 - 1)^2 * (q0 - q1)^2)
    b1 = (q0^2 * q1) / ((q0 - q1)^2 * (q1 - 1)^2)
    b2 = (q0^2 * q1^2) / ((q0 - 1)^2 * (q1 - 1)^2)

    return (a0, a1, a2), (b0, b1, b2)
end

function lmm_coefficients(::LithBoonkkampIJzerman{S,1}, xs, fs) where {S}
    error("not computed")
end

function lmm(L::LithBoonkkampIJzerman{S,1}, xs, fs, f′s) where {S}
    as, bs = lmm_coefficients(L, xs, fs)
    Fs = 1 ./ f′s # F = (g⁻¹)'
    h = -fs[S]

    -sum(as[i] * xs[i] for i in 1:S) + h * sum(bs[i] * Fs[i] for i in 1:S)
end

function lmm(::LithBoonkkampIJzerman{4,1}, xs, fs, f′s)
    x0, x1, x2, x3 = xs
    f0, f1, f2, f3 = fs
    f′0, f′1, f′2, f′3 = f′s

    # can get with script, but too long as found
    error("not implemented")
end

## d = 2; Halley-like
function lmm(::LithBoonkkampIJzerman{1,2}, xs, fs, f′s, f′′s)
    x0 = xs[1]
    f0 = fs[1]
    f′0 = f′s[1]
    f′′0 = f′′s[1]

    -f0^2 * f′′0 / (2 * f′0^3) - f0 / f′0 + x0
end

function lmm(::LithBoonkkampIJzerman{2,2}, xs, fs, f′s, f′′s)
    x0, x1 = xs
    f0, f1 = fs
    f′0, f′1 = f′s
    f′′0, f′′1 = f′′s

    (
        -f0^5 * f1^2 * f′0^3 * f′′1 / 2 - f0^5 * f1 * f′0^3 * f′1^2 +
        f0^5 * f′0^3 * f′1^3 * x1 +
        f0^4 * f1^3 * f′0^3 * f′′1 +
        f0^4 * f1^3 * f′1^3 * f′′0 / 2 +
        5 * f0^4 * f1^2 * f′0^3 * f′1^2 - 5 * f0^4 * f1 * f′0^3 * f′1^3 * x1 -
        f0^3 * f1^4 * f′0^3 * f′′1 / 2 - f0^3 * f1^4 * f′1^3 * f′′0 -
        4 * f0^3 * f1^3 * f′0^3 * f′1^2 +
        4 * f0^3 * f1^3 * f′0^2 * f′1^3 +
        10 * f0^3 * f1^2 * f′0^3 * f′1^3 * x1 +
        f0^2 * f1^5 * f′1^3 * f′′0 / 2 - 5 * f0^2 * f1^4 * f′0^2 * f′1^3 -
        10 * f0^2 * f1^3 * f′0^3 * f′1^3 * x0 +
        f0 * f1^5 * f′0^2 * f′1^3 +
        5 * f0 * f1^4 * f′0^3 * f′1^3 * x0 - f1^5 * f′0^3 * f′1^3 * x0
    ) / (
        f′0^3 *
        f′1^3 *
        (f0^5 - 5 * f0^4 * f1 + 10 * f0^3 * f1^2 - 10 * f0^2 * f1^3 + 5 * f0 * f1^4 - f1^5)
    )
end

function lmm(::LithBoonkkampIJzerman{3,2}, xs, fs, f′s, f′′s)
    x0, x1, x2 = xs
    f0, f1, f2 = fs
    f′0, f′1, f′2 = f′s
    f′′0, f′′1, f′′2 = f′′s

    ## can get from script, but too long for inclusion here
    error("not implemented")
end

function lmm(::LithBoonkkampIJzerman{4,2}, xs, fs, f′s, f′′s)
    error("Not computed")
end

## d = 3
function lmm(::LithBoonkkampIJzerman{1,3}, xs, fs, f′s, f′′s, f′′′s)
    x0 = xs[1]
    f0 = fs[1]
    f′0 = f′s[1]
    f′′0 = f′′s[1]
    f′′′0 = f′′′s[1]

    (
        f0^3 * (f′0 * f′′′0 - 3 * f′′0^2) / 6 - f0^2 * f′0^2 * f′′0 / 2 - f0 * f′0^4 +
        f′0^5 * x0
    ) / f′0^5
end

function lmm(::LithBoonkkampIJzerman{2,3}, xs, fs, f′s, f′′s, f′′′s)
    x0, x1 = xs
    f0, f1 = fs
    f′0, f′1 = f′s
    f′′0, f′′1 = f′′s
    f′′′0, f′′′1 = f′′′s

    (
        f0^7 * f1^3 * f′0^5 * f′1 * f′′′1 - 3 * f0^7 * f1^3 * f′0^5 * f′′1^2 -
        3 * f0^7 * f1^2 * f′0^5 * f′1^2 * f′′1 - 6 * f0^7 * f1 * f′0^5 * f′1^4 +
        6 * f0^7 * f′0^5 * f′1^5 * x1 - 3 * f0^6 * f1^4 * f′0^5 * f′1 * f′′′1 +
        9 * f0^6 * f1^4 * f′0^5 * f′′1^2 +
        f0^6 * f1^4 * f′0 * f′1^5 * f′′′0 - 3 * f0^6 * f1^4 * f′1^5 * f′′0^2 +
        21 * f0^6 * f1^3 * f′0^5 * f′1^2 * f′′1 +
        42 * f0^6 * f1^2 * f′0^5 * f′1^4 - 42 * f0^6 * f1 * f′0^5 * f′1^5 * x1 +
        3 * f0^5 * f1^5 * f′0^5 * f′1 * f′′′1 - 9 * f0^5 * f1^5 * f′0^5 * f′′1^2 -
        3 * f0^5 * f1^5 * f′0 * f′1^5 * f′′′0 + 9 * f0^5 * f1^5 * f′1^5 * f′′0^2 -
        33 * f0^5 * f1^4 * f′0^5 * f′1^2 * f′′1 - 15 * f0^5 * f1^4 * f′0^2 * f′1^5 * f′′0 -
        126 * f0^5 * f1^3 * f′0^5 * f′1^4 + 126 * f0^5 * f1^2 * f′0^5 * f′1^5 * x1 -
        f0^4 * f1^6 * f′0^5 * f′1 * f′′′1 +
        3 * f0^4 * f1^6 * f′0^5 * f′′1^2 +
        3 * f0^4 * f1^6 * f′0 * f′1^5 * f′′′0 - 9 * f0^4 * f1^6 * f′1^5 * f′′0^2 +
        15 * f0^4 * f1^5 * f′0^5 * f′1^2 * f′′1 +
        33 * f0^4 * f1^5 * f′0^2 * f′1^5 * f′′0 +
        90 * f0^4 * f1^4 * f′0^5 * f′1^4 - 90 * f0^4 * f1^4 * f′0^4 * f′1^5 -
        210 * f0^4 * f1^3 * f′0^5 * f′1^5 * x1 - f0^3 * f1^7 * f′0 * f′1^5 * f′′′0 +
        3 * f0^3 * f1^7 * f′1^5 * f′′0^2 - 21 * f0^3 * f1^6 * f′0^2 * f′1^5 * f′′0 +
        126 * f0^3 * f1^5 * f′0^4 * f′1^5 +
        210 * f0^3 * f1^4 * f′0^5 * f′1^5 * x0 +
        3 * f0^2 * f1^7 * f′0^2 * f′1^5 * f′′0 - 42 * f0^2 * f1^6 * f′0^4 * f′1^5 -
        126 * f0^2 * f1^5 * f′0^5 * f′1^5 * x0 +
        6 * f0 * f1^7 * f′0^4 * f′1^5 +
        42 * f0 * f1^6 * f′0^5 * f′1^5 * x0 - 6 * f1^7 * f′0^5 * f′1^5 * x0
    ) / (
        6 *
        f′0^5 *
        f′1^5 *
        (
            f0^7 - 7 * f0^6 * f1 + 21 * f0^5 * f1^2 - 35 * f0^4 * f1^3 + 35 * f0^3 * f1^4 -
            21 * f0^2 * f1^5 + 7 * f0 * f1^6 - f1^7
        )
    )
end

function lmm(::LithBoonkkampIJzerman{3,3}, xs, fs, f′s, f′′s, f′′′s)
    x0, x1, x2 = xs
    f0, f1, f2 = fs
    f′0, f′1, f′2 = f′s
    f′′0, f′′1, f′′2 = f′′s
    f′′′0, f′′′1, f′′′2 = f′′′s

    # can get from script, but too long for inclusion here
    error("not implemented")
end

function lmm(::LithBoonkkampIJzerman{4,3}, xs, fs, f′s, f′′s, f′′′s)
    x0, x1, x2, x3 = xs
    f0, f1, f2, f3 = fs
    f′0, f′1, f′2, f′3 = f′s
    f′′0, f′′1, f′′2, f′′3 = f′′s
    f′′′0, f′′′1, f′′′2, f′′′3 = f′′′s

    error("not computed")
end

## d = 4
function lmm(::LithBoonkkampIJzerman{1,4}, xs, fs, f′s, f′′s, f′′′s, f′′′′s)
    x0 = xs[1]
    f0 = fs[1]
    f′0 = f′s[1]
    f′′0 = f′′s[1]
    f′′′0 = f′′′s[1]
    f′′′′0 = f′′′′s[1]

    (
        -f0^4 * (f′0^2 * f′′′′0 + 10 * f′0 * f′′0 * f′′′0 + 15 * f′′0^2) / 24 +
        f0^3 * f′0^2 * (f′0 * f′′′0 - 3 * f′′0^2) / 6 - f0^2 * f′0^4 * f′′0 / 2 -
        f0 * f′0^6 + f′0^7 * x0
    ) / f′0^7
end

function lmm(::LithBoonkkampIJzerman{2,4}, xs, fs, f′s, f′′s, f′′′s, f′′′′s)
    x0, x1 = xs
    f0, f1 = fs
    f′0, f′1 = f′s
    f′′0, f′′1 = f′′s
    f′′′0, f′′′1 = f′′′s
    f′′′′0, f′′′′1 = f′′′′s

    (
        -f0^9 * f1^4 * f′0^7 * f′1^2 * f′′′′1 -
        10 * f0^9 * f1^4 * f′0^7 * f′1 * f′′1 * f′′′1 - 15 * f0^9 * f1^4 * f′0^7 * f′′1^2 +
        4 * f0^9 * f1^3 * f′0^7 * f′1^3 * f′′′1 -
        12 * f0^9 * f1^3 * f′0^7 * f′1^2 * f′′1^2 -
        12 * f0^9 * f1^2 * f′0^7 * f′1^4 * f′′1 - 24 * f0^9 * f1 * f′0^7 * f′1^6 +
        24 * f0^9 * f′0^7 * f′1^7 * x1 +
        4 * f0^8 * f1^5 * f′0^7 * f′1^2 * f′′′′1 +
        40 * f0^8 * f1^5 * f′0^7 * f′1 * f′′1 * f′′′1 +
        60 * f0^8 * f1^5 * f′0^7 * f′′1^2 +
        f0^8 * f1^5 * f′0^2 * f′1^7 * f′′′′0 +
        10 * f0^8 * f1^5 * f′0 * f′1^7 * f′′0 * f′′′0 +
        15 * f0^8 * f1^5 * f′1^7 * f′′0^2 - 36 * f0^8 * f1^4 * f′0^7 * f′1^3 * f′′′1 +
        108 * f0^8 * f1^4 * f′0^7 * f′1^2 * f′′1^2 +
        108 * f0^8 * f1^3 * f′0^7 * f′1^4 * f′′1 +
        216 * f0^8 * f1^2 * f′0^7 * f′1^6 - 216 * f0^8 * f1 * f′0^7 * f′1^7 * x1 -
        6 * f0^7 * f1^6 * f′0^7 * f′1^2 * f′′′′1 -
        60 * f0^7 * f1^6 * f′0^7 * f′1 * f′′1 * f′′′1 - 90 * f0^7 * f1^6 * f′0^7 * f′′1^2 -
        4 * f0^7 * f1^6 * f′0^2 * f′1^7 * f′′′′0 -
        40 * f0^7 * f1^6 * f′0 * f′1^7 * f′′0 * f′′′0 - 60 * f0^7 * f1^6 * f′1^7 * f′′0^2 +
        84 * f0^7 * f1^5 * f′0^7 * f′1^3 * f′′′1 -
        252 * f0^7 * f1^5 * f′0^7 * f′1^2 * f′′1^2 -
        24 * f0^7 * f1^5 * f′0^3 * f′1^7 * f′′′0 +
        72 * f0^7 * f1^5 * f′0^2 * f′1^7 * f′′0^2 -
        432 * f0^7 * f1^4 * f′0^7 * f′1^4 * f′′1 - 864 * f0^7 * f1^3 * f′0^7 * f′1^6 +
        864 * f0^7 * f1^2 * f′0^7 * f′1^7 * x1 +
        4 * f0^6 * f1^7 * f′0^7 * f′1^2 * f′′′′1 +
        40 * f0^6 * f1^7 * f′0^7 * f′1 * f′′1 * f′′′1 +
        60 * f0^6 * f1^7 * f′0^7 * f′′1^2 +
        6 * f0^6 * f1^7 * f′0^2 * f′1^7 * f′′′′0 +
        60 * f0^6 * f1^7 * f′0 * f′1^7 * f′′0 * f′′′0 +
        90 * f0^6 * f1^7 * f′1^7 * f′′0^2 - 76 * f0^6 * f1^6 * f′0^7 * f′1^3 * f′′′1 +
        228 * f0^6 * f1^6 * f′0^7 * f′1^2 * f′′1^2 +
        76 * f0^6 * f1^6 * f′0^3 * f′1^7 * f′′′0 -
        228 * f0^6 * f1^6 * f′0^2 * f′1^7 * f′′0^2 +
        588 * f0^6 * f1^5 * f′0^7 * f′1^4 * f′′1 +
        252 * f0^6 * f1^5 * f′0^4 * f′1^7 * f′′0 +
        2016 * f0^6 * f1^4 * f′0^7 * f′1^6 - 2016 * f0^6 * f1^3 * f′0^7 * f′1^7 * x1 -
        f0^5 * f1^8 * f′0^7 * f′1^2 * f′′′′1 -
        10 * f0^5 * f1^8 * f′0^7 * f′1 * f′′1 * f′′′1 - 15 * f0^5 * f1^8 * f′0^7 * f′′1^2 -
        4 * f0^5 * f1^8 * f′0^2 * f′1^7 * f′′′′0 -
        40 * f0^5 * f1^8 * f′0 * f′1^7 * f′′0 * f′′′0 - 60 * f0^5 * f1^8 * f′1^7 * f′′0^2 +
        24 * f0^5 * f1^7 * f′0^7 * f′1^3 * f′′′1 -
        72 * f0^5 * f1^7 * f′0^7 * f′1^2 * f′′1^2 -
        84 * f0^5 * f1^7 * f′0^3 * f′1^7 * f′′′0 +
        252 * f0^5 * f1^7 * f′0^2 * f′1^7 * f′′0^2 -
        252 * f0^5 * f1^6 * f′0^7 * f′1^4 * f′′1 -
        588 * f0^5 * f1^6 * f′0^4 * f′1^7 * f′′0 - 1344 * f0^5 * f1^5 * f′0^7 * f′1^6 +
        1344 * f0^5 * f1^5 * f′0^6 * f′1^7 +
        3024 * f0^5 * f1^4 * f′0^7 * f′1^7 * x1 +
        f0^4 * f1^9 * f′0^2 * f′1^7 * f′′′′0 +
        10 * f0^4 * f1^9 * f′0 * f′1^7 * f′′0 * f′′′0 +
        15 * f0^4 * f1^9 * f′1^7 * f′′0^2 +
        36 * f0^4 * f1^8 * f′0^3 * f′1^7 * f′′′0 -
        108 * f0^4 * f1^8 * f′0^2 * f′1^7 * f′′0^2 +
        432 * f0^4 * f1^7 * f′0^4 * f′1^7 * f′′0 - 2016 * f0^4 * f1^6 * f′0^6 * f′1^7 -
        3024 * f0^4 * f1^5 * f′0^7 * f′1^7 * x0 - 4 * f0^3 * f1^9 * f′0^3 * f′1^7 * f′′′0 +
        12 * f0^3 * f1^9 * f′0^2 * f′1^7 * f′′0^2 -
        108 * f0^3 * f1^8 * f′0^4 * f′1^7 * f′′0 +
        864 * f0^3 * f1^7 * f′0^6 * f′1^7 +
        2016 * f0^3 * f1^6 * f′0^7 * f′1^7 * x0 +
        12 * f0^2 * f1^9 * f′0^4 * f′1^7 * f′′0 - 216 * f0^2 * f1^8 * f′0^6 * f′1^7 -
        864 * f0^2 * f1^7 * f′0^7 * f′1^7 * x0 +
        24 * f0 * f1^9 * f′0^6 * f′1^7 +
        216 * f0 * f1^8 * f′0^7 * f′1^7 * x0 - 24 * f1^9 * f′0^7 * f′1^7 * x0
    ) / (
        24 *
        f′0^7 *
        f′1^7 *
        (
            f0^9 - 9 * f0^8 * f1 + 36 * f0^7 * f1^2 - 84 * f0^6 * f1^3 + 126 * f0^5 * f1^4 -
            126 * f0^4 * f1^5 + 84 * f0^3 * f1^6 - 36 * f0^2 * f1^7 + 9 * f0 * f1^8 - f1^9
        )
    )
end

function lmm(::LithBoonkkampIJzerman{3,4}, xs, fs, f′s, f′′s, f′′′s, f′′′′s)
    x0, x1, x2 = xs
    f0, f1, f2 = fs
    f′0, f′1, f′2 = f′s
    f′′0, f′′1, f′′2 = f′′s
    f′′′0, f′′′1, f′′′2 = f′′′s
    f′′′′0, f′′′′1, f′′′′2 = f′′′′s

    error("not computed")
end

# n = 5

function lmm(::LithBoonkkampIJzerman{1,5}, xs, fs, f′s, f′′s, f′′′s, f′′′′s, f′′′′′s)
    x0 = xs[1]
    f0 = fs[1]
    f′0 = f′s[1]
    f′′0 = f′′s[1]
    f′′′0 = f′′′s[1]
    f′′′′0 = f′′′′s[1]
    f′′′′′0 = f′′′′′s[1]

    (
        f0^5 * (
            f′0^3 * f′′′′′0 - 15 * f′0^2 * f′′0 * f′′′′0 - 10 * f′0^2 * f′′′0^2 +
            105 * f′0 * f′′0^2 * f′′′0 - 105 * f′′0^4
        ) - 5 * f0^4 * f′0^2 * (f′0^2 * f′′′′0 + 10 * f′0 * f′′0 * f′′′0 + 15 * f′′0^2) +
        20 * f0^3 * f′0^4 * (f′0 * f′′′0 - 3 * f′′0^2) - 60 * f0^2 * f′0^6 * f′′0 -
        120 * f0 * f′0^8 + 120 * f′0^9 * x0
    ) / (120 * f′0^9)
end

function lmm(::LithBoonkkampIJzerman{2,5}, xs, fs, f′s, f′′s, f′′′s, f′′′′s, f′′′′′s)
    x0, x1 = xs
    f0, f1 = fs
    f′0, f′1 = f′s
    f′′0, f′′1 = f′′s
    f′′′0, f′′′1 = f′′′s
    f′′′′0, f′′′′1 = f′′′′s
    f′′′′′0, f′′′′′1 = f′′′′′s

    (
        f0^11 * f1^5 * f′0^9 * f′1^3 * f′′′′′1 -
        15 * f0^11 * f1^5 * f′0^9 * f′1^2 * f′′1 * f′′′′1 -
        10 * f0^11 * f1^5 * f′0^9 * f′1^2 * f′′′1^2 +
        105 * f0^11 * f1^5 * f′0^9 * f′1 * f′′1^2 * f′′′1 -
        105 * f0^11 * f1^5 * f′0^9 * f′′1^4 - 5 * f0^11 * f1^4 * f′0^9 * f′1^4 * f′′′′1 -
        50 * f0^11 * f1^4 * f′0^9 * f′1^3 * f′′1 * f′′′1 -
        75 * f0^11 * f1^4 * f′0^9 * f′1^2 * f′′1^2 +
        20 * f0^11 * f1^3 * f′0^9 * f′1^5 * f′′′1 -
        60 * f0^11 * f1^3 * f′0^9 * f′1^4 * f′′1^2 -
        60 * f0^11 * f1^2 * f′0^9 * f′1^6 * f′′1 - 120 * f0^11 * f1 * f′0^9 * f′1^8 +
        120 * f0^11 * f′0^9 * f′1^9 * x1 - 5 * f0^10 * f1^6 * f′0^9 * f′1^3 * f′′′′′1 +
        75 * f0^10 * f1^6 * f′0^9 * f′1^2 * f′′1 * f′′′′1 +
        50 * f0^10 * f1^6 * f′0^9 * f′1^2 * f′′′1^2 -
        525 * f0^10 * f1^6 * f′0^9 * f′1 * f′′1^2 * f′′′1 +
        525 * f0^10 * f1^6 * f′0^9 * f′′1^4 +
        f0^10 * f1^6 * f′0^3 * f′1^9 * f′′′′′0 -
        15 * f0^10 * f1^6 * f′0^2 * f′1^9 * f′′0 * f′′′′0 -
        10 * f0^10 * f1^6 * f′0^2 * f′1^9 * f′′′0^2 +
        105 * f0^10 * f1^6 * f′0 * f′1^9 * f′′0^2 * f′′′0 -
        105 * f0^10 * f1^6 * f′1^9 * f′′0^4 +
        55 * f0^10 * f1^5 * f′0^9 * f′1^4 * f′′′′1 +
        550 * f0^10 * f1^5 * f′0^9 * f′1^3 * f′′1 * f′′′1 +
        825 * f0^10 * f1^5 * f′0^9 * f′1^2 * f′′1^2 -
        220 * f0^10 * f1^4 * f′0^9 * f′1^5 * f′′′1 +
        660 * f0^10 * f1^4 * f′0^9 * f′1^4 * f′′1^2 +
        660 * f0^10 * f1^3 * f′0^9 * f′1^6 * f′′1 +
        1320 * f0^10 * f1^2 * f′0^9 * f′1^8 - 1320 * f0^10 * f1 * f′0^9 * f′1^9 * x1 +
        10 * f0^9 * f1^7 * f′0^9 * f′1^3 * f′′′′′1 -
        150 * f0^9 * f1^7 * f′0^9 * f′1^2 * f′′1 * f′′′′1 -
        100 * f0^9 * f1^7 * f′0^9 * f′1^2 * f′′′1^2 +
        1050 * f0^9 * f1^7 * f′0^9 * f′1 * f′′1^2 * f′′′1 -
        1050 * f0^9 * f1^7 * f′0^9 * f′′1^4 - 5 * f0^9 * f1^7 * f′0^3 * f′1^9 * f′′′′′0 +
        75 * f0^9 * f1^7 * f′0^2 * f′1^9 * f′′0 * f′′′′0 +
        50 * f0^9 * f1^7 * f′0^2 * f′1^9 * f′′′0^2 -
        525 * f0^9 * f1^7 * f′0 * f′1^9 * f′′0^2 * f′′′0 +
        525 * f0^9 * f1^7 * f′1^9 * f′′0^4 - 170 * f0^9 * f1^6 * f′0^9 * f′1^4 * f′′′′1 -
        1700 * f0^9 * f1^6 * f′0^9 * f′1^3 * f′′1 * f′′′1 -
        2550 * f0^9 * f1^6 * f′0^9 * f′1^2 * f′′1^2 -
        35 * f0^9 * f1^6 * f′0^4 * f′1^9 * f′′′′0 -
        350 * f0^9 * f1^6 * f′0^3 * f′1^9 * f′′0 * f′′′0 -
        525 * f0^9 * f1^6 * f′0^2 * f′1^9 * f′′0^2 +
        1100 * f0^9 * f1^5 * f′0^9 * f′1^5 * f′′′1 -
        3300 * f0^9 * f1^5 * f′0^9 * f′1^4 * f′′1^2 -
        3300 * f0^9 * f1^4 * f′0^9 * f′1^6 * f′′1 - 6600 * f0^9 * f1^3 * f′0^9 * f′1^8 +
        6600 * f0^9 * f1^2 * f′0^9 * f′1^9 * x1 -
        10 * f0^8 * f1^8 * f′0^9 * f′1^3 * f′′′′′1 +
        150 * f0^8 * f1^8 * f′0^9 * f′1^2 * f′′1 * f′′′′1 +
        100 * f0^8 * f1^8 * f′0^9 * f′1^2 * f′′′1^2 -
        1050 * f0^8 * f1^8 * f′0^9 * f′1 * f′′1^2 * f′′′1 +
        1050 * f0^8 * f1^8 * f′0^9 * f′′1^4 +
        10 * f0^8 * f1^8 * f′0^3 * f′1^9 * f′′′′′0 -
        150 * f0^8 * f1^8 * f′0^2 * f′1^9 * f′′0 * f′′′′0 -
        100 * f0^8 * f1^8 * f′0^2 * f′1^9 * f′′′0^2 +
        1050 * f0^8 * f1^8 * f′0 * f′1^9 * f′′0^2 * f′′′0 -
        1050 * f0^8 * f1^8 * f′1^9 * f′′0^4 +
        230 * f0^8 * f1^7 * f′0^9 * f′1^4 * f′′′′1 +
        2300 * f0^8 * f1^7 * f′0^9 * f′1^3 * f′′1 * f′′′1 +
        3450 * f0^8 * f1^7 * f′0^9 * f′1^2 * f′′1^2 +
        145 * f0^8 * f1^7 * f′0^4 * f′1^9 * f′′′′0 +
        1450 * f0^8 * f1^7 * f′0^3 * f′1^9 * f′′0 * f′′′0 +
        2175 * f0^8 * f1^7 * f′0^2 * f′1^9 * f′′0^2 -
        2180 * f0^8 * f1^6 * f′0^9 * f′1^5 * f′′′1 +
        6540 * f0^8 * f1^6 * f′0^9 * f′1^4 * f′′1^2 +
        560 * f0^8 * f1^6 * f′0^5 * f′1^9 * f′′′0 -
        1680 * f0^8 * f1^6 * f′0^4 * f′1^9 * f′′0^2 +
        9900 * f0^8 * f1^5 * f′0^9 * f′1^6 * f′′1 +
        19800 * f0^8 * f1^4 * f′0^9 * f′1^8 - 19800 * f0^8 * f1^3 * f′0^9 * f′1^9 * x1 +
        5 * f0^7 * f1^9 * f′0^9 * f′1^3 * f′′′′′1 -
        75 * f0^7 * f1^9 * f′0^9 * f′1^2 * f′′1 * f′′′′1 -
        50 * f0^7 * f1^9 * f′0^9 * f′1^2 * f′′′1^2 +
        525 * f0^7 * f1^9 * f′0^9 * f′1 * f′′1^2 * f′′′1 -
        525 * f0^7 * f1^9 * f′0^9 * f′′1^4 - 10 * f0^7 * f1^9 * f′0^3 * f′1^9 * f′′′′′0 +
        150 * f0^7 * f1^9 * f′0^2 * f′1^9 * f′′0 * f′′′′0 +
        100 * f0^7 * f1^9 * f′0^2 * f′1^9 * f′′′0^2 -
        1050 * f0^7 * f1^9 * f′0 * f′1^9 * f′′0^2 * f′′′0 +
        1050 * f0^7 * f1^9 * f′1^9 * f′′0^4 - 145 * f0^7 * f1^8 * f′0^9 * f′1^4 * f′′′′1 -
        1450 * f0^7 * f1^8 * f′0^9 * f′1^3 * f′′1 * f′′′1 -
        2175 * f0^7 * f1^8 * f′0^9 * f′1^2 * f′′1^2 -
        230 * f0^7 * f1^8 * f′0^4 * f′1^9 * f′′′′0 -
        2300 * f0^7 * f1^8 * f′0^3 * f′1^9 * f′′0 * f′′′0 -
        3450 * f0^7 * f1^8 * f′0^2 * f′1^9 * f′′0^2 +
        1840 * f0^7 * f1^7 * f′0^9 * f′1^5 * f′′′1 -
        5520 * f0^7 * f1^7 * f′0^9 * f′1^4 * f′′1^2 -
        1840 * f0^7 * f1^7 * f′0^5 * f′1^9 * f′′′0 +
        5520 * f0^7 * f1^7 * f′0^4 * f′1^9 * f′′0^2 -
        12240 * f0^7 * f1^6 * f′0^9 * f′1^6 * f′′1 -
        5040 * f0^7 * f1^6 * f′0^6 * f′1^9 * f′′0 - 39600 * f0^7 * f1^5 * f′0^9 * f′1^8 +
        39600 * f0^7 * f1^4 * f′0^9 * f′1^9 * x1 - f0^6 * f1^10 * f′0^9 * f′1^3 * f′′′′′1 +
        15 * f0^6 * f1^10 * f′0^9 * f′1^2 * f′′1 * f′′′′1 +
        10 * f0^6 * f1^10 * f′0^9 * f′1^2 * f′′′1^2 -
        105 * f0^6 * f1^10 * f′0^9 * f′1 * f′′1^2 * f′′′1 +
        105 * f0^6 * f1^10 * f′0^9 * f′′1^4 +
        5 * f0^6 * f1^10 * f′0^3 * f′1^9 * f′′′′′0 -
        75 * f0^6 * f1^10 * f′0^2 * f′1^9 * f′′0 * f′′′′0 -
        50 * f0^6 * f1^10 * f′0^2 * f′1^9 * f′′′0^2 +
        525 * f0^6 * f1^10 * f′0 * f′1^9 * f′′0^2 * f′′′0 -
        525 * f0^6 * f1^10 * f′1^9 * f′′0^4 +
        35 * f0^6 * f1^9 * f′0^9 * f′1^4 * f′′′′1 +
        350 * f0^6 * f1^9 * f′0^9 * f′1^3 * f′′1 * f′′′1 +
        525 * f0^6 * f1^9 * f′0^9 * f′1^2 * f′′1^2 +
        170 * f0^6 * f1^9 * f′0^4 * f′1^9 * f′′′′0 +
        1700 * f0^6 * f1^9 * f′0^3 * f′1^9 * f′′0 * f′′′0 +
        2550 * f0^6 * f1^9 * f′0^2 * f′1^9 * f′′0^2 -
        560 * f0^6 * f1^8 * f′0^9 * f′1^5 * f′′′1 +
        1680 * f0^6 * f1^8 * f′0^9 * f′1^4 * f′′1^2 +
        2180 * f0^6 * f1^8 * f′0^5 * f′1^9 * f′′′0 -
        6540 * f0^6 * f1^8 * f′0^4 * f′1^9 * f′′0^2 +
        5040 * f0^6 * f1^7 * f′0^9 * f′1^6 * f′′1 +
        12240 * f0^6 * f1^7 * f′0^6 * f′1^9 * f′′0 +
        25200 * f0^6 * f1^6 * f′0^9 * f′1^8 - 25200 * f0^6 * f1^6 * f′0^8 * f′1^9 -
        55440 * f0^6 * f1^5 * f′0^9 * f′1^9 * x1 - f0^5 * f1^11 * f′0^3 * f′1^9 * f′′′′′0 +
        15 * f0^5 * f1^11 * f′0^2 * f′1^9 * f′′0 * f′′′′0 +
        10 * f0^5 * f1^11 * f′0^2 * f′1^9 * f′′′0^2 -
        105 * f0^5 * f1^11 * f′0 * f′1^9 * f′′0^2 * f′′′0 +
        105 * f0^5 * f1^11 * f′1^9 * f′′0^4 - 55 * f0^5 * f1^10 * f′0^4 * f′1^9 * f′′′′0 -
        550 * f0^5 * f1^10 * f′0^3 * f′1^9 * f′′0 * f′′′0 -
        825 * f0^5 * f1^10 * f′0^2 * f′1^9 * f′′0^2 -
        1100 * f0^5 * f1^9 * f′0^5 * f′1^9 * f′′′0 +
        3300 * f0^5 * f1^9 * f′0^4 * f′1^9 * f′′0^2 -
        9900 * f0^5 * f1^8 * f′0^6 * f′1^9 * f′′0 +
        39600 * f0^5 * f1^7 * f′0^8 * f′1^9 +
        55440 * f0^5 * f1^6 * f′0^9 * f′1^9 * x0 +
        5 * f0^4 * f1^11 * f′0^4 * f′1^9 * f′′′′0 +
        50 * f0^4 * f1^11 * f′0^3 * f′1^9 * f′′0 * f′′′0 +
        75 * f0^4 * f1^11 * f′0^2 * f′1^9 * f′′0^2 +
        220 * f0^4 * f1^10 * f′0^5 * f′1^9 * f′′′0 -
        660 * f0^4 * f1^10 * f′0^4 * f′1^9 * f′′0^2 +
        3300 * f0^4 * f1^9 * f′0^6 * f′1^9 * f′′0 - 19800 * f0^4 * f1^8 * f′0^8 * f′1^9 -
        39600 * f0^4 * f1^7 * f′0^9 * f′1^9 * x0 -
        20 * f0^3 * f1^11 * f′0^5 * f′1^9 * f′′′0 +
        60 * f0^3 * f1^11 * f′0^4 * f′1^9 * f′′0^2 -
        660 * f0^3 * f1^10 * f′0^6 * f′1^9 * f′′0 +
        6600 * f0^3 * f1^9 * f′0^8 * f′1^9 +
        19800 * f0^3 * f1^8 * f′0^9 * f′1^9 * x0 +
        60 * f0^2 * f1^11 * f′0^6 * f′1^9 * f′′0 - 1320 * f0^2 * f1^10 * f′0^8 * f′1^9 -
        6600 * f0^2 * f1^9 * f′0^9 * f′1^9 * x0 +
        120 * f0 * f1^11 * f′0^8 * f′1^9 +
        1320 * f0 * f1^10 * f′0^9 * f′1^9 * x0 - 120 * f1^11 * f′0^9 * f′1^9 * x0
    ) / (
        120 *
        f′0^9 *
        f′1^9 *
        (
            f0^11 - 11 * f0^10 * f1 + 55 * f0^9 * f1^2 - 165 * f0^8 * f1^3 +
            330 * f0^7 * f1^4 - 462 * f0^6 * f1^5 + 462 * f0^5 * f1^6 - 330 * f0^4 * f1^7 +
            165 * f0^3 * f1^8 - 55 * f0^2 * f1^9 + 11 * f0 * f1^10 - f1^11
        )
    )
end

function lmm(::LithBoonkkampIJzerman{3,5}, xs, fs, f′s, f′′s, f′′′s, f′′′′s, f′′′′′s)
    x0, x1, x2 = xs
    f0, f1, f2 = fs
    f′0, f′1, f′2 = f′s
    f′′0, f′′1, f′′2 = f′′s
    f′′′0, f′′′1, f′′′2 = f′′′s
    f′′′′0, f′′′′1, f′′′′2 = f′′′′s
    f′′′′′0, f′′′′′1, f′′′′′2 = f′′′′′s

    error("not computed")
end

## n = 6
function lmm(
    ::LithBoonkkampIJzerman{1,6},
    xs,
    fs,
    f′s,
    f′′s,
    f′′′s,
    f′′′′s,
    f′′′′′s,
    f′′′′′′s,
)
    x0 = xs[1]
    f0 = fs[1]
    f′0 = f′s[1]
    f′′0 = f′′s[1]
    f′′′0 = f′′′s[1]
    f′′′′0 = f′′′′s[1]
    f′′′′′0 = f′′′′′s[1]
    f′′′′′′0 = f′′′′′′s[1]

    (
        f0^6 * (
            -f′0^4 * f′′′′′′0 + 21 * f′0^3 * f′′0 * f′′′′′0 + 35 * f′0^3 * f′′′0 * f′′′′0 -
            210 * f′0^2 * f′′0^2 * f′′′′0 - 280 * f′0^2 * f′′0 * f′′′0^2 +
            1260 * f′0 * f′′0^3 * f′′′0 - 945 * f′′0^5
        ) +
        6 *
        f0^5 *
        f′0^2 *
        (
            f′0^3 * f′′′′′0 - 15 * f′0^2 * f′′0 * f′′′′0 - 10 * f′0^2 * f′′′0^2 +
            105 * f′0 * f′′0^2 * f′′′0 - 105 * f′′0^4
        ) - 30 * f0^4 * f′0^4 * (f′0^2 * f′′′′0 + 10 * f′0 * f′′0 * f′′′0 + 15 * f′′0^2) +
        120 * f0^3 * f′0^6 * (f′0 * f′′′0 - 3 * f′′0^2) - 360 * f0^2 * f′0^8 * f′′0 -
        720 * f0 * f′0^10 + 720 * f′0^11 * x0
    ) / (720 * f′0^11)
end

function lmm(
    ::LithBoonkkampIJzerman{2,6},
    xs,
    fs,
    f′s,
    f′′s,
    f′′′s,
    f′′′′s,
    f′′′′′s,
    f′′′′′′s,
)
    x0, x1 = xs
    f0, f1 = fs
    f′0, f′1 = f′s
    f′′0, f′′1 = f′′s
    f′′′0, f′′′1 = f′′′s
    f′′′′0, f′′′′1 = f′′′′s
    f′′′′′0, f′′′′′1 = f′′′′′s
    f′′′′′′0, f′′′′′′1 = f′′′′′′s

    (
        -f0^13 * f1^6 * f′0^11 * f′1^4 * f′′′′′′1 +
        21 * f0^13 * f1^6 * f′0^11 * f′1^3 * f′′1 * f′′′′′1 +
        35 * f0^13 * f1^6 * f′0^11 * f′1^3 * f′′′1 * f′′′′1 -
        210 * f0^13 * f1^6 * f′0^11 * f′1^2 * f′′1^2 * f′′′′1 -
        280 * f0^13 * f1^6 * f′0^11 * f′1^2 * f′′1 * f′′′1^2 +
        1260 * f0^13 * f1^6 * f′0^11 * f′1 * f′′1^3 * f′′′1 -
        945 * f0^13 * f1^6 * f′0^11 * f′′1^5 + 6 * f0^13 * f1^5 * f′0^11 * f′1^5 * f′′′′′1 -
        90 * f0^13 * f1^5 * f′0^11 * f′1^4 * f′′1 * f′′′′1 -
        60 * f0^13 * f1^5 * f′0^11 * f′1^4 * f′′′1^2 +
        630 * f0^13 * f1^5 * f′0^11 * f′1^3 * f′′1^2 * f′′′1 -
        630 * f0^13 * f1^5 * f′0^11 * f′1^2 * f′′1^4 -
        30 * f0^13 * f1^4 * f′0^11 * f′1^6 * f′′′′1 -
        300 * f0^13 * f1^4 * f′0^11 * f′1^5 * f′′1 * f′′′1 -
        450 * f0^13 * f1^4 * f′0^11 * f′1^4 * f′′1^2 +
        120 * f0^13 * f1^3 * f′0^11 * f′1^7 * f′′′1 -
        360 * f0^13 * f1^3 * f′0^11 * f′1^6 * f′′1^2 -
        360 * f0^13 * f1^2 * f′0^11 * f′1^8 * f′′1 - 720 * f0^13 * f1 * f′0^11 * f′1^10 +
        720 * f0^13 * f′0^11 * f′1^11 * x1 +
        6 * f0^12 * f1^7 * f′0^11 * f′1^4 * f′′′′′′1 -
        126 * f0^12 * f1^7 * f′0^11 * f′1^3 * f′′1 * f′′′′′1 -
        210 * f0^12 * f1^7 * f′0^11 * f′1^3 * f′′′1 * f′′′′1 +
        1260 * f0^12 * f1^7 * f′0^11 * f′1^2 * f′′1^2 * f′′′′1 +
        1680 * f0^12 * f1^7 * f′0^11 * f′1^2 * f′′1 * f′′′1^2 -
        7560 * f0^12 * f1^7 * f′0^11 * f′1 * f′′1^3 * f′′′1 +
        5670 * f0^12 * f1^7 * f′0^11 * f′′1^5 +
        f0^12 * f1^7 * f′0^4 * f′1^11 * f′′′′′′0 -
        21 * f0^12 * f1^7 * f′0^3 * f′1^11 * f′′0 * f′′′′′0 -
        35 * f0^12 * f1^7 * f′0^3 * f′1^11 * f′′′0 * f′′′′0 +
        210 * f0^12 * f1^7 * f′0^2 * f′1^11 * f′′0^2 * f′′′′0 +
        280 * f0^12 * f1^7 * f′0^2 * f′1^11 * f′′0 * f′′′0^2 -
        1260 * f0^12 * f1^7 * f′0 * f′1^11 * f′′0^3 * f′′′0 +
        945 * f0^12 * f1^7 * f′1^11 * f′′0^5 -
        78 * f0^12 * f1^6 * f′0^11 * f′1^5 * f′′′′′1 +
        1170 * f0^12 * f1^6 * f′0^11 * f′1^4 * f′′1 * f′′′′1 +
        780 * f0^12 * f1^6 * f′0^11 * f′1^4 * f′′′1^2 -
        8190 * f0^12 * f1^6 * f′0^11 * f′1^3 * f′′1^2 * f′′′1 +
        8190 * f0^12 * f1^6 * f′0^11 * f′1^2 * f′′1^4 +
        390 * f0^12 * f1^5 * f′0^11 * f′1^6 * f′′′′1 +
        3900 * f0^12 * f1^5 * f′0^11 * f′1^5 * f′′1 * f′′′1 +
        5850 * f0^12 * f1^5 * f′0^11 * f′1^4 * f′′1^2 -
        1560 * f0^12 * f1^4 * f′0^11 * f′1^7 * f′′′1 +
        4680 * f0^12 * f1^4 * f′0^11 * f′1^6 * f′′1^2 +
        4680 * f0^12 * f1^3 * f′0^11 * f′1^8 * f′′1 +
        9360 * f0^12 * f1^2 * f′0^11 * f′1^10 - 9360 * f0^12 * f1 * f′0^11 * f′1^11 * x1 -
        15 * f0^11 * f1^8 * f′0^11 * f′1^4 * f′′′′′′1 +
        315 * f0^11 * f1^8 * f′0^11 * f′1^3 * f′′1 * f′′′′′1 +
        525 * f0^11 * f1^8 * f′0^11 * f′1^3 * f′′′1 * f′′′′1 -
        3150 * f0^11 * f1^8 * f′0^11 * f′1^2 * f′′1^2 * f′′′′1 -
        4200 * f0^11 * f1^8 * f′0^11 * f′1^2 * f′′1 * f′′′1^2 +
        18900 * f0^11 * f1^8 * f′0^11 * f′1 * f′′1^3 * f′′′1 -
        14175 * f0^11 * f1^8 * f′0^11 * f′′1^5 -
        6 * f0^11 * f1^8 * f′0^4 * f′1^11 * f′′′′′′0 +
        126 * f0^11 * f1^8 * f′0^3 * f′1^11 * f′′0 * f′′′′′0 +
        210 * f0^11 * f1^8 * f′0^3 * f′1^11 * f′′′0 * f′′′′0 -
        1260 * f0^11 * f1^8 * f′0^2 * f′1^11 * f′′0^2 * f′′′′0 -
        1680 * f0^11 * f1^8 * f′0^2 * f′1^11 * f′′0 * f′′′0^2 +
        7560 * f0^11 * f1^8 * f′0 * f′1^11 * f′′0^3 * f′′′0 -
        5670 * f0^11 * f1^8 * f′1^11 * f′′0^5 +
        300 * f0^11 * f1^7 * f′0^11 * f′1^5 * f′′′′′1 -
        4500 * f0^11 * f1^7 * f′0^11 * f′1^4 * f′′1 * f′′′′1 -
        3000 * f0^11 * f1^7 * f′0^11 * f′1^4 * f′′′1^2 +
        31500 * f0^11 * f1^7 * f′0^11 * f′1^3 * f′′1^2 * f′′′1 -
        31500 * f0^11 * f1^7 * f′0^11 * f′1^2 * f′′1^4 -
        48 * f0^11 * f1^7 * f′0^5 * f′1^11 * f′′′′′0 +
        720 * f0^11 * f1^7 * f′0^4 * f′1^11 * f′′0 * f′′′′0 +
        480 * f0^11 * f1^7 * f′0^4 * f′1^11 * f′′′0^2 -
        5040 * f0^11 * f1^7 * f′0^3 * f′1^11 * f′′0^2 * f′′′0 +
        5040 * f0^11 * f1^7 * f′0^2 * f′1^11 * f′′0^4 -
        2340 * f0^11 * f1^6 * f′0^11 * f′1^6 * f′′′′1 -
        23400 * f0^11 * f1^6 * f′0^11 * f′1^5 * f′′1 * f′′′1 -
        35100 * f0^11 * f1^6 * f′0^11 * f′1^4 * f′′1^2 +
        9360 * f0^11 * f1^5 * f′0^11 * f′1^7 * f′′′1 -
        28080 * f0^11 * f1^5 * f′0^11 * f′1^6 * f′′1^2 -
        28080 * f0^11 * f1^4 * f′0^11 * f′1^8 * f′′1 -
        56160 * f0^11 * f1^3 * f′0^11 * f′1^10 +
        56160 * f0^11 * f1^2 * f′0^11 * f′1^11 * x1 +
        20 * f0^10 * f1^9 * f′0^11 * f′1^4 * f′′′′′′1 -
        420 * f0^10 * f1^9 * f′0^11 * f′1^3 * f′′1 * f′′′′′1 -
        700 * f0^10 * f1^9 * f′0^11 * f′1^3 * f′′′1 * f′′′′1 +
        4200 * f0^10 * f1^9 * f′0^11 * f′1^2 * f′′1^2 * f′′′′1 +
        5600 * f0^10 * f1^9 * f′0^11 * f′1^2 * f′′1 * f′′′1^2 -
        25200 * f0^10 * f1^9 * f′0^11 * f′1 * f′′1^3 * f′′′1 +
        18900 * f0^10 * f1^9 * f′0^11 * f′′1^5 +
        15 * f0^10 * f1^9 * f′0^4 * f′1^11 * f′′′′′′0 -
        315 * f0^10 * f1^9 * f′0^3 * f′1^11 * f′′0 * f′′′′′0 -
        525 * f0^10 * f1^9 * f′0^3 * f′1^11 * f′′′0 * f′′′′0 +
        3150 * f0^10 * f1^9 * f′0^2 * f′1^11 * f′′0^2 * f′′′′0 +
        4200 * f0^10 * f1^9 * f′0^2 * f′1^11 * f′′0 * f′′′0^2 -
        18900 * f0^10 * f1^9 * f′0 * f′1^11 * f′′0^3 * f′′′0 +
        14175 * f0^10 * f1^9 * f′1^11 * f′′0^5 -
        540 * f0^10 * f1^8 * f′0^11 * f′1^5 * f′′′′′1 +
        8100 * f0^10 * f1^8 * f′0^11 * f′1^4 * f′′1 * f′′′′1 +
        5400 * f0^10 * f1^8 * f′0^11 * f′1^4 * f′′′1^2 -
        56700 * f0^10 * f1^8 * f′0^11 * f′1^3 * f′′1^2 * f′′′1 +
        56700 * f0^10 * f1^8 * f′0^11 * f′1^2 * f′′1^4 +
        246 * f0^10 * f1^8 * f′0^5 * f′1^11 * f′′′′′0 -
        3690 * f0^10 * f1^8 * f′0^4 * f′1^11 * f′′0 * f′′′′0 -
        2460 * f0^10 * f1^8 * f′0^4 * f′1^11 * f′′′0^2 +
        25830 * f0^10 * f1^8 * f′0^3 * f′1^11 * f′′0^2 * f′′′0 -
        25830 * f0^10 * f1^8 * f′0^2 * f′1^11 * f′′0^4 +
        6060 * f0^10 * f1^7 * f′0^11 * f′1^6 * f′′′′1 +
        60600 * f0^10 * f1^7 * f′0^11 * f′1^5 * f′′1 * f′′′1 +
        90900 * f0^10 * f1^7 * f′0^11 * f′1^4 * f′′1^2 +
        1080 * f0^10 * f1^7 * f′0^6 * f′1^11 * f′′′′0 +
        10800 * f0^10 * f1^7 * f′0^5 * f′1^11 * f′′0 * f′′′0 +
        16200 * f0^10 * f1^7 * f′0^4 * f′1^11 * f′′0^2 -
        34320 * f0^10 * f1^6 * f′0^11 * f′1^7 * f′′′1 +
        102960 * f0^10 * f1^6 * f′0^11 * f′1^6 * f′′1^2 +
        102960 * f0^10 * f1^5 * f′0^11 * f′1^8 * f′′1 +
        205920 * f0^10 * f1^4 * f′0^11 * f′1^10 -
        205920 * f0^10 * f1^3 * f′0^11 * f′1^11 * x1 -
        15 * f0^9 * f1^10 * f′0^11 * f′1^4 * f′′′′′′1 +
        315 * f0^9 * f1^10 * f′0^11 * f′1^3 * f′′1 * f′′′′′1 +
        525 * f0^9 * f1^10 * f′0^11 * f′1^3 * f′′′1 * f′′′′1 -
        3150 * f0^9 * f1^10 * f′0^11 * f′1^2 * f′′1^2 * f′′′′1 -
        4200 * f0^9 * f1^10 * f′0^11 * f′1^2 * f′′1 * f′′′1^2 +
        18900 * f0^9 * f1^10 * f′0^11 * f′1 * f′′1^3 * f′′′1 -
        14175 * f0^9 * f1^10 * f′0^11 * f′′1^5 -
        20 * f0^9 * f1^10 * f′0^4 * f′1^11 * f′′′′′′0 +
        420 * f0^9 * f1^10 * f′0^3 * f′1^11 * f′′0 * f′′′′′0 +
        700 * f0^9 * f1^10 * f′0^3 * f′1^11 * f′′′0 * f′′′′0 -
        4200 * f0^9 * f1^10 * f′0^2 * f′1^11 * f′′0^2 * f′′′′0 -
        5600 * f0^9 * f1^10 * f′0^2 * f′1^11 * f′′0 * f′′′0^2 +
        25200 * f0^9 * f1^10 * f′0 * f′1^11 * f′′0^3 * f′′′0 -
        18900 * f0^9 * f1^10 * f′1^11 * f′′0^5 +
        510 * f0^9 * f1^9 * f′0^11 * f′1^5 * f′′′′′1 -
        7650 * f0^9 * f1^9 * f′0^11 * f′1^4 * f′′1 * f′′′′1 -
        5100 * f0^9 * f1^9 * f′0^11 * f′1^4 * f′′′1^2 +
        53550 * f0^9 * f1^9 * f′0^11 * f′1^3 * f′′1^2 * f′′′1 -
        53550 * f0^9 * f1^9 * f′0^11 * f′1^2 * f′′1^4 -
        510 * f0^9 * f1^9 * f′0^5 * f′1^11 * f′′′′′0 +
        7650 * f0^9 * f1^9 * f′0^4 * f′1^11 * f′′0 * f′′′′0 +
        5100 * f0^9 * f1^9 * f′0^4 * f′1^11 * f′′′0^2 -
        53550 * f0^9 * f1^9 * f′0^3 * f′1^11 * f′′0^2 * f′′′0 +
        53550 * f0^9 * f1^9 * f′0^2 * f′1^11 * f′′0^4 -
        7590 * f0^9 * f1^8 * f′0^11 * f′1^6 * f′′′′1 -
        75900 * f0^9 * f1^8 * f′0^11 * f′1^5 * f′′1 * f′′′1 -
        113850 * f0^9 * f1^8 * f′0^11 * f′1^4 * f′′1^2 -
        4590 * f0^9 * f1^8 * f′0^6 * f′1^11 * f′′′′0 -
        45900 * f0^9 * f1^8 * f′0^5 * f′1^11 * f′′0 * f′′′0 -
        68850 * f0^9 * f1^8 * f′0^4 * f′1^11 * f′′0^2 +
        60600 * f0^9 * f1^7 * f′0^11 * f′1^7 * f′′′1 -
        181800 * f0^9 * f1^7 * f′0^11 * f′1^6 * f′′1^2 -
        14400 * f0^9 * f1^7 * f′0^7 * f′1^11 * f′′′0 +
        43200 * f0^9 * f1^7 * f′0^6 * f′1^11 * f′′0^2 -
        257400 * f0^9 * f1^6 * f′0^11 * f′1^8 * f′′1 -
        514800 * f0^9 * f1^5 * f′0^11 * f′1^10 +
        514800 * f0^9 * f1^4 * f′0^11 * f′1^11 * x1 +
        6 * f0^8 * f1^11 * f′0^11 * f′1^4 * f′′′′′′1 -
        126 * f0^8 * f1^11 * f′0^11 * f′1^3 * f′′1 * f′′′′′1 -
        210 * f0^8 * f1^11 * f′0^11 * f′1^3 * f′′′1 * f′′′′1 +
        1260 * f0^8 * f1^11 * f′0^11 * f′1^2 * f′′1^2 * f′′′′1 +
        1680 * f0^8 * f1^11 * f′0^11 * f′1^2 * f′′1 * f′′′1^2 -
        7560 * f0^8 * f1^11 * f′0^11 * f′1 * f′′1^3 * f′′′1 +
        5670 * f0^8 * f1^11 * f′0^11 * f′′1^5 +
        15 * f0^8 * f1^11 * f′0^4 * f′1^11 * f′′′′′′0 -
        315 * f0^8 * f1^11 * f′0^3 * f′1^11 * f′′0 * f′′′′′0 -
        525 * f0^8 * f1^11 * f′0^3 * f′1^11 * f′′′0 * f′′′′0 +
        3150 * f0^8 * f1^11 * f′0^2 * f′1^11 * f′′0^2 * f′′′′0 +
        4200 * f0^8 * f1^11 * f′0^2 * f′1^11 * f′′0 * f′′′0^2 -
        18900 * f0^8 * f1^11 * f′0 * f′1^11 * f′′0^3 * f′′′0 +
        14175 * f0^8 * f1^11 * f′1^11 * f′′0^5 -
        246 * f0^8 * f1^10 * f′0^11 * f′1^5 * f′′′′′1 +
        3690 * f0^8 * f1^10 * f′0^11 * f′1^4 * f′′1 * f′′′′1 +
        2460 * f0^8 * f1^10 * f′0^11 * f′1^4 * f′′′1^2 -
        25830 * f0^8 * f1^10 * f′0^11 * f′1^3 * f′′1^2 * f′′′1 +
        25830 * f0^8 * f1^10 * f′0^11 * f′1^2 * f′′1^4 +
        540 * f0^8 * f1^10 * f′0^5 * f′1^11 * f′′′′′0 -
        8100 * f0^8 * f1^10 * f′0^4 * f′1^11 * f′′0 * f′′′′0 -
        5400 * f0^8 * f1^10 * f′0^4 * f′1^11 * f′′′0^2 +
        56700 * f0^8 * f1^10 * f′0^3 * f′1^11 * f′′0^2 * f′′′0 -
        56700 * f0^8 * f1^10 * f′0^2 * f′1^11 * f′′0^4 +
        4590 * f0^8 * f1^9 * f′0^11 * f′1^6 * f′′′′1 +
        45900 * f0^8 * f1^9 * f′0^11 * f′1^5 * f′′1 * f′′′1 +
        68850 * f0^8 * f1^9 * f′0^11 * f′1^4 * f′′1^2 +
        7590 * f0^8 * f1^9 * f′0^6 * f′1^11 * f′′′′0 +
        75900 * f0^8 * f1^9 * f′0^5 * f′1^11 * f′′0 * f′′′0 +
        113850 * f0^8 * f1^9 * f′0^4 * f′1^11 * f′′0^2 -
        48600 * f0^8 * f1^8 * f′0^11 * f′1^7 * f′′′1 +
        145800 * f0^8 * f1^8 * f′0^11 * f′1^6 * f′′1^2 +
        48600 * f0^8 * f1^8 * f′0^7 * f′1^11 * f′′′0 -
        145800 * f0^8 * f1^8 * f′0^6 * f′1^11 * f′′0^2 +
        297000 * f0^8 * f1^7 * f′0^11 * f′1^8 * f′′1 +
        118800 * f0^8 * f1^7 * f′0^8 * f′1^11 * f′′0 +
        926640 * f0^8 * f1^6 * f′0^11 * f′1^10 -
        926640 * f0^8 * f1^5 * f′0^11 * f′1^11 * x1 -
        f0^7 * f1^12 * f′0^11 * f′1^4 * f′′′′′′1 +
        21 * f0^7 * f1^12 * f′0^11 * f′1^3 * f′′1 * f′′′′′1 +
        35 * f0^7 * f1^12 * f′0^11 * f′1^3 * f′′′1 * f′′′′1 -
        210 * f0^7 * f1^12 * f′0^11 * f′1^2 * f′′1^2 * f′′′′1 -
        280 * f0^7 * f1^12 * f′0^11 * f′1^2 * f′′1 * f′′′1^2 +
        1260 * f0^7 * f1^12 * f′0^11 * f′1 * f′′1^3 * f′′′1 -
        945 * f0^7 * f1^12 * f′0^11 * f′′1^5 -
        6 * f0^7 * f1^12 * f′0^4 * f′1^11 * f′′′′′′0 +
        126 * f0^7 * f1^12 * f′0^3 * f′1^11 * f′′0 * f′′′′′0 +
        210 * f0^7 * f1^12 * f′0^3 * f′1^11 * f′′′0 * f′′′′0 -
        1260 * f0^7 * f1^12 * f′0^2 * f′1^11 * f′′0^2 * f′′′′0 -
        1680 * f0^7 * f1^12 * f′0^2 * f′1^11 * f′′0 * f′′′0^2 +
        7560 * f0^7 * f1^12 * f′0 * f′1^11 * f′′0^3 * f′′′0 -
        5670 * f0^7 * f1^12 * f′1^11 * f′′0^5 +
        48 * f0^7 * f1^11 * f′0^11 * f′1^5 * f′′′′′1 -
        720 * f0^7 * f1^11 * f′0^11 * f′1^4 * f′′1 * f′′′′1 -
        480 * f0^7 * f1^11 * f′0^11 * f′1^4 * f′′′1^2 +
        5040 * f0^7 * f1^11 * f′0^11 * f′1^3 * f′′1^2 * f′′′1 -
        5040 * f0^7 * f1^11 * f′0^11 * f′1^2 * f′′1^4 -
        300 * f0^7 * f1^11 * f′0^5 * f′1^11 * f′′′′′0 +
        4500 * f0^7 * f1^11 * f′0^4 * f′1^11 * f′′0 * f′′′′0 +
        3000 * f0^7 * f1^11 * f′0^4 * f′1^11 * f′′′0^2 -
        31500 * f0^7 * f1^11 * f′0^3 * f′1^11 * f′′0^2 * f′′′0 +
        31500 * f0^7 * f1^11 * f′0^2 * f′1^11 * f′′0^4 -
        1080 * f0^7 * f1^10 * f′0^11 * f′1^6 * f′′′′1 -
        10800 * f0^7 * f1^10 * f′0^11 * f′1^5 * f′′1 * f′′′1 -
        16200 * f0^7 * f1^10 * f′0^11 * f′1^4 * f′′1^2 -
        6060 * f0^7 * f1^10 * f′0^6 * f′1^11 * f′′′′0 -
        60600 * f0^7 * f1^10 * f′0^5 * f′1^11 * f′′0 * f′′′0 -
        90900 * f0^7 * f1^10 * f′0^4 * f′1^11 * f′′0^2 +
        14400 * f0^7 * f1^9 * f′0^11 * f′1^7 * f′′′1 -
        43200 * f0^7 * f1^9 * f′0^11 * f′1^6 * f′′1^2 -
        60600 * f0^7 * f1^9 * f′0^7 * f′1^11 * f′′′0 +
        181800 * f0^7 * f1^9 * f′0^6 * f′1^11 * f′′0^2 -
        118800 * f0^7 * f1^8 * f′0^11 * f′1^8 * f′′1 -
        297000 * f0^7 * f1^8 * f′0^8 * f′1^11 * f′′0 -
        570240 * f0^7 * f1^7 * f′0^11 * f′1^10 +
        570240 * f0^7 * f1^7 * f′0^10 * f′1^11 +
        1235520 * f0^7 * f1^6 * f′0^11 * f′1^11 * x1 +
        f0^6 * f1^13 * f′0^4 * f′1^11 * f′′′′′′0 -
        21 * f0^6 * f1^13 * f′0^3 * f′1^11 * f′′0 * f′′′′′0 -
        35 * f0^6 * f1^13 * f′0^3 * f′1^11 * f′′′0 * f′′′′0 +
        210 * f0^6 * f1^13 * f′0^2 * f′1^11 * f′′0^2 * f′′′′0 +
        280 * f0^6 * f1^13 * f′0^2 * f′1^11 * f′′0 * f′′′0^2 -
        1260 * f0^6 * f1^13 * f′0 * f′1^11 * f′′0^3 * f′′′0 +
        945 * f0^6 * f1^13 * f′1^11 * f′′0^5 +
        78 * f0^6 * f1^12 * f′0^5 * f′1^11 * f′′′′′0 -
        1170 * f0^6 * f1^12 * f′0^4 * f′1^11 * f′′0 * f′′′′0 -
        780 * f0^6 * f1^12 * f′0^4 * f′1^11 * f′′′0^2 +
        8190 * f0^6 * f1^12 * f′0^3 * f′1^11 * f′′0^2 * f′′′0 -
        8190 * f0^6 * f1^12 * f′0^2 * f′1^11 * f′′0^4 +
        2340 * f0^6 * f1^11 * f′0^6 * f′1^11 * f′′′′0 +
        23400 * f0^6 * f1^11 * f′0^5 * f′1^11 * f′′0 * f′′′0 +
        35100 * f0^6 * f1^11 * f′0^4 * f′1^11 * f′′0^2 +
        34320 * f0^6 * f1^10 * f′0^7 * f′1^11 * f′′′0 -
        102960 * f0^6 * f1^10 * f′0^6 * f′1^11 * f′′0^2 +
        257400 * f0^6 * f1^9 * f′0^8 * f′1^11 * f′′0 -
        926640 * f0^6 * f1^8 * f′0^10 * f′1^11 -
        1235520 * f0^6 * f1^7 * f′0^11 * f′1^11 * x0 -
        6 * f0^5 * f1^13 * f′0^5 * f′1^11 * f′′′′′0 +
        90 * f0^5 * f1^13 * f′0^4 * f′1^11 * f′′0 * f′′′′0 +
        60 * f0^5 * f1^13 * f′0^4 * f′1^11 * f′′′0^2 -
        630 * f0^5 * f1^13 * f′0^3 * f′1^11 * f′′0^2 * f′′′0 +
        630 * f0^5 * f1^13 * f′0^2 * f′1^11 * f′′0^4 -
        390 * f0^5 * f1^12 * f′0^6 * f′1^11 * f′′′′0 -
        3900 * f0^5 * f1^12 * f′0^5 * f′1^11 * f′′0 * f′′′0 -
        5850 * f0^5 * f1^12 * f′0^4 * f′1^11 * f′′0^2 -
        9360 * f0^5 * f1^11 * f′0^7 * f′1^11 * f′′′0 +
        28080 * f0^5 * f1^11 * f′0^6 * f′1^11 * f′′0^2 -
        102960 * f0^5 * f1^10 * f′0^8 * f′1^11 * f′′0 +
        514800 * f0^5 * f1^9 * f′0^10 * f′1^11 +
        926640 * f0^5 * f1^8 * f′0^11 * f′1^11 * x0 +
        30 * f0^4 * f1^13 * f′0^6 * f′1^11 * f′′′′0 +
        300 * f0^4 * f1^13 * f′0^5 * f′1^11 * f′′0 * f′′′0 +
        450 * f0^4 * f1^13 * f′0^4 * f′1^11 * f′′0^2 +
        1560 * f0^4 * f1^12 * f′0^7 * f′1^11 * f′′′0 -
        4680 * f0^4 * f1^12 * f′0^6 * f′1^11 * f′′0^2 +
        28080 * f0^4 * f1^11 * f′0^8 * f′1^11 * f′′0 -
        205920 * f0^4 * f1^10 * f′0^10 * f′1^11 -
        514800 * f0^4 * f1^9 * f′0^11 * f′1^11 * x0 -
        120 * f0^3 * f1^13 * f′0^7 * f′1^11 * f′′′0 +
        360 * f0^3 * f1^13 * f′0^6 * f′1^11 * f′′0^2 -
        4680 * f0^3 * f1^12 * f′0^8 * f′1^11 * f′′0 +
        56160 * f0^3 * f1^11 * f′0^10 * f′1^11 +
        205920 * f0^3 * f1^10 * f′0^11 * f′1^11 * x0 +
        360 * f0^2 * f1^13 * f′0^8 * f′1^11 * f′′0 - 9360 * f0^2 * f1^12 * f′0^10 * f′1^11 -
        56160 * f0^2 * f1^11 * f′0^11 * f′1^11 * x0 +
        720 * f0 * f1^13 * f′0^10 * f′1^11 +
        9360 * f0 * f1^12 * f′0^11 * f′1^11 * x0 - 720 * f1^13 * f′0^11 * f′1^11 * x0
    ) / (
        720 *
        f′0^11 *
        f′1^11 *
        (
            f0^13 - 13 * f0^12 * f1 + 78 * f0^11 * f1^2 - 286 * f0^10 * f1^3 +
            715 * f0^9 * f1^4 - 1287 * f0^8 * f1^5 + 1716 * f0^7 * f1^6 -
            1716 * f0^6 * f1^7 + 1287 * f0^5 * f1^8 - 715 * f0^4 * f1^9 +
            286 * f0^3 * f1^10 - 78 * f0^2 * f1^11 + 13 * f0 * f1^12 - f1^13
        )
    )
end

function lmm(
    ::LithBoonkkampIJzerman{3,6},
    xs,
    fs,
    f′s,
    f′′s,
    f′′′s,
    f′′′′s,
    f′′′′′s,
    f′′′′′′s,
)
    error("not computed")
end


text:
#############################
# Algorithms from Section 3 #
#############################
export MultivariateVandermondeMatrix, FindEquations

"""
FindEquations(data::Array{T,2},
              alg::Symbol,
              d::Int64,
              homogeneous_equations::Bool)
              where  {T<:Number}

Finds equations.
* data is a matrix whose colums are the data points in Ω.
* alg is the algorithm that should be used (one of :with_svd, :with_qr, :with_rref).
* d is the degree of the equations.
* homogeneous_equations = true restricts the search space to homogeneous polynomials.
* homogeneous_equations = false computes all polynomials of degree at most d.

Alternatively, you can specifiy the exponents of the monomials.

FindEquations(data::Array{T,2},
              alg::Symbol,
              exponents::Array{Array{Int64,1},1})
              where  {T<:Number}

Here, exponents is an array of exponents.
"""
function FindEquations(data::Array{T,2}, alg::Symbol, exponents::Array{Array{Int64,1},1}) where  {T<:Number}
    M = MultivariateVandermondeMatrix(data, exponents)
    FindEquations(M,alg)
end
function FindEquations(data::Array{T,2}, alg::Symbol, degree::Int64, homogeneous_equations::Bool) where  {T<:Number}
    @assert degree > 0 "The degree must be positive."
    M = MultivariateVandermondeMatrix(data, degree, homogeneous_equations)
    FindEquations(M, alg)
end

"""
MultivariateVandermondeMatrix(data::Array{T},
                              d::Int64,
                              homogeneous_equations::Bool)

Creates a multivariate Vandermonde.
* data is a matrix whose colums are the data points in Ω.
* d is the degree of the monomials.
* homogeneous_equations = true restricts the space of monomials to monomials of degree d.
* homogeneous_equations = false computes all monomials of degree at most d.

Alternatively, you can specifiy the exponents of the monomials.

FindEquations(data::Array{T,2},
              alg::Symbol,
              exponents::Array{Array{Int64,1},1})
              where  {T<:Number}

Here, exponents is an array of exponents.
"""
# MultivariateVandermondeMatrix struct
struct MultivariateVandermondeMatrix
    Vandermonde::Array
    exponents::Vector

    function MultivariateVandermondeMatrix(data::Array{T}, exponents::Vector) where {T<:Number}

        @assert length(unique(length.(exponents))) == 1 "Error: Exponents differ in size."

        n,m = size(data)
        N = length(exponents)
        v = veronese(exponents,T, n)
        U = zeros(T,m,N)
        for i=1:m
            U[i,:] = v(data[:,i])
        end
        new(U, exponents)
    end

    function MultivariateVandermondeMatrix(data::Array{T}, d::Int64,  homogeneous_equations::Bool) where {T<:Number}
        n=size(data,1)
        if homogeneous_equations
            exponents = collect(Combinatorics.multiexponents(n, d))
        else
            exponents = vcat(map(i -> collect(Combinatorics.multiexponents(n,-i)), -d:0)...)
        end
        MultivariateVandermondeMatrix(data, exponents)
    end
end

# This function creates a function v.
# v(x) is the array with all the monomials in the entries of x of degree d
function veronese(exponents::Vector, ::Type{T}, n::Int)  where {T<:Number}
    v = map(exponents) do e
        FP.Polynomial(reshape(e, n, 1), [one(T)])
    end
    function (x::Vector)
        FP.evaluate(v,x)
    end
end

"""
FindEquations(M::MultivariateVandermondeMatrix,
              alg::Symbol)

Finds equations on input a multivariate Vandermonde matrix M.
* alg is one of :with_svd, :with_qr, :with_rref.

You can specifiy the tolerance.

FindEquations(M::MultivariateVandermondeMatrix,
             alg::Symbol,
             tol::Float64)
"""
# function that gets the equations from a MultivariateVandermondeMatrix
function FindEquations(M::MultivariateVandermondeMatrix, alg::Symbol)

    if alg == :with_svd
        m, N = size(M.Vandermonde)
        s = LinearAlgebra.svd(M.Vandermonde, full = true)
        tol = max(m,N) * maximum(s.S) * eps(Float64)
        rk = sum(s.S .> tol)
        return Polynomials_from_coefficients(s.Vt[rk + 1:end,:], M.exponents)
    elseif alg == :with_qr
        return Polynomials_from_coefficients(with_qr(M), M.exponents)
    elseif alg == :with_rref
        return Polynomials_from_coefficients(with_rref(M), M.exponents)
    else
        println("Method $(alg) not known.")
    end

end

function FindEquations(M::MultivariateVandermondeMatrix, alg::Symbol, tol::Float64)

    if alg == :with_svd
        m, N = size(M.Vandermonde)
        s = LinearAlgebra.svd(M.Vandermonde, full = true)
        rk = sum(s.S .> tol)
        return Polynomials_from_coefficients(s.Vt[rk + 1:end,:], M.exponents)
    elseif alg == :with_qr
        return Polynomials_from_coefficients(with_qr(M, tol), M.exponents)
    elseif alg == :with_rref
        return Polynomials_from_coefficients(with_rref(M, tol), M.exponents)
    else
        println("Method $(alg) not known.")
    end

end

# Creates a polynomial from a coefficient vector
function Polynomials_from_coefficients(kernel::Matrix{T}, exponents::Vector) where {T<:Number}
    tol = 1e-10
    l = size(kernel,1)
    nvars = length(exponents[1])
    DynamicPolynomials.@polyvar x_[1:nvars]

    if l == 0
        return [0]
    else
        map(1:l) do i
            non_zero_coeffs = findall(x -> abs(x) > tol, kernel[i,:])
            if length(non_zero_coeffs) > 0
                monomial = map(c -> prod(map(i -> x_[i]^exponents[c][i], 1:nvars)), non_zero_coeffs)
                return LinearAlgebra.dot(kernel[i,non_zero_coeffs],  monomial)
            else
                return 0.0*x[1]
            end
        end
    end
end

function with_qr(M::MultivariateVandermondeMatrix)
    R = LinearAlgebra.qr(M.Vandermonde).R
    m, N = size(R)
    s = LinearAlgebra.svdvals(R)
    tol = max(m,N) * maximum(s) * eps(Float64)
    return kernel_qr(R, tol)
end

function with_qr(M::MultivariateVandermondeMatrix, tol::Float64)
    R = LinearAlgebra.qr(M.Vandermonde).R
    return kernel_qr(R, tol)
end

function kernel_qr(R::Array{T,2}, tol::Float64) where {T <: Number}
    n,m = size(R)

    if n>=m
        index = findall(x -> abs(x) < tol, [R[i,i] for i in 1:m])
        index2 = setdiff([i for i in 1:m], index)
        R_small = R[:,index2]

        kernel = zeros(eltype(R), length(index), m)

        for i = 1:length(index)
            kernel[i,index[i]] = one(eltype(R))
            kernel[i,index2] =  (-1) .* R_small\R[:,index[i]]
        end
        return kernel
    else
        index = findall(x -> abs(x) < tol, [R[i,i] for i in 1:n])
        index2 = setdiff([i for i in 1:m], index)
        R_small = R[:,index2]

        kernel = zeros(eltype(R), length(index), m)

        for i = 1:length(index)
            kernel[i,index[i]] = one(eltype(R))
            kernel[i,index2] =  (-1) .* R_small\R[:,index[i]]
        end
        R_new = [R; LinearAlgebra.pinv(kernel) * kernel]
        s = LinearAlgebra.svd(R_new, full = true)
        rk = sum(s.S .> tol)
        return  [kernel; s.Vt[rk + 1:end,:]]
    end
end

function with_rref(M::MultivariateVandermondeMatrix)
    R = RowEchelon.rref(M.Vandermonde)
    m, N = size(R)
    s = LinearAlgebra.svdvals(R)
    tol = max(m,N) * maximum(s) * eps(Float64)
    return kernel_rref(R, tol)
end
function with_rref(M::MultivariateVandermondeMatrix, tol::Float64)
    R = RowEchelon.rref(M.Vandermonde)
    return kernel_rref(R, tol)
end
function kernel_rref(R::Array{T,2}, tol::Float64) where {T <: Number}
    n,m = size(R)
    R = R[findall([LinearAlgebra.norm(R[i,:]) for i in 1:n] .> sqrt(m) * tol),:]
    rk = size(R,1)
    index = zeros(Int64, rk, 2)
    for i = 1:rk
        where_are_the_ones = findall(abs.(R[i,:]).> tol)
        index[i,:] = [i where_are_the_ones[1]]
    end
    pivots = setdiff([i for i in 1:m], index[:,2])

    kernel = zeros(eltype(R), m-rk, m)
    for i = 1:(m-rk)
        t = findall(index[:,2] .< pivots[i])
        kernel[i,pivots[i]] = 1
        for j in t
            kernel[i, index[j,2]] = - R[index[j,1],pivots[i]]
        end
    end
    return kernel
end

ext: jl


text:
# Skoro mamy podobne czasy dla stopow = {32, 64, 128} sprawdzmy jakie sa bledy przy kazdym ze stopow
# rozmiary = indeksy, iteracjie = 20, rysuje wykres err * 10^18, wyniki zapisane w BladInvStrassStop2.png , z tego  BladStrass_zaleznosc_od_zmiany.png było dobrze też ale źle zapisałem wyniki :/
# te rzeczy wychodzą zupełnie z czapy
indeksy = [128, 256, 512, 1024, 2048]
a = [3.28352e-24 8.94986e-23 3.08026e-21 6.2254e-19 1.33283e-19; 1.1329e-24 1.15167e-22 3.28867e-21 3.57518e-19 5.32419e-20; 1.30968e-24 3.06569e-22 2.94942e-21 1.52337e-19 1.14285e-19]
ai = [9.8256e-25 8.17142e-23 1.28578e-21 1.47343e-19 1.43099e-19; 6.47901e-25 6.32194e-23 1.46342e-22 1.88038e-19 1.39254e-19; 6.50748e-25 3.91214e-23 8.40869e-23 3.92073e-19 1.47434e-19]

max_stars_repo_path: II rok/I semestr/ANM/Pracownia 1/prog/wynikibladInvstrassstop.jl

max_stars_repo_name: antonitomaszewski/Studia

max_issues_repo_name: antonitomaszewski/Studia

max_forks_repo_name: antonitomaszewski/Studia


text:
using Intervals
using Dates
using LibPQ

# Determine if postgres is installed
postgres_installed = try
    run(`createdb test-yesql`)
    true
catch
    @info "Not running LibPQ tests because PostgreSQL is not installed/configured"
    false
end

if postgres_installed
    @testset "Test LibPQ query execution" begin
        queries = load_path(joinpath(@__DIR__, "sql"))

        conn = LibPQ.Connection("""dbname=test-yesql user=$(ENV["USER"]) host=localhost""")

        # create a table
        queries[:create_table](conn)

        # insert a musician with dict params
        adrian = Dict(:first_name => "Adrian",
                      :last_name => "Belew",
                      :instrument =>"guitar",
                      :active => Date(1977, 1,1)..Date(1978,12,31))
        (id, ) = queries[:insert_musician_id](conn, adrian)
        @test id == 1

        id = queries[:get_musician_id](conn, adrian)
        @test id == 1

        # test returning value
        queries[:delete_musician](conn, adrian)
        id = queries[:insert_musician_id_value](conn, adrian)
        @test id == 2

        # select data about a musician
        (first_name, last_name) = queries[:get_musician_name](conn, (musician_id=id,))
        @test first_name == adrian[:first_name]
        @test last_name == adrian[:last_name]

        # insert multiple musicians
        guitarists = (first_name=["Steve", "Mike"], last_name=["Vai", "Keneally"],
                      instrument=["guitar" for _ in 1:2],
                      active=[Date(1981, 1, 1)..Date(1987,12,31), Date(1988, 1, 1)..Date(1992,12,31)])
        queries[:insert_musicians](conn, guitarists)

        queries[:insert_musician](conn, (first_name="Pierre", last_name="Boulez", instrument="Conductor",
                                         active=Date(1984,1,1)..Date(1985,1,1)))

        # count, delete and count again
        @test queries[:count_musicians](conn) == 4
        queries[:delete_musician](conn, adrian)
        @test queries[:count_musicians](conn) == 3

        # test that adrian is deleted and nothing is returned
        id = queries[:get_musician_id](conn, adrian)
        @test id === missing
        row = queries[:get_musician_name](conn, (musician_id=1,))
        @test row === missing

        # test arrays of ranges
        queries[:create_test_types](conn)

        d1 = Date(2001,1,1)
        d2 = Date(2002,1,1)
        params = (int_val=4,
                  real_val=4.0,
                  double_val=4.0,
                  ts_val=now(),
                  intrange_val=4..10,
                  numrange_val=3.0..5.0,
                  daterange_val=d1..d2,
                  tsrange_val=(now()-Day(1))..now(),
                  int_array=[1,2],
                  real_array=Float32[1,2,3],
                  intrange_array=[4..10, Interval{Closed,Open}(4,10)],
                  daterange_array=[d1..d2, Interval{Closed,Open}(d1,d2)])

        queries[:insert_test_types](conn, params)

        close(conn)
        run(`dropdb test-yesql`)
    end
end


text:
export hs0, hs

"""
`hs0(x)` returns the Rufus Isaacs derived half sine
that is valid on [-π,π].
"""
hs0(x) = (x>=0) ? -x : sin(-x)

"""
`hs(x)` is a continuous half-sine function created by extending
`hs0`.
"""
function hs(x::T) where T<:Real
    x = _reset(x)
    return hs0(x)
end

"""
`_reset(x)` is `asin(sin(x))` to map `x`
into the domain of `h0`.
"""
function _reset(x::T) where T<:Real
    x = asin(sin(x))
end

function make_picture()
    plot(hs,-4pi,4pi,aspect_ratio=1,legend=false,color=:black)
    draw_xaxis()
    draw_yaxis()

    finish()
end

function make_picture_1()
    x = -pi/2 : 0.05 : pi/2
    y = hs.(x)

    plot(x,y,aspect_ratio=1,legend=false, color=:black, linewidth=1.5)
    draw_vector(-2,0,color=:black,linewidth=0.5)
    draw_vector(2,0,color=:black,linewidth=0.5)
    draw_vector(0,1.2,color=:black,linewidth=0.5)
    draw_vector(0,-1.7,color=:black,linewidth=0.5)

    # draw_segment(-pi/2,-0.1,-pi/2,0.1,color=:black,linewidth=0.5)
    # draw_segment(pi/2,-0.1,pi/2,0.1,color=:black,linewidth=0.5)

    draw_xaxis()
    draw_yaxis()

    finish()
end


text:
module QXContexts

using Reexport

include("logger.jl")
include("parameters.jl")
include("compute_graph/compute_graph.jl")
include("contexts/contexts.jl")
include("sampling.jl")
include("execution.jl")
include("sysimage/sysimage.jl")

@reexport using QXContexts.Logger
@reexport using QXContexts.Param
@reexport using QXContexts.ComputeGraphs
@reexport using QXContexts.Contexts
@reexport using QXContexts.Sampling

end


text:
using LightGraphs

function sparse_beta_loglik(g::Graph, mu, beta)
    tot = 0
    for i in 1:n
        tot -= degree(g)[i] * beta[i]
        for j in i:n
            tot += log(1 + exp(mu + beta[i] + beta[j]))
        end
    end
    return -(g.ne * mu) + tot
end

function sparse_beta_conditional_loglik(g::Graph, mu, beta, i,)
    beta_i = beta[i]
    tot = degree(g)[i] * beta_i
    for b in beta[1:end .!= i]
        tot -= log1pexp(mu + beta_i + b)
    end
    return tot
end

function sparse_beta_conditional_loglik(g::Graph, mu, beta_i, beta_mi, deg)
    tot = deg * beta_i
    for b in beta[1:end .!= i]
        tot -= log1pexp(mu + beta_i + b)
    end
    return tot
end

max_stars_repo_name: jacobusmmsmit/sparse-beta

max_issues_repo_name: jacobusmmsmit/sparse-beta

max_forks_repo_name: jacobusmmsmit/sparse-beta


text:
using StatsBase
using RDatasets
a = [123,4234,23423,1231231,1432432423,1341413413]
var(a)
var(a, 2)
std(a)
mean_and_var(a)
mean_and_std(a)
skewness(a)
kurtosis(a)
moment(a,3)
span(a)
variation(a)

max_stars_repo_path: Chapter 3/5.jl

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_path: Chapter 3/5.jl

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_path: Chapter 3/5.jl

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
module PPPack

include("../deps/deps.jl")

function interv(xt::Vector{Float32},
                lxt::Integer, x::Float32,
                left::Integer,mflag::Integer)
    left = Int32[left]
    mflag = Int32[mflag]
    ccall((:interv_, _jl_libpppack),
          Void,
          (Ref{Float32}, Ref{Int32}, Ref{Float32}, Ref{Int32}, Ref{Int32}),
          xt, Int32(lxt), x, left, mflag)
    left[1],mflag[1]
end

function bsplvb(t::Vector{Float32}, jhigh::Integer, index::Integer,
                  x::Float32, left::Integer, biatx::AbstractVector{Float32})
    ccall((:bsplvb_, _jl_libpppack),
          Void,
          (Ref{Float32}, Ref{Int32}, Ref{Int32}, Ref{Float32}, Ref{Int32}, Ref{Float32}),
          t, Int32(jhigh), Int32(index), x, Int32(left), biatx)
end

function bsplvd(t::Vector{Float32}, k::Integer, x::Float32,
                left::Integer, a::Matrix{Float32},
                dbiatx::Matrix{Float32},
                nderiv::Integer)
    ccall((:bsplvd_, _jl_libpppack),
          Void,
          (Ref{Float32}, Ref{Int32}, Ref{Float32},
           Ref{Int32}, Ref{Float32},
           Ref{Float32},
           Ref{Int32}),
          t, Int32(k), x, Int32(left), a, dbiatx, nderiv)
end

pppbindir = abspath(joinpath(dirname(_jl_libpppack), "../../builds/pppack/bin"))
example(exe) = `$(joinpath(pppbindir, exe))`

export interv, bsplvb, bsplvd

end # module


text:
# Abstract interpreter semantics

"Abstract (Cartesian) PDDL interpreter."
@kwdef struct AbstractInterpreter <: Interpreter
    abstractions::Dict{Symbol,Any} = DEFAULT_ABSTRACTIONS
    autowiden::Bool = true
end

include("utils.jl")
include("domain.jl")
include("interface.jl")
include("satisfy.jl")
include("initstate.jl")
include("update.jl")

"""
    abstracted(domain; options...)

Construct an abstract domain from a concrete domain.
"""
abstracted(domain::GenericDomain; options...) =
    AbstractedDomain(domain; options...)

"""
    abstracted(domain, state; options...)
    abstracted(domain, problem; options...)

Construct an abstract domain and state from a concrete `domain` and `state`.
A `problem` can be provided instead of a `state`.
"""
function abstracted(domain::GenericDomain, state::GenericState; options...)
    absdom = abstracted(domain; options...)
    return (absdom, abstractstate(absdom, state))
end

abstracted(domain::GenericDomain, problem::GenericProblem) =
    abstracted(domain, initstate(domain, problem))

"""
    abstractstate(domain, state)

Construct a state in an abstract `domain` from a concrete `state.
"""
function abstractstate(domain::AbstractedDomain, state::GenericState)
    # Copy over facts
    abs_state = GenericState(copy(state.types), copy(state.facts))
    # Abstract non-Boolean values if necessary
    funcsigs = get_functions(domain)
    if isempty(funcsigs) return abs_state end
    abstypes = domain.interpreter.abstractions
    for (term, val) in get_fluents(state)
        if is_pred(term, domain) continue end
        type = funcsigs[term.name].type
        val = get(abstypes, type, identity)(val)
        abs_state[term] = val
    end
    return abs_state
end

abstractstate(domain::AbstractedDomain, problem::GenericProblem) =
    initstate(domain, problem)


text:
# Example for basic 2D modeling:
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

using JUDI, SegyIO, LinearAlgebra, PyPlot

# Set up model structure
n = (100, 100, 80)   # (x,y,z) or (x,z)
d = (10., 10., 10.)
o = (0., 0., 0.)

# Velocity [km/s]
v = ones(Float32,n) .+ 0.4f0
v[:, :, Int(round(end/2)):end] .= 3f0
v0 = ones(Float32,n) .+ 0.4f0

# Slowness squared [s^2/km^2]
m0 = (1f0 ./ v0).^2
m = (1f0 ./ v).^2
dm = vec(m - m0)

# Setup info and model structure
nsrc = 1	# number of sources
model0 = Model(n, d, o, m0)
model = Model(n, d, o, m)

# Receiver geometry
nxrec = 120
nyrec = 100
xrec = range(50f0, stop=950f0, length=nxrec)
yrec = range(100f0, stop=900f0, length=nyrec)
zrec = 50f0

# Construct 3D grid from basis vectors
(xrec, yrec, zrec) = setup_3D_grid(xrec, yrec, zrec)

# receiver sampling and recording time
time = 50f0   # receiver recording time [ms]
dt = 4f0    # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dt, t=time, nsrc=nsrc)

# Source wavelet
f0 = 0.01f0     # kHz
wavelet = ricker_wavelet(time, dt, f0)

# Set up info structure for linear operators
ntComp = get_computational_nt(recGeometry, model)
info = Info(prod(n), nsrc, ntComp)

###################################################################################################

# Write shots as segy files to disk
opt = Options(return_array=false)

# Setup operators
Pr = judiProjection(info, recGeometry)
F = judiModeling(info, model; options=opt)

# Random weights (size of the model)
weights = randn(Float32, model.n)
w = judiWeights(weights)

# Create operator for injecting the weights, multiplied by the provided wavelet(s)
Pw = judiLRWF(info, wavelet)

# Model observed data w/ extended source
F = Pr*F*adjoint(Pw)

# Simultaneous observed data
d_sim = F*w
dw = adjoint(F)*d_sim

# Jacobian
J = judiJacobian(F, w)
d_lin = J*dm
g = adjoint(J)*d_lin


text:
# ---
# title: 437. Path Sum III
# id: problem437
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree
# link: <https://leetcode.com/problems/path-sum-iii/description/>
# hidden: true
# ---
#
# You are given a binary tree in which each node contains an integer value.
#
# Find the number of paths that sum to a given value.
#
# The path does not need to start or end at the root or a leaf, but it must go
# downwards (traveling only from parent nodes to child nodes).
#
# The tree has no more than 1,000 nodes and the values are in the range
# -1,000,000 to 1,000,000.
#
# **Example:**
#
#
#
#     root = [10,5,-3,3,2,null,11,3,-2,null,1], sum = 8
#
#           10
#          /  \
#         **5**   **-3**
#        **/** **\**    **\**
#       **3**   **2**   **11**
#      / \   **\**
#     3  -2   **1**
#
#     Return 3. The paths that sum to 8 are:
#
#     1.  5 -> 3
#     2.  5 -> 2 -> 1
#     3. -3 -> 11
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "OpenJpeg"
version = v"2.4.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/uclouvain/openjpeg/archive/v$(version)/openjpeg-$(version).tar.gz",
                  "8702ba68b442657f11aaeb2b338443ca8d5fb95b0d845757968a7be31ef7f16d"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/openjpeg-*/
for f in ${WORKSPACE}/srcdir/patches/*.patch; do
    atomic_patch -p1 ${f}
done
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DBUILD_STATIC_LIBS=OFF
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libopenjp2", :libopenjp2),
    ExecutableProduct("opj_decompress", :opj_decompress),
    ExecutableProduct("opj_dump", :opj_dump),
    ExecutableProduct("opj_compress", :opj_compress)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="LittleCMS_jll", uuid="d3a379c0-f9a3-5b72-a4c0-6bf4d2e8af0f"))
    Dependency(PackageSpec(name="libpng_jll", uuid="b53b4c65-9356-5827-b1ea-8c7a1a84506f"))
    Dependency("Libtiff_jll"; compat="4.3.0")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using ComplexVisual
@ComplexVisual.import_huge

const fontface = cv_fontface("sans-serif")

function axis1(layout, can_math_l, attach)
    t1, t2 = collect(-2:2), collect(-1.5:1.0:1.5)
    t3 = setdiff(setdiff(collect(-2:0.1:2), t1), t2)

    ticks, ticksh, tickss = "%.0f" ⇒ t1, "%.1f" ⇒ t2, "" ⇒ t3

    app = CV_TickLabelAppearance(; tick_length=10,
        label_style=cv_color(0, 0, 1) → fontface → cv_fontsize(20))
    apph = CV_TickLabelAppearance(; tick_length=6,
        label_style=cv_color(.4, .4, 0) → fontface → cv_fontsize(10))
    apps = CV_TickLabelAppearance(; tick_length=3)
    return cv_ticks_labels(layout, can_math_l, attach,
        (app ↦ ticks, apph ↦ ticksh, apps ↦ tickss))
end

function main()
    layout = CV_2DLayout()

    can_math = CV_Math2DCanvas(-2.0 + 2.0im, 2.0 - 2.0im, 100)
    can_math_l = cv_add_canvas!(layout, can_math,
        cv_anchor(can_math, :center), (0,0))
    cc_can_math = cv_create_context(can_math)

    a1_l = axis1(layout, can_math_l, cv_south)
    a2_l = axis1(layout, can_math_l, cv_north)
    a3_l = axis1(layout, can_math_l, cv_east)
    a4_l = axis1(layout, can_math_l, cv_west)

    cv_add_padding!(layout, 30)

    fill_red = cv_color(0.7,0.1,0.1) ↦ CV_FillPainter()

    can_layout = cv_canvas_for_layout(layout)
    cc_can_layout = cv_create_context(can_layout)

    cv_paint(cc_can_math, fill_red)
    can_math_l(cc_can_layout)
    a1_l(cc_can_layout)
    a2_l(cc_can_layout)
    a3_l(cc_can_layout)
    a4_l(cc_can_layout)

    cv_save_image(can_layout, "./Axis01.png")
end

main()

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4: