abstract type DAEBDFMutableCache <: OrdinaryDiffEqMutableCache end
function get_fsalfirstlast(cache::DAEBDFMutableCache, u)
    return (cache.fsalfirst, du_alias_or_new(cache.nlsolver, cache.fsalfirst))
end

@cache mutable struct DImplicitEulerCache{uType, rateType, uNoUnitsType, N} <:
    DAEBDFMutableCache
    u::uType
    uprev::uType
    uprev2::uType
    atmp::uNoUnitsType
    k₁::rateType
    k₂::rateType
    nlsolver::N
end

# Not FSAL
get_fsalfirstlast(cache::DImplicitEulerCache, u) = (nothing, nothing)

mutable struct DImplicitEulerConstantCache{N} <: OrdinaryDiffEqConstantCache
    nlsolver::N
end

function alg_cache(
        alg::DImplicitEuler, du, u, res_prototype, rate_prototype,
        ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits},
        ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
        ::Val{false}, verbose
    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
    γ, c = 1, 1
    α = 1
    nlsolver = build_nlsolver(
        alg, u, uprev, p, t, dt, f, res_prototype, uEltypeNoUnits,
        uBottomEltypeNoUnits, tTypeNoUnits, γ, c, α, Val(false), verbose
    )

    return DImplicitEulerConstantCache(nlsolver)
end

function alg_cache(
        alg::DImplicitEuler, du, u, res_prototype, rate_prototype,
        ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits},
        ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
        ::Val{true}, verbose
    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
    γ, c = 1, 1
    α = 1
    k₁ = zero(rate_prototype)
    k₂ = zero(rate_prototype)
    nlsolver = build_nlsolver(
        alg, u, uprev, p, t, dt, f, res_prototype, uEltypeNoUnits,
        uBottomEltypeNoUnits, tTypeNoUnits, γ, c, α, Val(true), verbose
    )

    atmp = similar(u, uEltypeNoUnits)
    recursivefill!(atmp, false)

    return DImplicitEulerCache(u, uprev, uprev2, atmp, k₁, k₂, nlsolver)
end

@cache mutable struct DABDF2ConstantCache{N, dtType, rate_prototype} <:
    OrdinaryDiffEqConstantCache
    nlsolver::N
    eulercache::DImplicitEulerConstantCache
    dtₙ₋₁::dtType
    fsalfirstprev::rate_prototype
end

function alg_cache(
        alg::DABDF2, du, u, res_prototype, rate_prototype,
        ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits},
        ::Type{tTypeNoUnits},
        uprev, uprev2, f, t, dt, reltol, p, calck,
        ::Val{false}, verbose
    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
    γ, c = Int64(1) // 1, 1
    α = Int64(1) // 1
    nlsolver = build_nlsolver(
        alg, u, uprev, p, t, dt, f, res_prototype, uEltypeNoUnits,
        uBottomEltypeNoUnits, tTypeNoUnits, γ, c, α, Val(false), verbose
    )
    eulercache = DImplicitEulerConstantCache(nlsolver)

    dtₙ₋₁ = one(dt)
    fsalfirstprev = rate_prototype

    return DABDF2ConstantCache(nlsolver, eulercache, dtₙ₋₁, fsalfirstprev)
end

@cache mutable struct DABDF2Cache{uType, rateType, uNoUnitsType, N, dtType} <:
    DAEBDFMutableCache
    uₙ::uType
    uₙ₋₁::uType
    uₙ₋₂::uType
    fsalfirst::rateType
    fsalfirstprev::rateType
    atmp::uNoUnitsType
    nlsolver::N
    eulercache::DImplicitEulerCache
    dtₙ₋₁::dtType
end

