function largest_plus(a::AbstractMatrix{Bool})::Int
    is_still_plus(i, j, n) = (i+n) <= size(a, 1) && (i-n) >= 1 && (j+n) <= size(a, 2) && (j-n) >= 1 && a[i,j] && a[i+n,j] && a[i-n,j] && a[i,j+n] && a[i,j-n]
    is_plus_center(i, j) = is_still_plus(i, j, 1)

    initial_centers = [(i, j) for i = 2:size(a, 1)-1 for j = 2:size(a, 2)-1 if is_plus_center(i, j)]

    centers = initial_centers
    for n = 2:size(a, 1)
        new_centers = [(i, j) for (i, j) = centers if is_still_plus(i, j, n)]
        if isempty(new_centers)
            return n-1
        end
        centers = new_centers
    end
    -1
end

"""
ACDB is interleaving of AB and CD
 
ADEBCF is interleaving of ABC and DEF
 
ACBCD is interleaving of ABC and CD
 
ACDABC is interleaving of ABC and ACD
"""
function string_is_interleaving(a::AbstractString, b::AbstractString, c::AbstractString)::Bool
    if length(a) == 0 && length(b) == 0 && length(c) == 0
        true
    else
        if !isempty(a) && a[begin] == c[begin]
            if length(a) == 1 && length(b) == 0
                length(c) == 1
            else
                string_is_interleaving(a[2:end], b, c[2:end])
            end
        elseif !isempty(b) && b[begin] == c[begin]
            if length(b) == 1 && length(a) == 0
                length(c) == 1
            else
                string_is_interleaving(a, b[2:end], c[2:end])
            end
        else
            false
        end
    end
end

function partition_3(s::AbstractSet{T}, s1=nothing, s2=nothing, s3=nothing, n=length(s)) where {T<:Integer}
    if isnothing(s1) || isnothing(s2) || isnothing(s3)
        s1, s2, s3 = Set{T}(), Set{T}(), Set{T}()
    end
    if length(s1) > 0 && length(s2) > 0 && length(s3) > 0 && length(s1)+length(s2)+length(s3) == n && sum(s1) == sum(s2) && sum(s2) == sum(s3)
        return (true, (s1, s2, s3))
    elseif isempty(s)
        return (false, (s1,s2,s3))
    end

    e = pop!(s)

    push!(s1, e)
    ok, sets = partition_3(s, s1, s2, s3, n)
    if ok
        return ok, sets
    end

    pop!(s1, e)
    push!(s2, e)
    ok, sets = partition_3(s, s1, s2, s3, n)
    if ok
        return ok, sets
    end

    pop!(s2, e)
    push!(s3, e)
    ok, sets = partition_3(s, s1, s2, s3, n)
    if ok
        return ok, sets
    end

    pop!(s3, e)
    push!(s, e)

    return false, (s1, s2, s3)
end

function total_possible_solutions(coefficients::AbstractVector{<:Integer}, rhs::Integer)::Int
    if isempty(coefficients)
        return rhs == 0 ? 1 : 0
    end
    c = 0
    for v = 0:rhs
        c += total_possible_solutions((@view coefficients[2:end]), rhs-coefficients[begin]*v)
    end
    c
end

function recursive_liss(v::AbstractVector{<:Integer})
    N = length(v)
    prev, L = zeros(Int, N), zeros(Int, N)
    Lmax = recursive_liss(v, 1, 1, 0, prev, L)

    # Reconstruct:
    i = L[Lmax]
    l = Lmax
    liss = zeros(Int, Lmax)
    while i > 0
        liss[l] = v[i]
        i = prev[i]
        l -= 1
    end
    liss
end

function recursive_liss(v::AbstractVector{<:Integer}, last=1, i=1, l=0, prev=zeros(Int, length(v)), L=zeros(Int, length(v)), Lmax=[0])::Int
    if i > length(v)
        return l
    else
        if v[i] >= v[last]
            incl = recursive_liss(v, i, i+1, l+1, prev, L, Lmax)
            excl = recursive_liss(v, last, i+1, l, prev, L, Lmax)
            if incl >= excl
                Lmax[1] = max(Lmax[1], incl)
            end
            if incl == Lmax[1]
                L[l+1] = i
                prev[i] = last
            end
            max(incl, excl)
        else # v[i] < v[last]
            incl = recursive_liss(v, i, i+1, 1, prev, L, Lmax)
            excl = recursive_liss(v, last, i+1, l, prev, L, Lmax)
            max(incl, excl)
        end
    end
end

function number_of_pattern_times(hs::AbstractString, pat::AbstractString)::Int
    if isempty(pat) || isempty(hs)
        0
    else
        if length(pat) == 1
            c = hs[begin] == pat[begin] ? 1 : 0
            c + number_of_pattern_times((@view hs[2:end]), pat)
        else
            if hs[begin] == pat[begin]
                number_of_pattern_times((@view hs[2:end]), (@view pat[2:end])) + number_of_pattern_times((@view hs[2:end]), pat)
            else
                number_of_pattern_times((@view hs[2:end]), pat)
            end
        end
    end
