function testsquare1()::Matrix{Int}
    [1 0 4 0;
     0 2 0 3;
     0 0 3 0;
     0 0 0 4]
end
function testsquare2()::Matrix{Int}
    [4 0 0 2;
     0 3 0 0;
     0 0 2 0;
     3 0 0 1]
end

function neighbors_4x4(m::Matrix{Int}, i::Int, j::Int)::Vector{Int}
    q, r = div(i-1, 2), div(j-1, 2)
    s = union(Set(m[i,:]), Set(m[:,j]), Set([m[2q+1, 2r+1], m[2q+1, 2r+2], m[2q+2, 2r+1], m[2q+2, 2r+2]]))
    pop!(s, 0)
    collect(s)
end

function candidates_4x4(m::Matrix{Int}, i::Int, j::Int)::Vector{Int}
    n = neighbors_4x4(m, i, j)
    c = setdiff(1:4, n)
    collect(c)
end

function solve_4x4(m::Matrix{Int})::Union{Some{Matrix{Int}}, Nothing}
    @assert all(m .<= 4) && all(m .>= 0)

    m = copy(m)
    n = 1
    d = Matrix{Union{Some{Vector{Int}}, Nothing}}(nothing, 4, 4)
    # 1) Find candidates for all positions.
    while n > 0
        n = 0
        for i = 1:4
            for j = 1:4
                if m[i,j] == 0
                    d[i,j] = Some(candidates_4x4(m, i, j))
                end
            end
        end
        # 2) Fix those cells where only one candidate matches
        for i = 1:4
            for j = 1:4
                if m[i,j] == 0 && length(something(d[i,j])) == 1
                    n += 1
                    m[i,j] = something(d[i,j])[1]
                    d[i,j] = nothing
                elseif m[i,j] == 0 && length(something(d[i,j])) == 0
                    display(m)
                    error("Unsolveable magic square! at $((i,j))")
                end
            end
        end
    end
   
    if all(m .> 0)
        return Some(m)
    end
    # 3) Backtracking step: Select candidate, try to solve.
    for i = 1:4
        for j = 1:4
            if !isnothing(d[i,j])
                for cand in something(d[i,j])
                    m[i,j] = cand
                    m_ = solve_4x4(m)
                    if !isnothing(m_)
                        return m_
                    end
                end
                return nothing
            end
        end
    end
end