function alg_cache(
        alg::DABDF2, du, u, res_prototype, rate_prototype,
        ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits},
        ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
        ::Val{true}, verbose
    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
    γ, c = Int64(1) // 1, 1
    α = Int64(1) // 1
    nlsolver = build_nlsolver(
        alg, u, uprev, p, t, dt, f, res_prototype, uEltypeNoUnits,
        uBottomEltypeNoUnits, tTypeNoUnits, γ, c, α, Val(true), verbose
    )
    fsalfirst = zero(rate_prototype)

    fsalfirstprev = zero(rate_prototype)
    atmp = similar(u, uEltypeNoUnits)
    recursivefill!(atmp, false)

    k₁ = zero(rate_prototype)
    k₂ = zero(rate_prototype)

    eulercache = DImplicitEulerCache(u, uprev, uprev2, atmp, k₁, k₂, nlsolver)

    dtₙ₋₁ = one(dt)

    return DABDF2Cache(
        u, uprev, uprev2, fsalfirst, fsalfirstprev, atmp,
        nlsolver, eulercache, dtₙ₋₁
    )
end

@cache mutable struct DFBDFConstantCache{
        MO, N, tsType, tType, uType, uuType, coeffType,
        EEstType, rType, wType, fdWeightsType,
    } <:
    OrdinaryDiffEqConstantCache
    nlsolver::N
    ts::tsType
    ts_tmp::tsType
    t_old::tType
    u_history::uuType
    order::Int
    prev_order::Int
    u₀::uType
    u_corrector::uuType
    bdf_coeffs::coeffType
    max_order::Val{MO}
    nconsteps::Int
    consfailcnt::Int
    qwait::Int # countdown to next order change consideration (CVODE-style)
    terkm2::EEstType
    terkm1::EEstType
    terk::EEstType
    terkp1::EEstType
    r::rType
    weights::wType
    iters_from_event::Int
    fd_weights::fdWeightsType
end

function alg_cache(
        alg::DFBDF{MO}, du, u, res_prototype, rate_prototype, uEltypeNoUnits,
        uBottomEltypeNoUnits, tTypeNoUnits,
        uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{false}, verbose
    ) where {MO}
    γ, c = 1.0, 1.0
    max_order = MO
    nlsolver = build_nlsolver(
        alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits,
        uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(false), verbose
    )
    bdf_coeffs = _make_bdf_coeffs_dfbdf()
    ts = zero(Vector{typeof(t)}(undef, max_order + 2)) #ts is the successful past points, it will be updated after successful step
    ts_tmp = similar(ts)

    if u isa Number
        u_history = zeros(eltype(u), max_order + 2)
        u_corrector = zeros(eltype(u), max_order + 2)
    else
        u_history = [zero(u) for _ in 1:(max_order + 2)]
        u_corrector = [zero(u) for _ in 1:(max_order + 2)]
    end
    order = 1
    prev_order = 1
    terkm2 = tTypeNoUnits(1)
    terkm1 = tTypeNoUnits(1)
    terk = tTypeNoUnits(1)
    terkp1 = tTypeNoUnits(1)
    r = zero(Vector{typeof(t)}(undef, max_order + 2))
    weights = zero(Vector{typeof(t)}(undef, max_order + 2))
    weights[1] = 1
    nconsteps = 0
    consfailcnt = 0
    qwait = 3 # order + 2, matching nconsteps >= order + 2 for failure-free runs
    t_old = zero(t)
    iters_from_event = 0
    u₀ = zero(u)

    fd_weights = zeros(typeof(t), max_order + 1, max_order + 1)

    return DFBDFConstantCache(
        nlsolver, ts, ts_tmp, t_old, u_history, order, prev_order, u₀,
        u_corrector, bdf_coeffs, Val(MO), nconsteps, consfailcnt, qwait, terkm2,
        terkm1, terk, terkp1, r, weights, iters_from_event, fd_weights
    )
end