end

function binpack_simple(items::AbstractVector{<:Pair}, caps::AbstractVector{<:Integer}=Int[], i=1, bincap=5, bincost=1)
    @show (i, caps)
    if i > length(items)
        return 0
    end
    size, profit = items[i]

    bin_candidates = zeros(Int, length(caps))
    for (b, c) = enumerate(caps)
        if c >= size
            caps[b] -= size
            bin_candidates[b] = profit + binpack_simple(items, caps, i+1, bincap, bincost)
            caps[b] += size
        else
            continue
        end
    end

    maxbincand = isempty(bin_candidates) ? 0 : maximum(bin_candidates)
    newbin = profit + binpack_simple(items, vcat(caps, [bincap-size]), i+1, bincap, bincost) - bincost

    if maxbincand > newbin
        maxbincand
    else
        newbin
    end

end

function knapsack_one_track(cap::Integer, items::AbstractVector{<:Pair}, i=1, map::Vector{Bool}=zeros(Bool, length(items)))
    if i > length(items) || cap <= 0
        return zeros(Bool, length(items)), 0
    end

    current_size, current_profit = items[i]
    incl_map, incl_item = if current_size > cap
        zeros(Bool, length(items)), 0
    else
        map2, nxp = knapsack_one_track(cap - current_size, items, i+1, map)
        map2, current_profit + nxp
    end
    excl_map, excl_item = knapsack_one_track(cap, items, i+1, map)

    if incl_item > excl_item
        incl_map[i] = true
        (incl_map, incl_item)
    else
        (excl_map, excl_item)
    end
end

function knapsack_one(cap::Integer, sizes::AbstractVector{<:Integer}, profits::AbstractVector{<:Integer})::Int
    if isempty(sizes) || cap <= 0
        return 0
    end
    # Include first item:
    with_current = if cap >= sizes[begin]
            profits[begin] + knapsack_one(cap - sizes[begin], (@view sizes[2:end]), (@view profits[2:end]))
        else
            0
        end
    # or don't
    without_current = knapsack_one(cap, (@view sizes[2:end]), (@view profits[2:end]))
    max(without_current, with_current)
end

"""Maximum Matrix Points
https://www.techiedelight.com/collect-maximum-points-matrix-satisfying-given-constraints/
"""
function max_matrix_points(mat::AbstractMatrix{Int}, start=(1,1), P=0, Pmax=[0], Ps=fill((0,0), 100))::Int
    row, col = start
    if any(start .> size(mat))
        return 0
    end
    opts = if iseven(row)
        [(0, -1), (1, 0)]
    else
        [(0, 1), (1, 0)]
    end
    @show (start, opts)

    max_pts = mat[start...]
    for opt = opts
        ix = start .+ opt
        if any(ix .< (1,1)) || any(ix .> size(mat)) || mat[ix...] < 0
            continue
        end
        cand = mat[start...] + max_matrix_points(mat, ix, P+mat[start...], Pmax, Ps)
        max_pts = max(max_pts, cand)
        Pmax[1] = max(Pmax[1], P+max_pts)
        if max_pts == cand && (P+max_pts) == Pmax[1]
            Ps[P+1] = start
        end
    end
    if max_pts == mat[start...]
        Pmax[1] = max(Pmax[1], P+max_pts)
        if (P+max_pts) == Pmax[1]
            Ps[P+1] = start
        end
    end

    max_pts
end

function matmul_operations(a, b, c)::Int
    2 * (a*c) * b
end

function matmul_opt(matrices::AbstractVector{<:Pair}, startix=1, endix=length(matrices))::Pair{Int,Tuple}
    if length(matrices) <= 2
        a, b = (matrices[begin])
        b, c = (matrices[end])
        matmul_operations(a, b, c) => tuple(matrices...)
    elseif endix - startix == 1
        a, b = (matrices[startix])
        b, c = (matrices[endix])
        matmul_operations(a, b, c) => tuple(matrices[startix:endix]...)
    elseif endix == startix
        0 => tuple(matrices[startix])
    else
        for i = startix+1:endix
            @assert matrices[i-1][2] == matrices[i][1]
        end
        first_axis(p::Pair) = p[1]
        first_axis(t::Tuple) = first_axis(t[1])
        last_axis(p::Pair) = p[2]
        last_axis(t::Tuple) = last_axis(t[end])

        # Choices: first * (rest), (first * second) * rest
        minops, minsplit = 1e21, nothing
        # split *after* i
        for i = startix:endix-1
            leftops, left = matmul_opt(matrices, startix, i)
            rightops, right = matmul_opt(matrices, i+1, endix)
            ops = leftops + rightops + matmul_operations(first_axis(left), last_axis(left), last_axis(right))
            if leftops+rightops < ops
                minsplit = (left, right)
                minops = ops
            end
        end
        minops => minsplit
    end
end