@cache mutable struct DFBDFCache{
        MO, N, rateType, uNoUnitsType, tsType, tType, uType,
        uuType, coeffType, EEstType, rType, wType, fdWeightsType,
    } <:
    DAEBDFMutableCache
    fsalfirst::rateType
    nlsolver::N
    ts::tsType
    ts_tmp::tsType
    t_old::tType
    u_history::uuType
    order::Int
    prev_order::Int
    u_corrector::uuType
    u₀::uType
    bdf_coeffs::coeffType
    max_order::Val{MO}
    nconsteps::Int
    consfailcnt::Int
    qwait::Int # countdown to next order change consideration (CVODE-style)
    tmp::uType
    atmp::uNoUnitsType
    terkm2::EEstType
    terkm1::EEstType
    terk::EEstType
    terkp1::EEstType
    terk_tmp::uType
    terkp1_tmp::uType
    r::rType
    weights::wType
    equi_ts::tsType
    iters_from_event::Int
    dense::Vector{uType}
    fd_weights::fdWeightsType
end

function alg_cache(
        alg::DFBDF{MO}, du, u, res_prototype, rate_prototype, uEltypeNoUnits,
        uBottomEltypeNoUnits,
        tTypeNoUnits, uprev, uprev2, f, t, dt, reltol, p, calck,
        ::Val{true}, verbose
    ) where {MO}
    γ, c = 1.0, 1.0
    fsalfirst = zero(rate_prototype)
    max_order = MO
    nlsolver = build_nlsolver(
        alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits,
        uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(true), verbose
    )
    bdf_coeffs = _make_bdf_coeffs_dfbdf()
    ts = Vector{typeof(t)}(undef, max_order + 2) #ts is the successful past points, it will be updated after successful step
    u_history = [zero(u) for _ in 1:(max_order + 2)]
    order = 1
    prev_order = 1
    u_corrector = [zero(u) for _ in 1:(max_order + 2)]
    recursivefill!(ts, zero(t))
    terkm2 = tTypeNoUnits(1)
    terkm1 = tTypeNoUnits(1)
    terk = tTypeNoUnits(1)
    terkp1 = tTypeNoUnits(1)
    terk_tmp = similar(u)
    terkp1_tmp = similar(u)
    r = Vector{typeof(t)}(undef, max_order + 2)
    weights = Vector{typeof(t)}(undef, max_order + 2)
    recursivefill!(r, zero(t))
    recursivefill!(weights, zero(t))
    weights[1] = 1
    nconsteps = 0
    consfailcnt = 0
    qwait = 3 # order + 2, matching nconsteps >= order + 2 for failure-free runs
    t_old = zero(t)
    atmp = similar(u, uEltypeNoUnits)
    recursivefill!(atmp, zero(uEltypeNoUnits))
    u₀ = similar(u)
    equi_ts = similar(ts)
    tmp = similar(u)
    ts_tmp = similar(ts)
    iters_from_event = 0

    dense = [zero(u) for _ in 1:(2 * (max_order + 1))]  # first half for integrator.k, second half as scratch

    fd_weights = zeros(typeof(t), max_order + 1, max_order + 1)

    return DFBDFCache(
        fsalfirst, nlsolver, ts, ts_tmp, t_old, u_history, order, prev_order,
        u_corrector, u₀, bdf_coeffs, Val(MO), nconsteps, consfailcnt, qwait, tmp, atmp,
        terkm2, terkm1, terk, terkp1, terk_tmp, terkp1_tmp, r, weights, equi_ts,
        iters_from_event, dense, fd_weights
    )
end

############################################ DNordsieckBDF
@cache mutable struct DNordsieckBDFCache{
        MO, N, rateType, uNoUnitsType, uType, tType, coeffType, staldType,
    } <: BDFMutableCache
    fsalfirst::rateType
    nlsolver::N
    zn::Vector{uType}
    u₀::uType
    ypred::uType
    acor::uType
    tempv::uType
    tmp::uType
    atmp::uNoUnitsType
    l::coeffType
    tau::coeffType
    tq::coeffType
    order::Int
    qprime::Int
    qwait::Int
    nst::Int
    nef::Int
    ncf::Int
    index_acor::Int
    max_order::Val{MO}
    max_order_int::Int
    hscale::tType
    eta::tType
    etamax::tType
    etaq::tType
    etaqm1::tType
    etaqp1::tType
    saved_tq5::tType
    predicted::Bool
    stald::staldType
end

@truncate_stacktrace DNordsieckBDFCache 1

mutable struct DNordsieckBDFConstantCache{MO, N, uType, tType, coeffType, staldType} <:
    OrdinaryDiffEqConstantCache
    nlsolver::N
    zn::Vector{uType}
    u₀::uType
    ypred::uType
    acor::uType
    l::coeffType
    tau::coeffType
    tq::coeffType
    order::Int
    qprime::Int
    qwait::Int
    nst::Int
    nef::Int
    ncf::Int
    index_acor::Int
    max_order::Val{MO}
    max_order_int::Int
    hscale::tType
    eta::tType
    etamax::tType
    etaq::tType
    etaqm1::tType
    etaqp1::tType
    saved_tq5::tType
    predicted::Bool
    stald::staldType
end

const NordsieckCaches = Union{
    NordsieckBDFCache, NordsieckBDFConstantCache,
    DNordsieckBDFCache, DNordsieckBDFConstantCache,
}

function alg_cache(
        alg::DNordsieckBDF{MO}, du, u, res_prototype, rate_prototype,
        ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
        uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}, verbose
    ) where {MO, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
    γ, c = one(tTypeNoUnits), one(tTypeNoUnits)
    nlsolver = build_nlsolver(
        alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits,
        uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(true), verbose
    )
    zn = [zero(u) for _ in 1:(MO + 1)]
    coeffs() = zeros(typeof(t), MO + 3)
    stald = StabilityLimitDetectionState(real(uBottomEltypeNoUnits); enabled = alg.stald)
    return DNordsieckBDFCache{
        MO, typeof(nlsolver), typeof(rate_prototype),
        typeof(similar(u, uEltypeNoUnits)), typeof(u), typeof(t), typeof(coeffs()),
        typeof(stald),
    }(
        zero(rate_prototype), nlsolver, zn, zero(u), zero(u), zero(u), zero(u),
        zero(u), similar(u, uEltypeNoUnits), coeffs(), coeffs(), zeros(typeof(t), 6),
        1, 1, 2, 0, 0, 0, MO, Val(MO), MO,
        zero(t), one(t), typeof(t)(NORD_ETA_MAX_FS), one(t), zero(t), zero(t),
        zero(t), false, stald
    )
end

function alg_cache(
        alg::DNordsieckBDF{MO}, du, u, res_prototype, rate_prototype,
        ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
        uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{false}, verbose
    ) where {MO, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
    γ, c = one(tTypeNoUnits), one(tTypeNoUnits)
    nlsolver = build_nlsolver(
        alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits,
        uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(false), verbose
    )
    zn = [zero(u) for _ in 1:(MO + 1)]
    coeffs() = zeros(typeof(t), MO + 3)
    stald = StabilityLimitDetectionState(real(uBottomEltypeNoUnits); enabled = alg.stald)
    return DNordsieckBDFConstantCache{
        MO, typeof(nlsolver), typeof(u), typeof(t), typeof(coeffs()), typeof(stald),
    }(
        nlsolver, zn, zero(u), zero(u), zero(u), coeffs(), coeffs(),
        zeros(typeof(t), 6), 1, 1, 2, 0, 0, 0, MO, Val(MO), MO,
        zero(t), one(t), typeof(t)(NORD_ETA_MAX_FS), one(t), zero(t), zero(t),
        zero(t), false, stald
    )
end

# With `adaptive = false` the integrator never calls `step_accept_controller!`, so
# the array has to be committed here instead. The order still ramps up to
# `max_order` through the usual `qwait` countdown, which keeps the saved
# correction column valid for the order increase.
