text:
# Floating point arithmetic functions

for op ∈ (
    :(+), #:(-),
    :abs, :floor, :ceil, :round,
    :inv,# :log10, :log2, :exp2, :sin, :cos,
    :sqrt, :trunc
    # :exp, :log
    )
    rename = VECTOR_SYMBOLS[op]
    @eval begin
        @inline $rename(s1::FloatingTypes) = @fastmath $op(s1)

        @vectordef $rename function Base.$op(v1) where {W,T<:FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1))
        end
    end
end

@inline vsign(s1::FloatingTypes) = sign(s1)
@vectordef vsign function Base.sign(v1) where {W,T<:FloatingTypes}
    vifelse(
        visequal(extract_data(v1), vbroadcast(Vec{W,T}, zero(T))),
        vbroadcast(Vec{W,T}, zero(T)),
        vcopysign(vbroadcast(Vec{W,T}, one(T)), extract_data(v1))
    )
end
# @inline function vsign(v1::SVec{W,T}) where {W,T<:FloatingTypes}
#     SVec(vsign(extract_data(v1)))
# end
# @inline function Base.sign(v1::SVec{W,T}) where {W,T<:FloatingTypes}
#     SVec(vsign(extract_data(v1)))
# end

for op ∈ (
        :(+), :(-), :(*), :(/), :(%),# :(^),
        :copysign#, :max, :min
    )
    rename = VECTOR_SYMBOLS[op]
    # exact / explicit version
    erename = Symbol(:e, rename)
    @eval begin
        @inline $rename(s1::FloatingTypes, s2::FloatingTypes) = @fastmath $op(s1, s2)

        @vectordef $rename function Base.$op(v1, v2) where {W,T <: FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1), extract_data(v2))
        end

        @inline $erename(s1::FloatingTypes, s2::FloatingTypes) = $op(s1, s2)

        @evectordef $erename function Base.$op(v1, v2) where {W,T <: FloatingTypes}
            # llvmwrap_notfast(Val{$(QuoteNode(op))}(), extract_data(v1), extract_data(v2))
            llvmwrap_notfast(Val{$(QuoteNode(op))}(), v1, v2)
        end

        @inline function Base.$op(v1::AbstractStructVec{W,T1}, v2::AbstractStructVec{W,T2}) where {W,T1,T2}
            T = promote_type(T1, T2)
            $op(vconvert(SVec{W,T}, v1), vconvert(SVec{W,T}, v2))
        end
        @inline function $rename(v1::V1, v2::V2) where {W,V1<:AbstractSIMDVector{W},V2<:AbstractSIMDVector{W}}
            V = promote_vtype(V1, V2)
            $rename(vconvert(V, v1), vconvert(V, v2))
        end
    end
end
# Julia 1.0 was bad at resovling dispatch (it liked to pick the wrong method)
# On Julia 1.0, this was less specific:
# ::SVec{W,Float64}, ::SVec{W,Float64}
# Than:
# ::Union{Vec{W,T},AbstractStructVec{W,T}}, ::Union{Vec{W,T},AbstractStructVec{W,T}} where {T <: Union{Float64,Float32,...}
#
# Note that SVec{W,Float64} <: AbstractStructVec{W,Float64}.
# We have to make the type fully concrete (e.g., SVec{8,Float64}) before
# to get Julia 1.0 and 1.1 to dispatch on it.
# This is definitely not needed on Julia 1.4, so I place the bound there. I have not actually tested Julia 1.2 and Julia 1.3.
@static if VERSION ≤ v"1.4"
    for T ∈ [Float32, Float64, Int16, Int32, Int64, UInt16, UInt32, UInt64]
        for op ∈ (
            :(+), :(-), :(*), :(/), :(%),# :(^),
            :copysign#, :max, :min
        )
            rename = VECTOR_SYMBOLS[op]
            for W ∈ 1:16
                @eval @inline $rename(v1::SVec{$W,$T},v2::SVec{$W,$T}) = SVec($rename(extract_data(v1), extract_data(v2)))
            end
        end
        for W ∈ 1:16
            U = VectorizationBase.mask_type(W)
            @eval @inline vadd(m::Mask{$W,$U}, v::SVec{$W,$T}) = vifelse(m, vadd(v, one(v)), v)
            @eval @inline vadd(v::SVec{$W,$T}, m::Mask{$W,$U}) = vifelse(m, vadd(v, one(v)), v)
            @eval @inline vsub(m::Mask{$W,$U}, v::SVec{$W,$T}) = vsub(vifelse(m, one(v), zero(v)), v)
            @eval @inline vsub(v::SVec{$W,$T}, m::Mask{$W,$U}) = vifelse(m, vsub(v, one(v)), v)
            @eval @inline vmul(m::Mask{$W,$U}, v::SVec{$W,$T}) = vifelse(m, v, vzero(SVec{$W,$T}))
            @eval @inline vmul(v::SVec{$W,$T}, m::Mask{$W,$U}) = vifelse(m, v, vzero(SVec{$W,$T}))
        end
    end
end
@inline vfdiv(v1::_Vec{W,Int32}, v2::_Vec{W,Int32}) where {W} = vfdiv(vconvert(_Vec{W,Float32}, v1), vconvert(_Vec{W,Float32}, v2))
@inline vfdiv(v1::_Vec{W,Int64}, v2::_Vec{W,Int64}) where {W} = vfdiv(vconvert(_Vec{W,Float64}, v1), vconvert(_Vec{W,Float64}, v2))
@inline function Base.:(/)(v1::SVec{W,I},v2::SVec{W,I}) where {W,I<:Integer}
    T = sizeequivalentfloat(I)
    SVec(vfdiv(vconvert(Vec{W,T}, v1), vconvert(Vec{W,T}, v2)))
end
@inline Base.:(/)(v::SVec{W,I}, s::Integer) where {W,I} = v / vbroadcast(SVec{W,I}, s % I)
@inline Base.:(/)(s::Integer, v::SVec{W,I}) where {W,I} = vbroadcast(SVec{W,I}, s % I) / v
@inline vmax(x::Number, y::Number) = Base.FastMath.max_fast(x,y)
@vectordef vmax function Base.max(v1, v2) where {W,T<:FloatingTypes}
    vifelse(vgreater(extract_data(v1), extract_data(v2)), extract_data(v1), extract_data(v2))
end
@inline vmin(x::Number, y::Number) = Base.FastMath.min_fast(x,y)
@vectordef vmin function Base.min(v1, v2) where {W,T<:FloatingTypes}
    vifelse(vless(extract_data(v1), extract_data(v2)), extract_data(v1), extract_data(v2))
end

@inline vflipsign(s1::FloatingTypes, s2::FloatingTypes) = flipsign(s1, s2)

@vectordef vflipsign function Base.flipsign(v1, v2) where {W,T<:FloatingTypes}
    vifelse(vsignbit(extract_data(v2)), vsub(extract_data(v1)), extract_data(v1))
end

@inline vcopysign(v1::_Vec{W,T}, v2::_Vec{W,U}) where {W,T,U<:Unsigned} = vcopysign(v1, vreinterpret(_Vec{W,T}, v2))

for op ∈ (:fma, :muladd)
# let op = :fma
    rename = VECTOR_SYMBOLS[op]
    # if op == :muladd
        # rename = Symbol(:e,rename)
    # end
    @eval begin

        @vectordef $rename function Base.$op(v1, v2, v3) where {W,T<:FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1), extract_data(v2), extract_data(v3))
        end
        @vectordef $rename function Base.$op(s1::T, v2, v3) where {W,T<:FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), vbroadcast(Vec{W,T}, s1), extract_data(v2), extract_data(v3))
        end
        @vectordef $rename function Base.$op(v1, s2::T, v3) where {W,T<:FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1), vbroadcast(Vec{W,T}, s2), extract_data(v3))
        end
        @vectordef $rename function Base.$op(v1, v2, s3::T) where {W,T<:FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1), extract_data(v2), vbroadcast(Vec{W,T}, s3))
        end
        @vectordef $rename function Base.$op(s1::T, s2::T, v3) where {W,T<:FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), vbroadcast(Vec{W,T}, s1), vbroadcast(Vec{W,T}, s2), extract_data(v3))
        end
        @vectordef $rename function Base.$op(s1::T, v2, s3::T) where {W,T<:FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), vbroadcast(Vec{W,T}, s1), extract_data(v2), vbroadcast(Vec{W,T}, s3))
        end
        @vectordef $rename function Base.$op(v1, s2::T, s3::T) where {W,T<:FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1), vbroadcast(Vec{W,T}, s2), vbroadcast(Vec{W,T}, s3))
        end
        @vpromote $rename 3
    end
end

# Type promotion

# Promote scalars of all ScalarTypes to vectors of FloatingTypes, leaving the
# vector type unchanged

for op ∈ (
        :(==), :(!=), :(<), :(<=), :(>), :(>=),
        :+, :-, :*, :/,# :^,
        :copysign, :flipsign, :max, :min, :%
    )
    rename = VECTOR_SYMBOLS[op]
    @eval begin
        @vectordef $rename function Base.$op(s1::ScalarTypes, v2) where {W,T<:FloatingTypes}
            ev2 = extract_data(v2)
            $rename(vbroadcast(typeof(ev2), s1), ev2)
        end
        @vectordef $rename function Base.$op(v1, s2::ScalarTypes) where {W,T<:FloatingTypes}
            ev1 = extract_data(v1)
            $rename(ev1, vbroadcast(typeof(ev1), s2))
        end
        @inline function $rename(s1::T2, v2::_Vec{W,T}) where {W,T<:Integer,T2<:FloatingTypes}
            $rename(vbroadcast(_Vec{W,T2}, s1), v2)
        end
        @inline function $rename(v1::_Vec{W,T}, s2::T2) where {W,T<:Integer,T2<:FloatingTypes}
            $rename(v1, vbroadcast(_Vec{W,T2}, s2))
        end
        @vpromote $rename 2
    end
end
for op ∈ (:(==), :(!=), :(<), :(<=), :(>), :(>=))
    rename = VECTOR_SYMBOLS[op]
    @eval begin
        @inline function Base.$op(s1::T2, v2::SVec{W,T}) where {W,T<:Integer,T2<:FloatingTypes}
            Mask{W}($rename(vbroadcast(Vec{W,T2}, s1), vconvert(Vec{W,T2}, extract_data(v2))))
        end
        @inline function Base.$op(v1::SVec{W,T}, s2::T2) where {W,T<:Integer,T2<:FloatingTypes}
            Mask{W}($rename(vconvert(Vec{W,T2}, extract_data(v1)), vbroadcast(Vec{W,T2}, s2)))
        end
        @inline function $rename(s1::T2, v2::SVec{W,T}) where {W,T<:Integer,T2<:FloatingTypes}
            Mask{W}($rename(vbroadcast(Vec{W,T2}, s1), extract_data(v2)))
        end
        @inline function $rename(v1::SVec{W,T}, s2::T2) where {W,T<:Integer,T2<:FloatingTypes}
            Mask{W}($rename(extract_data(v1), vbroadcast(Vec{W,T2}, s2)))
        end
    end
end
for op ∈ (:+, :-, :*, :/, :copysign, :flipsign, :max, :min, :%)
    rename = VECTOR_SYMBOLS[op]
    @eval begin
        @inline function Base.$op(s1::T2, v2::SVec{W,T}) where {W,T<:Integer,T2<:FloatingTypes}
            SVec($rename(vbroadcast(Vec{W,T2}, s1), vconvert(Vec{W,T2}, extract_data(v2))))
        end
        @inline function Base.$op(v1::SVec{W,T}, s2::T2) where {W,T<:Integer,T2<:FloatingTypes}
            SVec($rename(vconvert(Vec{W,T2}, extract_data(v1)), vbroadcast(Vec{W,T2}, s2)))
        end
        @inline function $rename(s1::T2, v2::SVec{W,T}) where {W,T<:Integer,T2<:FloatingTypes}
            SVec($rename(vbroadcast(Vec{W,T2}, s1), extract_data(v2)))
        end
        @inline function $rename(v1::SVec{W,T}, s2::T2) where {W,T<:Integer,T2<:FloatingTypes}
            SVec($rename(extract_data(v1), vbroadcast(Vec{W,T2}, s2)))
        end
    end
end

for op ∈ (
        :(+), :(-), :(*), :(/), :(%),# :(^),
        :copysign#, :max, :min
    )
    # exact / explicit version
    erename = Symbol(:e, VECTOR_SYMBOLS[op])
    @eval begin

        @evectordef $erename function Base.$op(v1, s2::ScalarTypes) where {W,T <: FloatingTypes}
            $erename(extract_data(v1), vbroadcast(Vec{W,T}, s2))
        end
        @evectordef $erename function Base.$op(s1::ScalarTypes, v2) where {W,T <: FloatingTypes}
            $erename(vbroadcast(Vec{W,T}, s1), extract_data(v2))
        end
    end
end

# Reduction operations

# TODO: map, mapreduce

function getneutral(op::Symbol, ::Type{T}) where {T}
    if op == :&
        return ~zero(T)
    elseif op == :|
        return zero(T)
    elseif op == :max
        return typemin(T)
    elseif op == :min
        return typemax(T)
    elseif op == :+
        return zero(T)
    elseif op == :*
        return one(T)
    end
    throw("Op $op not recognized.")
end

@static if Base.libllvm_version >= v"9"
    @generated function vsum(v::_Vec{_W,T}) where {_W,T<:FloatingTypes}
        W = _W + 1
        instrs = String[]
        typ = llvmtype(T)
        vtyp = "<$W x $typ>"
        bits = 8sizeof(T)
        ins = "@llvm.experimental.vector.reduce.v2.fadd.f$(bits).v$(W)f$(bits)"
        decl = "declare $(typ) $(ins)($(typ), $(vtyp))"
        push!(instrs, "%res = call $(fastflags(T)) $typ $ins($typ 0.0, $vtyp %0)")
        push!(instrs, "ret $typ %res")
        quote
            $(Expr(:meta, :inline))
            llvmcall(
                $((decl, join(instrs, "\n"))),
                $T, Tuple{Vec{$W,$T}}, v
            )
        end
    end
    @generated function vprod(v::_Vec{_W,T}) where {_W,T<:FloatingTypes}
        W = _W + 1
        instrs = String[]
        typ = llvmtype(T)
        vtyp = "<$W x $typ>"
        bits = 8sizeof(T)
        ins = "@llvm.experimental.vector.reduce.v2.fmul.f$(bits).v$(W)f$(bits)"
        decl = "declare $(typ) $(ins)($(typ), $(vtyp))"
        push!(instrs, "%res = call $(fastflags(T)) $typ $ins($typ 1.0, $vtyp %0)")
        push!(instrs, "ret $typ %res")
        quote
            $(Expr(:meta, :inline))
            llvmcall(
                $((decl, join(instrs, "\n"))),
                $T, Tuple{Vec{$W,$T}}, v
            )
        end
    end
    @generated function reduced_add(v::_Vec{_W,T}, s::T) where {_W,T<:FloatingTypes}
        W = _W + 1
        instrs = String[]
        typ = llvmtype(T)
        vtyp = "<$W x $typ>"
        bits = 8sizeof(T)
        ins = "@llvm.experimental.vector.reduce.v2.fadd.f$(bits).v$(W)f$(bits)"
        decl = "declare $(typ) $(ins)($(typ), $(vtyp))"
        push!(instrs, "%res = call $(fastflags(T)) $typ $ins($typ %1, $vtyp %0)")
        push!(instrs, "ret $typ %res")
        quote
            $(Expr(:meta, :inline))
            llvmcall(
                $((decl, join(instrs, "\n"))),
                $T, Tuple{Vec{$W,$T},$T}, v, s
            )
        end
    end
    @generated function reduced_prod(v::_Vec{_W,T}, s::T) where {_W,T<:FloatingTypes}
        W = _W + 1
        instrs = String[]
        typ = llvmtype(T)
        vtyp = "<$W x $typ>"
        bits = 8sizeof(T)
        ins = "@llvm.experimental.vector.reduce.v2.fmul.f$(bits).v$(W)f$(bits)"
        decl = "declare $(typ) $(ins)($(typ), $(vtyp))"
        push!(instrs, "%res = call $(fastflags(T)) $typ $ins($typ %1, $vtyp %0)")
        push!(instrs, "ret $typ %res")
        quote
            $(Expr(:meta, :inline))
            llvmcall(
                $((decl, join(instrs, "\n"))),
                $T, Tuple{Vec{$W,$T},$T}, v, s
            )
        end
    end
    @generated function vsum(v::_Vec{_W,T}) where {_W,T<:IntegerTypes}
        W = _W + 1
        instrs = String[]
        typ = llvmtype(T)
        vtyp = "<$W x $typ>"
        bits = 8sizeof(T)
        ins = "@llvm.experimental.vector.reduce.add.v$(W)i$(bits)"
        decl = "declare $(typ) $(ins)($(vtyp))"
        push!(instrs, "%res = call $typ $ins($vtyp %0)")
        push!(instrs, "ret $typ %res")
        quote
            $(Expr(:meta, :inline))
            llvmcall(
                $((decl, join(instrs, "\n"))),
                $T, Tuple{Vec{$W,$T}}, v
            )
        end
    end
    @generated function vprod(v::_Vec{_W,T}) where {_W,T<:IntegerTypes}
        W = _W + 1
        instrs = String[]
        typ = llvmtype(T)
        vtyp = "<$W x $typ>"
        bits = 8sizeof(T)
        ins = "@llvm.experimental.vector.reduce.mul.v$(W)i$(bits)"
        decl = "declare $(typ) $(ins)($(vtyp))"
        push!(instrs, "%res = call $typ $ins($vtyp %0)")
        push!(instrs, "ret $typ %res")
        quote
            $(Expr(:meta, :inline))
            llvmcall(
                $((decl, join(instrs, "\n"))),
                $T, Tuple{Vec{$W,$T}}, v
            )
        end
    end
else
    @generated function llvmwrapreduce(::Val{Op}, v::_Vec{_W,T}) where {Op,_W,T}
        @assert isa(Op, Symbol)
        W = _W + 1
        z = getneutral(Op, T)
        typ = llvmtype(T)
        decls = String[]
        instrs = String[]
        n = W
        nam = "%0"
        nold,n = n, VectorizationBase.nextpow2(n)
        if n > nold
            namold,nam = nam,"%vec_$n"
            append!(instrs,
                    extendvector(namold, nold, typ, n, n-nold,
                                 llvmtypedconst(T,z), nam))
        end
        while n > 1
            nold,n = n, div(n, 2)
            namold,nam = nam,"%vec_$n"
            vtyp = "<$n x $typ>"
            ins = llvmins(Op, n, T)
            append!(instrs, subvector(namold, nold, typ, "$(nam)_1", n, 0))
            append!(instrs, subvector(namold, nold, typ, "$(nam)_2", n, n))
            if ins[1] == '@'
                push!(decls, "declare $vtyp $ins($vtyp, $vtyp)")
                push!(instrs,
                      "$nam = call $vtyp $ins($vtyp $(nam)_1, $vtyp $(nam)_2)")
            else
                push!(instrs, "$nam = $ins $vtyp $(nam)_1, $(nam)_2")
            end
        end
        push!(instrs, "%res = extractelement <$n x $typ> $nam, i32 0")
        push!(instrs, "ret $typ %res")
        quote
            $(Expr(:meta, :inline))
            llvmcall(
                $((join(decls, "\n"), join(instrs, "\n"))),
                $T, Tuple{Vec{$W,$T}}, v
            )
        end
    end

    for (name, rename, op) ∈ [(:(Base.sum),:vsum,:+), (:(Base.prod),:vprod,:*)]
        @eval begin
            @inline $rename(v::AbstractSIMDVector{W}) where {W} = llvmwrapreduce(Val{$(QuoteNode(op))}(), extract_data(v))
        end
    end

end
@static if Base.libllvm_version >= v"8"

    @generated function vsub(v::_Vec{_W,T}) where {_W,T<:FloatingTypes}
        W = _W + 1
        typ = llvmtype(T)
        vtyp = "<$W x $typ>"
        instrs = "%res = fneg $(fastflags(T)) $vtyp %0\nret $vtyp %res"
        quote
            $(Expr(:meta, :inline))
            llvmcall( $instrs, Vec{$W,$T}, Tuple{Vec{$W,$T}}, v )
        end
    end
    @inline Base.:(-)(v::SVec{W,T}) where {W,T} = SVec(vsub(extract_data(v)))
    @inline vsub(v::SVec{W,T}) where {W,T} = SVec(vsub(extract_data(v)))
else
    let rename = VECTOR_SYMBOLS[:-]
        @eval begin
            @vectordef $rename function Base.:(-)(v1) where {W,T<:FloatingTypes}
                llvmwrap(Val{:(-)}(), extract_data(v1))
            end
        end
    end

end
@inline vsub(x::FloatingTypes) = Base.FastMath.sub_fast(x)

@static if Base.libllvm_version >= v"6"
    @generated function vmaximum(v::_Vec{_W,T}) where {_W,T<:ScalarTypes}
        W = _W + 1
        instrs = String[]
        typ = llvmtype(T)
        vtyp = "<$W x $typ>"
        bits = 8sizeof(T)
        if T <: FloatingTypes
            prefix = prefix2 = 'f'
        else
            prefix2 = 'i'
            prefix = T <: Signed ? 's' : 'u'
        end
        ins = "@llvm.experimental.vector.reduce.$(prefix)max.v$(W)$(prefix2)$(bits)"
        decl = "declare $(typ) $(ins)($(vtyp))"
        push!(instrs, "%res = call $typ $ins($vtyp %0)")
        push!(instrs, "ret $typ %res")
        quote
            $(Expr(:meta, :inline))
            llvmcall(
                $((decl, join(instrs, "\n"))),
                $T, Tuple{Vec{$W,$T}}, v
            )
        end
    end
    @generated function vminimum(v::_Vec{_W,T}) where {_W,T<:ScalarTypes}
        W = _W + 1
        instrs = String[]
        typ = llvmtype(T)
        vtyp = "<$W x $typ>"
        bits = 8sizeof(T)
        if T <: FloatingTypes
            prefix = prefix2 = 'f'
        else
            prefix2 = 'i'
            prefix = T <: Signed ? 's' : 'u'
        end
        ins = "@llvm.experimental.vector.reduce.$(prefix)min.v$(W)$(prefix2)$(bits)"
        decl = "declare $(typ) $(ins)($(vtyp))"
        push!(instrs, "%res = call $typ $ins($vtyp %0)")
        push!(instrs, "ret $typ %res")
        quote
            $(Expr(:meta, :inline))
            llvmcall(
                $((decl, join(instrs, "\n"))),
                $T, Tuple{Vec{$W,$T}}, v
            )
        end
    end
else
    for (name, rename, op) ∈ [(:(Base.maximum), :vmaximum, :max), (:(Base.minimum), :vminimum, :min)]

        @eval begin
            @inline $rename(v::AbstractSIMDVector{W}) where {W} = llvmwrapreduce(Val{$(QuoteNode(op))}(), extract_data(v))
        end
    end

end
for (name, rename, op) ∈ [(:(Base.maximum), :vmaximum, :max), (:(Base.minimum), :vminimum, :min), (:(Base.sum),:vsum,:+), (:(Base.prod),:vprod,:*)]
    @eval begin
        @inline $name(v::SVec{W,T}) where {W,T} = $rename(extract_data(v))
        @inline $rename(v::SVec{W,T}) where {W,T} = $rename(extract_data(v))
        # @inline $name(v::SVec{W,T}) where {W,T} = llvmwrapreduce(Val{$(QuoteNode(op))}(), extract_data(v))
    end
end

# TODO: Handle cases with vectors of different lengths correctly!
@inline vmul(x::Number, y::Number) = Base.FastMath.mul_fast(x, y)
@inline vadd(x::Number, y::Number) = Base.FastMath.add_fast(x, y)
@inline vmul(x,y,z...) = vmul(x,vmul(y,z...))
@inline vadd(x,y,z...) = vadd(x,vadd(y,z...))
@inline vadd(a,b,c,d) = evadd(vadd(a,b), vadd(c,d))
@inline vadd(a,b,c,d,e,f) = evadd(evadd(vadd(a,b), vadd(c,d)),vad(e,f))
@inline vadd(a,b,c,d,e,f,g,h) = vadd(evadd(vadd(a,b), vadd(c,d)),evadd(vadd(e,f),vadd(g,h)))
@inline Base.:(+)(a::SVec, b, c, d) = evadd(vadd(a,b), vadd(c,d))
@inline vmuladd(a::Number, b::Number, c::Number) = muladd(a, b, c)
# These intrinsics are not FastMath.

@generated function vfmadd(v1::_Vec{_W,T}, v2::_Vec{_W,T}, v3::_Vec{_W,T}) where {_W,T<:FloatingTypes}
    W = _W + 1
    typ = llvmtype(T)
    vtyp = "<$W x $typ>"
    # ins = "@llvm.fma.v$(W)f$(8*sizeof(T))"
    ins = "@llvm.fmuladd.v$(W)f$(8*sizeof(T))"
    decl = "declare $vtyp $ins($vtyp, $vtyp, $vtyp)"
    instrs = "%res = call $vtyp $ins($vtyp %0, $vtyp %1, $vtyp %2)\nret $vtyp %res"
    quote
        $(Expr(:meta, :inline))
        llvmcall(
            $((decl,instrs)),
            Vec{$W,$T}, Tuple{Vec{$W,$T}, Vec{$W,$T}, Vec{$W,$T}},
            v1, v2, v3
        )
    end
end
vfmadd(a::Number, b::Number, c::Number) = muladd(a, b, c)
vfnmadd(a::Number, b::Number, c::Number) = muladd(-a, b, c)
vfmsub(a::Number, b::Number, c::Number) = muladd(a, b, -c)
vfnmsub(a::Number, b::Number, c::Number) = -muladd(a, b, c)
@inline vfmadd(a::_Vec{W,T}, b::_Vec{W,T}, c::_Vec{W,T}) where {W,T <: Integer} = vmuladd(a, b, c)
@inline vfnmadd(a::_Vec{W,T}, b::_Vec{W,T}, c::_Vec{W,T}) where {W,T} = vfmadd(vsub(a), b, c)
@inline vfmsub(a::_Vec{W,T}, b::_Vec{W,T}, c::_Vec{W,T}) where {W,T} = vfmadd(a, b, vsub(c))
@inline vfnmsub(a::_Vec{W,T}, b::_Vec{W,T}, c::_Vec{W,T}) where {W,T} = vsub(vfmadd(a, b, c))
@vpromote vfmadd 3
@vpromote vfnmadd 3
@vpromote vfmsub 3
@vpromote vfnmsub 3

@generated function vfmadd_fast(v1::_Vec{_W,T}, v2::_Vec{_W,T}, v3::_Vec{_W,T}) where {_W,T<:FloatingTypes}
    W = _W + 1
    typ = llvmtype(T)
    vtyp = "<$W x $typ>"
    # ins = "@llvm.fmuladd.v$(W)f$(8*sizeof(T))"
    # decls = "declare $vtyp $ins($vtyp, $vtyp, $vtyp)"
    # instrs = "%res = call fast $vtyp $ins($vtyp %0, $vtyp %1, $vtyp %2)\nret $vtyp %res"
    # I don't really want reassoc on the fmul part, so this is what I'm going with.
    instrs = """
    %prod = fmul nnan ninf nsz arcp contract $vtyp %0, %1
    %res = fadd $(fastflags(T)) $vtyp %prod, %2
    ret $vtyp %res
    """
    quote
        $(Expr(:meta, :inline))
        llvmcall(
            $instrs, Vec{$W,$T}, Tuple{Vec{$W,$T}, Vec{$W,$T}, Vec{$W,$T}}, v1, v2, v3
        )
    end
end
vfmadd_fast(a::Number, b::Number, c::Number) = Base.FastMath.add_fast(Base.FastMath.mul_fast(a, b), c)
vfnmadd_fast(a::Number, b::Number, c::Number) = Base.FastMath.sub_fast(c, Base.FastMath.mul_fast(a, b))
vfmsub_fast(a::Number, b::Number, c::Number) = Base.FastMath.sub_fast(Base.FastMath.mul_fast(a, b), c)
vfnmsub_fast(a::Number, b::Number, c::Number) = Base.FastMath.sub_fast(Base.FastMath.add_fast(Base.FastMath.mul_fast(a, b), c))
@inline vfmadd_fast(a::_Vec{W,T}, b::_Vec{W,T}, c::_Vec{W,T}) where {W,T <: Integer} = vmuladd(a, b, c)
@inline vfnmadd_fast(a::_Vec{W,T}, b::_Vec{W,T}, c::_Vec{W,T}) where {W,T} = vfmadd_fast(vsub(a), b, c)
@inline vfmsub_fast(a::_Vec{W,T}, b::_Vec{W,T}, c::_Vec{W,T}) where {W,T} = vfmadd_fast(a, b, vsub(c))
@inline vfnmsub_fast(a::_Vec{W,T}, b::_Vec{W,T}, c::_Vec{W,T}) where {W,T} = vsub(vfmadd_fast(a, b, c))
@inline vfmadd_fast(m::AbstractMask{W}, b::AbstractSIMDVector{W}, c::AbstractSIMDVector{W}) where {W} = vifelse(m, vadd(b, c), c)
@inline vfmadd_fast(b::AbstractSIMDVector{W}, m::AbstractMask{W}, c::AbstractSIMDVector{W}) where {W} = vifelse(m, vadd(b, c), c)
@inline vfnmadd_fast(m::AbstractMask{W}, b::AbstractSIMDVector{W}, c::AbstractSIMDVector{W}) where {W} = vifelse(m, vsub(c, b), c)
@inline vfnmadd_fast(b::AbstractSIMDVector{W}, m::AbstractMask{W}, c::AbstractSIMDVector{W}) where {W} = vifelse(m, vsub(c, b), c)
@inline vfmsub_fast(m::AbstractMask{W}, b::V, c::AbstractSIMDVector{W}) where {W,V<:AbstractSIMDVector{W}} = vsub(vifelse(m, b, vzero(V)), c)
@inline vfmsub_fast(b::V, m::AbstractMask{W}, c::AbstractSIMDVector{W}) where {W,V<:AbstractSIMDVector{W}} = vsub(vifelse(m, b, vzero(V)), c)
@inline vfnmsub_fast(m::AbstractMask{W}, b::AbstractSIMDVector{W}, c::AbstractSIMDVector{W}) where {W} = vsub(vfmadd_fast(m, b, c))
@inline vfnmsub_fast(b::AbstractSIMDVector{W}, m::AbstractMask{W}, c::AbstractSIMDVector{W}) where {W} = vsub(vfmadd_fast(m, b, c))
@inline vfmadd(m::AbstractMask{W}, b::AbstractSIMDVector{W}, c::AbstractSIMDVector{W}) where {W} = vifelse(m, vadd(b, c), c)
@inline vfmadd(b::AbstractSIMDVector{W}, m::AbstractMask{W}, c::AbstractSIMDVector{W}) where {W} = vifelse(m, vadd(b, c), c)
@inline vfnmadd(m::AbstractMask{W}, b::AbstractSIMDVector{W}, c::AbstractSIMDVector{W}) where {W} = vifelse(m, vsub(c, b), c)
@inline vfnmadd(b::AbstractSIMDVector{W}, m::AbstractMask{W}, c::AbstractSIMDVector{W}) where {W} = vifelse(m, vsub(c, b), c)
@inline vfmsub(m::AbstractMask{W}, b::V, c::AbstractSIMDVector{W}) where {W,V<:AbstractSIMDVector{W}} = vsub(vifelse(m, b, vzero(V)), c)
@inline vfmsub(b::V, m::AbstractMask{W}, c::AbstractSIMDVector{W}) where {W,V<:AbstractSIMDVector{W}} = vsub(vifelse(m, b, vzero(V)), c)
@inline vfnmsub(m::AbstractMask{W}, b::AbstractSIMDVector{W}, c::AbstractSIMDVector{W}) where {W} = vsub(vfmadd(m, b, c))
@inline vfnmsub(b::AbstractSIMDVector{W}, m::AbstractMask{W}, c::AbstractSIMDVector{W}) where {W} = vsub(vfmadd(m, b, c))
@inline vfmadd231(m::AbstractMask{W}, b::AbstractSIMDVector{W}, c::AbstractSIMDVector{W}) where {W} = vifelse(m, vadd(b, c), c)
@inline vfmadd231(b::AbstractSIMDVector{W}, m::AbstractMask{W}, c::AbstractSIMDVector{W}) where {W} = vifelse(m, vadd(b, c), c)
@inline vfnmadd231(m::AbstractMask{W}, b::AbstractSIMDVector{W}, c::AbstractSIMDVector{W}) where {W} = vifelse(m, vsub(c, b), c)
@inline vfnmadd231(b::AbstractSIMDVector{W}, m::AbstractMask{W}, c::AbstractSIMDVector{W}) where {W} = vifelse(m, vsub(c, b), c)
@inline vfmsub231(m::AbstractMask{W}, b::V, c::AbstractSIMDVector{W}) where {W,V<:AbstractSIMDVector{W}} = vsub(vifelse(m, b, vzero(V)), c)
@inline vfmsub231(b::V, m::AbstractMask{W}, c::AbstractSIMDVector{W}) where {W,V<:AbstractSIMDVector{W}} = vsub(vifelse(m, b, vzero(V)), c)
@inline vfnmsub231(m::AbstractMask{W}, b::AbstractSIMDVector{W}, c::AbstractSIMDVector{W}) where {W} = vsub(vfmadd231(m, b, c))
@inline vfnmsub231(b::AbstractSIMDVector{W}, m::AbstractMask{W}, c::AbstractSIMDVector{W}) where {W} = vsub(vfmadd231(m, b, c))
@vpromote vfmadd_fast 3
@vpromote vfnmadd_fast 3
@vpromote vfmsub_fast 3
@vpromote vfnmsub_fast 3

@inline Base.:(+)(m::AbstractMask{W}, v::AbstractStructVec{W}) where {W} = vifelse(m, vadd(v, one(v)), v)
@inline Base.:(+)(v::AbstractStructVec{W}, m::AbstractMask{W}) where {W} = vifelse(m, vadd(v, one(v)), v)
@inline vadd(m::AbstractMask{W}, v::AbstractSIMDVector{W}) where {W} = vifelse(m, vadd(v, one(v)), v)
@inline vadd(v::AbstractSIMDVector{W}, m::AbstractMask{W}) where {W} = vifelse(m, vadd(v, one(v)), v)
@inline Base.:(-)(m::AbstractMask{W}, v::AbstractStructVec{W}) where {W} = vsub(vifelse(m, one(v), zero(v)), v)
@inline Base.:(-)(v::AbstractStructVec{W}, m::AbstractMask{W}) where {W} = vifelse(m, vsub(v, one(v)), v)
@inline vsub(m::AbstractMask{W}, v::AbstractSIMDVector{W}) where {W} = vsub(vifelse(m, one(v), zero(v)), v)
@inline vsub(v::AbstractSIMDVector{W}, m::AbstractMask{W}) where {W} = vifelse(m, vsub(v, one(v)), v)
@inline Base.:(*)(v::AbstractStructVec{W,T}, m::Mask{W}) where {W,T} = SVec(vifelse(m, extract_data(v), vzero(Vec{W,T})))
@inline Base.:(*)(m::Mask{W}, v::AbstractStructVec{W,T}) where {W,T} = SVec(vifelse(m, extract_data(v), vzero(Vec{W,T})))
@inline vmul(v::AbstractStructVec{W,T}, m::Mask{W}) where {W,T} = SVec(vifelse(m, extract_data(v), vzero(Vec{W,T})))
@inline vmul(m::Mask{W}, v::AbstractStructVec{W,T}) where {W,T} = SVec(vifelse(m, extract_data(v), vzero(Vec{W,T})))
@inline Base.:(*)(v::AbstractStructVec{W,T}, m::SVec{W,Bool}) where {W,T} = SVec(vifelse(extract_data(m), extract_data(v), vzero(Vec{W,T})))
@inline Base.:(*)(m::SVec{W,Bool}, v::AbstractStructVec{W,T}) where {W,T} = SVec(vifelse(extract_data(m), extract_data(v), vzero(Vec{W,T})))
@inline vmul(v::AbstractStructVec{W,T}, m::SVec{W,Bool}) where {W,T} = SVec(vifelse(extract_data(m), extract_data(v), vzero(Vec{W,T})))
@inline vmul(m::SVec{W,Bool}, v::AbstractStructVec{W,T}) where {W,T} = SVec(vifelse(extract_data(m), extract_data(v), vzero(Vec{W,T})))

# Lowers to same split mul-add llvm as the
# @inline vfmadd(a, b, c) = vadd(vmul( a, b), c)
# definition, so I wont bother implementing
# vfnmadd, vfmsub, and vfnmsub
# in this manner.

# @inline Base.:*(a::IntegerTypes, b::SVec{W,T}) where {W,T} = SVec{W,T}(a) * b
# @inline Base.:*(a::T, b::SVec{W,<:IntegerTypes}) where {W,T<:FloatingTypes} = SVec(vmul(vbroadcast(Vec{W,T}, a), vconvert(Vec{W,T}, extract_data(b))))

# const Vec{W,T} = NTuple{W,Core.VecElement{T}}
@inline vfmadd231(a::_Vec, b, c) = vfmadd(a, b, c)
@inline vfnmadd231(a::_Vec, b, c) = vfnmadd(a, b, c)
@inline vfmsub231(a::_Vec, b, c) = vfmsub(a, b, c)
@inline vfnmsub231(a::_Vec, b, c) = vfnmsub(a, b, c)
@inline vfmadd231(a::Number, b::Number, c::Number) = Base.FastMath.add_fast(Base.FastMath.mul_fast(a, b), c)
@inline vfnmadd231(a::Number, b::Number, c::Number) = Base.FastMath.sub_fast(c, Base.FastMath.mul_fast(a, b))
@inline vfmsub231(a::Number, b::Number, c::Number) = Base.FastMath.sub_fast(Base.FastMath.mul_fast(a, b), c)
@inline vfnmsub231(a::Number, b::Number, c::Number) = Base.FastMath.sub_fast(Base.FastMath.sub_fast(c), Base.FastMath.mul_fast(a, b))
@inline function vifelse(::typeof(vfmadd231), m::AbstractMask, a, b, c)
    vifelse(m, vfmadd(a, b, c), c)
end
@inline function vifelse(::typeof(vfnmadd231), m::AbstractMask, a, b, c)
    vifelse(m, vfnmadd(a, b, c), c)
end
@inline function vifelse(::typeof(vfmsub231), m::AbstractMask, a, b, c)
    vifelse(m, vfmsub(a, b, c), c)
end
@inline function vifelse(::typeof(vfnmsub231), m::AbstractMask, a, b, c)
    vifelse(m, vfnmsub(a, b, c), c)
end
@inline function vifelse(::typeof(vfmadd_fast), m::AbstractMask, a, b, c)
    vifelse(m, vfmadd(a, b, c), c)
end
@inline function vifelse(::typeof(vfnmadd_fast), m::AbstractMask, a, b, c)
    vifelse(m, vfnmadd(a, b, c), c)
end
@inline function vifelse(::typeof(vfmsub_fast), m::AbstractMask, a, b, c)
    vifelse(m, vfmsub(a, b, c), c)
end
@inline function vifelse(::typeof(vfnmsub_fast), m::AbstractMask, a, b, c)
    vifelse(m, vfnmsub(a, b, c), c)
end

@vpromote vfmadd231 3
@vpromote vfnmadd231 3
@vpromote vfmsub231 3
@vpromote vfnmsub231 3

if VectorizationBase.FMA3
    for T ∈ [Float32,Float64]
        W = 16 ÷ sizeof(T)
        local suffix = T == Float32 ? "ps" : "pd"
        typ = llvmtype(T)
        while W <= VectorizationBase.REGISTER_SIZE ÷ sizeof(T)
            vfmadd_str = """%res = call <$W x $(typ)> asm "vfmadd231$(suffix) \$3, \$2, \$1", "=v,0,v,v"(<$W x $(typ)> %2, <$W x $(typ)> %1, <$W x $(typ)> %0)
                ret <$W x $(typ)> %res"""
            vfnmadd_str = """%res = call <$W x $(typ)> asm "vfnmadd231$(suffix) \$3, \$2, \$1", "=v,0,v,v"(<$W x $(typ)> %2, <$W x $(typ)> %1, <$W x $(typ)> %0)
                ret <$W x $(typ)> %res"""
            vfmsub_str = """%res = call <$W x $(typ)> asm "vfmsub231$(suffix) \$3, \$2, \$1", "=v,0,v,v"(<$W x $(typ)> %2, <$W x $(typ)> %1, <$W x $(typ)> %0)
                ret <$W x $(typ)> %res"""
            vfnmsub_str = """%res = call <$W x $(typ)> asm "vfnmsub231$(suffix) \$3, \$2, \$1", "=v,0,v,v"(<$W x $(typ)> %2, <$W x $(typ)> %1, <$W x $(typ)> %0)
                ret <$W x $(typ)> %res"""
            @eval begin
                @inline function vfmadd231(a::Vec{$W,$T}, b::Vec{$W,$T}, c::Vec{$W,$T})
                    llvmcall($vfmadd_str, Vec{$W,$T}, Tuple{Vec{$W,$T},Vec{$W,$T},Vec{$W,$T}}, a, b, c)
                end
                @inline function vfnmadd231(a::Vec{$W,$T}, b::Vec{$W,$T}, c::Vec{$W,$T})
                    llvmcall($vfnmadd_str, Vec{$W,$T}, Tuple{Vec{$W,$T},Vec{$W,$T},Vec{$W,$T}}, a, b, c)
                end
                @inline function vfmsub231(a::Vec{$W,$T}, b::Vec{$W,$T}, c::Vec{$W,$T})
                    llvmcall($vfmsub_str, Vec{$W,$T}, Tuple{Vec{$W,$T},Vec{$W,$T},Vec{$W,$T}}, a, b, c)
                end
                @inline function vfnmsub231(a::Vec{$W,$T}, b::Vec{$W,$T}, c::Vec{$W,$T})
                    llvmcall($vfnmsub_str, Vec{$W,$T}, Tuple{Vec{$W,$T},Vec{$W,$T},Vec{$W,$T}}, a, b, c)
                end
            end
            if VectorizationBase.AVX512BW && W ≥ 8
                vfmaddmask_str = """%res = call <$W x $(typ)> asm "vfmadd231$(suffix) \$3, \$2, \$1 {\$4}", "=v,0,v,v,^Yk"(<$W x $(typ)> %2, <$W x $(typ)> %1, <$W x $(typ)> %0, i$W %3)
                    ret <$W x $(typ)> %res"""
                vfnmaddmask_str = """%res = call <$W x $(typ)> asm "vfnmadd231$(suffix) \$3, \$2, \$1 {\$4}", "=v,0,v,v,^Yk"(<$W x $(typ)> %2, <$W x $(typ)> %1, <$W x $(typ)> %0, i$W %3)
                    ret <$W x $(typ)> %res"""
                vfmsubmask_str = """%res = call <$W x $(typ)> asm "vfmsub231$(suffix) \$3, \$2, \$1 {\$4}", "=v,0,v,v,^Yk"(<$W x $(typ)> %2, <$W x $(typ)> %1, <$W x $(typ)> %0, i$W %3)
                    ret <$W x $(typ)> %res"""
                vfnmsubmask_str = """%res = call <$W x $(typ)> asm "vfnmsub231$(suffix) \$3, \$2, \$1 {\$4}", "=v,0,v,v,^Yk"(<$W x $(typ)> %2, <$W x $(typ)> %1, <$W x $(typ)> %0, i$W %3)
                    ret <$W x $(typ)> %res"""
                U = VectorizationBase.mask_type(W)
                @eval begin
                    @inline function vifelse(::typeof(vfmadd231), m::Mask{$W,$U}, a::SVec{$W,$T}, b::SVec{$W,$T}, c::SVec{$W,$T})
                        SVec(llvmcall($vfmaddmask_str, Vec{$W,$T}, Tuple{Vec{$W,$T},Vec{$W,$T},Vec{$W,$T},$U}, extract_data(a), extract_data(b), extract_data(c), extract_data(m)))
                    end
                    @inline function vifelse(::typeof(vfnmadd231), m::Mask{$W,$U}, a::SVec{$W,$T}, b::SVec{$W,$T}, c::SVec{$W,$T})
                        SVec(llvmcall($vfnmaddmask_str, Vec{$W,$T}, Tuple{Vec{$W,$T},Vec{$W,$T},Vec{$W,$T},$U}, extract_data(a), extract_data(b), extract_data(c), extract_data(m)))
                    end
                    @inline function vifelse(::typeof(vfmsub231), m::Mask{$W,$U}, a::SVec{$W,$T}, b::SVec{$W,$T}, c::SVec{$W,$T})
                        SVec(llvmcall($vfmsubmask_str, Vec{$W,$T}, Tuple{Vec{$W,$T},Vec{$W,$T},Vec{$W,$T},$U}, extract_data(a), extract_data(b), extract_data(c), extract_data(m)))
                    end
                    @inline function vifelse(::typeof(vfnmsub231), m::Mask{$W,$U}, a::SVec{$W,$T}, b::SVec{$W,$T}, c::SVec{$W,$T})
                        SVec(llvmcall($vfnmsubmask_str, Vec{$W,$T}, Tuple{Vec{$W,$T},Vec{$W,$T},Vec{$W,$T},$U}, extract_data(a), extract_data(b), extract_data(c), extract_data(m)))
                    end
                end
            end
            W += W
        end
    end
end

# @generated function rsqrt_fast(x::NTuple{16,Core.VecElement{Float32}})
#     if VectorizationBase.REGISTER_SIZE == 64
#         return quote
#             $(Expr(:meta,:inline))
#             llvmcall(
#             """ %rs = call <16 x float> asm "vrsqrt14ps \$1, \$0", "=x,x"(<16 x float> %0)
#                 ret <16 x float> %rs""",
#             NTuple{16,Core.VecElement{Float32}}, Tuple{NTuple{16,Core.VecElement{Float32}}}, x)
#         end
#     else
#         return quote
#             vinv(vsqrt(x))
#         end
#     end
# end
# @inline function rsqrt(x::NTuple{16,Core.VecElement{Float32}})
#     r = rsqrt_fast(x)
#     # Performs a Newton step to increase accuracy.
#     # ns = vmuladd(vmul(-0.5f0, x), vmul(r, r), 1.5f0)
#     # vmul(r, ns)
#     ns = vfma(vmul(r,r), x, -3.0f0)
#     vmul(vmul(-0.5f0, r), ns)
# end
@inline rsqrt_fast(v) = vinv(vsqrt(v))
@inline rsqrt_fast(x::SVec) = SVec(rsqrt_fast(extract_data(x)))
@inline rsqrt(x::SVec) = SVec(rsqrt(extract_data(x)))
@inline rsqrt(x) = vinv(vsqrt(x))
@inline vinv(x::AbstractStructVec) = SVec(vinv(extract_data(x)))
@inline vinv(x::IntegerTypes) = vinv(float(x))
@inline vinv(x::_Vec{W,I}) where {W, I <: Union{Int64,UInt64}} = evfdiv(vone(_Vec{W,Float64}), vconvert(_Vec{W,Float64}, x))
@inline vinv(x::AbstractStructVec{W,I}) where {W, I <: Union{Int64,UInt64}} = evfdiv(vone(SVec{W,Float64}), vconvert(SVec{W,Float64}, x))
@inline vinv(x::_Vec{W,I}) where {W, I <: Union{Int32,UInt32}} = evfdiv(vone(_Vec{W,Float32}), vconvert(_Vec{W,Float32}, x))
@inline vinv(x::AbstractStructVec{W,I}) where {W, I <: Union{Int32,UInt32}} = evfdiv(vone(SVec{W,Float32}), vconvert(SVec{W,Float32}, x))
@inline vinv(x::_Vec{W,I}) where {W, I <: Union{Int16,UInt16}} = evfdiv(vone(_Vec{W,Float32}), vconvert(_Vec{W,Float32}, x))
@inline vinv(x::AbstractStructVec{W,I}) where {W, I <: Union{Int16,UInt16}} = evfdiv(vone(SVec{W,Float32}), vconvert(SVec{W,Float32}, x))
@inline Base.inv(v::AbstractStructVec{W,I}) where {W, I<: IntegerTypes} = vinv(v)
@inline Base.sqrt(v::AbstractStructVec{W,I}) where {W, I<: IntegerTypes} = SVec(vsqrt(extract_data(float(v))))

# for accumulating vector results of different sizes.
@generated function vadd(v1::_Vec{_W1,T}, v2::_Vec{_W2,T}) where {_W1,_W2,T}
    W1 = _W1 + 1; W2 = _W2 + 1
    @assert ispow2(W1)
    @assert ispow2(W2)
    W3 = W1
    W4 = W2
    v1e = :v1
    while W3 < W2
        W3 <<= 1
        v1e = Expr(:call, :zeropad, v1e)
    end
    v2e = :v2
    while W4 < W1
        W4 <<= 1
        v2e = Expr(:call, :zeropad, v2e)
    end
    Expr(:block, Expr(:meta, :inline), Expr(:call, :vadd, v1e, v2e))
end
@inline Base.:(+)(x::AbstractStructVec{W1,T}, y::AbstractStructVec{W2,T}) where {W1,W2,T} = SVec(vadd(extract_data(x), extract_data(y)))

@inline Base.abs2(v::SVec) = vmul(v,v)
@inline vabs2(v) = vmul(v,v)
@inline vsum(s::FloatingTypes) = s
@inline vprod(s::FloatingTypes) = s

@inline reduced_add(v::_Vec{W,T}, s::T) where {W,T<:IntegerTypes} = Base.FastMath.add_fast(s, vsum(v))
@inline reduced_add(v::AbstractStructVec{W,T}, s::T) where {W,T<:IntegerTypes} = Base.FastMath.add_fast(s, vsum(v))
@inline reduced_add(s::T, v::_Vec{W,T}) where {W,T<:IntegerTypes} = Base.FastMath.add_fast(s, vsum(v))
@inline reduced_add(s::T, v::AbstractStructVec{W,T}) where {W,T<:IntegerTypes} = Base.FastMath.add_fast(s, vsum(v))

@inline reduced_add(v::AbstractStructVec{W,T}, s::T) where {W,T<:FloatingTypes} = reduced_add(extract_data(v), s)
@inline reduced_add(s::T, v::_Vec{W,T}) where {W,T<:FloatingTypes} = reduced_add(v, s)
@inline reduced_add(s::T, v::AbstractStructVec{W,T}) where {W,T<:FloatingTypes} = reduced_add(extract_data(v), s)

# @inline reduced_add(v1::AbstractSIMDVector{W,T}, v2::AbstractSIMDVector{W,T}) where {W,T} = vadd(v1, v2)
@inline reduced_prod(v::_Vec{W,T}, s::T) where {W,T<:IntegerTypes} = Base.FastMath.mul_fast(s, vprod(v))
@inline reduced_prod(s::T, v::_Vec{W,T}) where {W,T<:IntegerTypes} = Base.FastMath.mul_fast(s, vprod(v))
@inline reduced_prod(v::AbstractStructVec{W,T}, s::T) where {W,T<:IntegerTypes} = Base.FastMath.mul_fast(s, vprod(v))
@inline reduced_prod(s::T, v::AbstractStructVec{W,T}) where {W,T<:IntegerTypes} = Base.FastMath.mul_fast(s, vprod(v))

@inline reduced_prod(s::T, v::_Vec{W,T}) where {W,T<:FloatingTypes} = reduced_prod(v, s)
@inline reduced_prod(v::AbstractStructVec{W,T}, s::T) where {W,T<:FloatingTypes} = reduced_prod(extract_data(v), s)
@inline reduced_prod(s::T, v::AbstractStructVec{W,T}) where {W,T<:FloatingTypes} = reduced_prod(extract_data(v), s)
# @inline reduced_prod(v1::AbstractSIMDVector{W,T}, v2::AbstractSIMDVector{W,T}) where {W,T} = vmul(v1, v2)

@inline reduced_add(v::SVec{W,T}, s::T) where {W,T<:IntegerTypes} = reduced_add(extract_data(v), s)
@inline reduced_prod(v::SVec{W,T}, s::T) where {W,T<:IntegerTypes} = reduced_prod(extract_data(v), s)
@inline reduced_add(v::SVec{W,T}, s::T) where {W,T<:FloatingTypes} = reduced_add(extract_data(v), s)
@inline reduced_prod(v::SVec{W,T}, s::T) where {W,T<:FloatingTypes} = reduced_prod(extract_data(v), s)

@inline reduced_add(v::AbstractSIMDVector{W}, s::T2) where {W,T2} = Base.FastMath.add_fast(s, convert(T2,vsum(v)))
@inline reduced_add(s::T2, v::AbstractSIMDVector{W}) where {W,T2} = Base.FastMath.add_fast(s, convert(T2,vsum(v)))
@inline reduced_add(v1::AbstractSIMDVector{W}, v2::AbstractSIMDVector{W}) where {W} = vadd(v1, v2)
@inline reduced_add(v1::T, v2::T) where {T<:Number} = Base.FastMath.add_fast(v1,v2)
# @inline function reduced_add(v1::T, v2::T) where {T<:Number}
#     @show v1, v2 typeof(v1), typeof(v2)
#     vadd(v1,v2)
# end
@inline reduced_prod(v::AbstractSIMDVector{W}, s::T2) where {W,T2} = Base.FastMath.mul_fast(s, convert(T2,vprod(v)))
@inline reduced_prod(s::T2, v::AbstractSIMDVector{W}) where {W,T2} = Base.FastMath.mul_fast(s, convert(T2,vprod(v)))
@inline reduced_prod(v1::AbstractSIMDVector{W}, v2::AbstractSIMDVector{W}) where {W} = vmul(v1, v2)
@inline reduced_prod(v1::T, v2::T) where {T<:Number} = Base.FastMath.mul_fast(v1,v2)

@inline reduce_to_add(v::_Vec{W,T}, ::T) where {W,T} = vsum(v)
@inline reduce_to_prod(v::_Vec{W,T}, ::T) where {W,T} = vprod(v)
@inline reduce_to_add(v::AbstractStructVec{W,T}, ::T) where {W,T} = vsum(v)
@inline reduce_to_prod(v::AbstractStructVec{W,T}, ::T) where {W,T} = vprod(v)
@inline reduce_to_add(v::AbstractSIMDVector{W}, ::T2) where {W,T2} = convert(T2,vsum(v))
@inline reduce_to_add(v::T, ::T) where {T<:Number} = v
@inline reduce_to_prod(v::AbstractSIMDVector{W}, ::T2) where {W,T2} = convert(T2,vprod(v))
@inline reduce_to_add(v::AbstractSIMDVector{W}, ::AbstractSIMDVector{W}) where {W} = v
@inline reduce_to_prod(v::AbstractSIMDVector{W}, ::AbstractSIMDVector{W}) where {W} = v
@inline reduce_to_prod(s::T, ::T) where {T<:Number} = v

@inline reduced_max(v::_Vec{W,T}, s::T) where {W,T} = max(vmaximum(v), s)
@inline reduced_min(v::_Vec{W,T}, s::T) where {W,T} = min(vminimum(v), s)
@inline reduced_max(v1::_Vec{W,T}, v2::_Vec{W,T}) where {W,T} = vmax(v1, v2)
@inline reduced_min(v1::_Vec{W,T}, v2::_Vec{W,T}) where {W,T} = vmin(v1, v2)
@inline reduce_to_max(v::_Vec{W,T}, s::T) where {W,T} = vmaximum(v)
@inline reduce_to_min(v::_Vec{W,T}, s::T) where {W,T} = vminimum(v)
@inline reduce_to_max(v::_Vec{W,T}, ::_Vec{W,T}) where {W,T} = v
@inline reduce_to_min(v::_Vec{W,T}, ::_Vec{W,T}) where {W,T} = v

@inline reduced_max(v::AbstractStructVec{W,T}, s::T) where {W,T} = max(vmaximum(v), s)
@inline reduced_min(v::AbstractStructVec{W,T}, s::T) where {W,T} = min(vminimum(v), s)
@inline reduced_max(v1::AbstractStructVec{W,T}, v2::AbstractStructVec{W,T}) where {W,T} = vmax(v1, v2)
@inline reduced_min(v1::AbstractStructVec{W,T}, v2::AbstractStructVec{W,T}) where {W,T} = vmin(v1, v2)
@inline reduce_to_max(v::AbstractStructVec{W,T}, s::T) where {W,T} = vmaximum(v)
@inline reduce_to_min(v::AbstractStructVec{W,T}, s::T) where {W,T} = vminimum(v)
@inline reduce_to_max(v::AbstractStructVec{W,T}, ::AbstractStructVec{W,T}) where {W,T} = v
@inline reduce_to_min(v::AbstractStructVec{W,T}, ::AbstractStructVec{W,T}) where {W,T} = v

@inline reduced_max(v::AbstractSIMDVector{W}, s::T2) where {W,T2} = max(convert(T2,vmaximum(v)), s)
@inline reduced_min(v::AbstractSIMDVector{W}, s::T2) where {W,T2} = min(convert(T2,vminimum(v)), s)
@inline reduced_max(v1::AbstractSIMDVector{W}, v2::AbstractStructVec{W,T2}) where {W,T2} = vmax(vconvert(SVec{W,T2},v1), v2)
@inline reduced_min(v1::AbstractSIMDVector{W,}, v2::AbstractStructVec{W,T2}) where {W,T2} = vmin(vconvert(SVec{W,T2},v1), v2)
@inline reduce_to_max(v::AbstractSIMDVector{W}, s::T2) where {W,T2} = convert(T2,vmaximum(v))
@inline reduce_to_min(v::AbstractSIMDVector{W}, s::T2) where {W,T2} = convert(T2,vminimum(v))
@inline reduce_to_max(v::AbstractSIMDVector{W}, ::AbstractStructVec{W,T2}) where {W,T2} = vconvert(SVec{W,T2}, v)
@inline reduce_to_min(v::AbstractSIMDVector{W}, ::AbstractStructVec{W,T2}) where {W,T2} = vconvert(SVec{W,T2}, v)

@inline vnmul(x,y) = vsub(vmul(x, y))
@inline vnsub(x,y) = vsub(y, x)
@inline vmul2(x) = vadd(x,x)
@inline vmul3(x::T) where {T <: Number} = Base.FastMath.mul_fast(T(3), x)
@inline vmul3(v::V) where {W, V <: AbstractSIMDVector{W}} = vmul(vbroadcast(V, T(3)), v)
@inline vadd1(x::T) where {T <: Number} = Base.FastMath.add_fast(x, one(x))
@inline vadd1(v::V) where {W, V <: AbstractSIMDVector{W}} = vadd(vone(V), v)

@generated function addscalar(v::_Vec{_W,T}, s::T) where {_W, T <: Integer}
    W = _W + 1
    typ = "i$(8sizeof(T))"
    vtyp = "<$W x $typ>"
    instrs = String[]
    push!(instrs, "%ie = insertelement $vtyp zeroinitializer, $typ %1, i32 0")
    push!(instrs, "%v = add $vtyp %0, %ie")
    push!(instrs, "ret $vtyp %v")
    quote
        $(Expr(:meta,:inline))
        llvmcall( $(join(instrs,"\n")), NTuple{$W,Core.VecElement{$T}}, Tuple{NTuple{$W,Core.VecElement{$T}},$T}, v, s )
    end
end
@generated function addscalar(v::_Vec{_W,T}, s::T) where {_W, T <: Union{Float16,Float32,Float64}}
    W = _W + 1
    typ = llvmtype(T)
    vtyp = "<$W x $typ>"
    instrs = String[]
    push!(instrs, "%ie = insertelement $vtyp zeroinitializer, $typ %1, i32 0")
    push!(instrs, "%v = fadd $(fastflags(T)) $vtyp %0, %ie")
    push!(instrs, "ret $vtyp %v")
    quote
        $(Expr(:meta,:inline))
        llvmcall( $(join(instrs,"\n")), NTuple{$W,Core.VecElement{$T}}, Tuple{NTuple{$W,Core.VecElement{$T}},$T}, v, s )
    end
end
@inline addscalar(v::SVec, s) = SVec(addscalar(extract_data(v), s))
@inline addscalar(s::T, v::_Vec{W,T}) where {W,T} = addscalar(v, s)
@inline addscalar(s::T, v::AbstractStructVec{W,T}) where {W,T} = addscalar(v, s)
@inline addscalar(a, b) = vadd(a, b)

@generated function mulscalar(v::_Vec{_W,T}, s::T) where {_W, T <: Integer}
    W = _W + 1
    typ = "i$(8sizeof(T))"
    vtyp = "<$W x $typ>"
    instrs = String[]
    push!(instrs, "%ie = insertelement $vtyp $(llvmconst(W, T, 1)), $typ %1, i32 0")
    push!(instrs, "%v = mul $vtyp %0, %ie")
    push!(instrs, "ret $vtyp %v")
    quote
        $(Expr(:meta,:inline))
        llvmcall( $(join(instrs,"\n")), NTuple{$W,Core.VecElement{$T}}, Tuple{NTuple{$W,Core.VecElement{$T}},$T}, v, s )
    end
end
@generated function mulscalar(v::_Vec{_W,T}, s::T) where {_W, T <: Union{Float16,Float32,Float64}}
    W = _W + 1
    typ = llvmtype(T)
    vtyp = "<$W x $typ>"
    instrs = String[]
    push!(instrs, "%ie = insertelement $vtyp $(llvmconst(W, T, 1.0)), $typ %1, i32 0")
    push!(instrs, "%v = fmul $(fastflags(T)) $vtyp %0, %ie")
    push!(instrs, "ret $vtyp %v")
    quote
        $(Expr(:meta,:inline))
        llvmcall( $(join(instrs,"\n")), NTuple{$W,Core.VecElement{$T}}, Tuple{NTuple{$W,Core.VecElement{$T}},$T}, v, s )
    end
end
@inline mulscalar(v::SVec, s) = SVec(mulscalar(extract_data(v), s))
@inline mulscalar(s::T, v::_Vec{W,T}) where {W,T} = mulscalar(v, s)
@inline mulscalar(s::T, v::AbstractStructVec{W,T}) where {W,T} = mulscalar(v, s)
@inline mulscalar(a, b) = vmul(a, b)

function scalar_maxmin(W, ::Type{T}, ismax) where {T}
    comp = if T <: Signed
        LLVM_INS_Int[ismax ? :(>) : :(<)]
    elseif T <: Unsigned
        LLVM_INS_UInt[ismax ? :(>) : :(<)]
    elseif ismax
        "fcmp ogt"
    else
        "fcmp olt"
    end
    basevalue = if T <: Integer
        minmaxzero = ifelse(ismax, typemin(T), typemax(T))
        llvmconst(W, T, minmaxzero)
    else
        llvmconst(W, T, repr(reinterpret(UInt64, ifelse(ismax, -Inf, Inf))))
    end
    typ = llvmtype(T)
    vtyp = "<$W x $typ>"
    instrs = String[]
    push!(instrs, "%ie = insertelement $vtyp $(basevalue), $typ %1, i32 0")
    push!(instrs, "%selection = $comp $vtyp %0, %ie")
    push!(instrs, "%v = select <$W x i1> %selection, $vtyp %0, $vtyp %ie")
    push!(instrs, "ret $vtyp %v")
    instrs
end
@generated function maxscalar(v::_Vec{_W,T}, s::T) where {_W, T}
    W = _W + 1
    instrs = scalar_maxmin(W, T, true)
    quote
        $(Expr(:meta,:inline))
        llvmcall( $(join(instrs,"\n")), NTuple{$W,Core.VecElement{$T}}, Tuple{NTuple{$W,Core.VecElement{$T}},$T}, v, s )
    end
end
@generated function minscalar(v::_Vec{_W,T}, s::T) where {_W, T}
    W = _W + 1
    instrs = scalar_maxmin(W, T, false)
    quote
        $(Expr(:meta,:inline))
        llvmcall( $(join(instrs,"\n")), NTuple{$W,Core.VecElement{$T}}, Tuple{NTuple{$W,Core.VecElement{$T}},$T}, v, s )
    end
end
@inline maxscalar(v::SVec, s) = SVec(maxscalar(extract_data(v), s))
@inline maxscalar(s::T, v::_Vec{W,T}) where {W,T} = maxscalar(v, s)
@inline maxscalar(s::T, v::AbstractStructVec{W,T}) where {W,T} = maxscalar(v, s)
@inline maxscalar(a, b) = vmax(a, b)
@inline minscalar(v::SVec, s) = SVec(minscalar(extract_data(v), s))
@inline minscalar(s::T, v::_Vec{W,T}) where {W,T} = minscalar(v, s)
@inline minscalar(s::T, v::AbstractStructVec{W,T}) where {W,T} = minscalar(v, s)
@inline minscalar(a, b) = vmin(a, b)

@inline Base.literal_pow(::typeof(^), x::AbstractStructVec, ::Val{-2}) = (xi = vinv(x); vmul(xi, xi))
@inline Base.literal_pow(::typeof(^), x::AbstractStructVec, ::Val{-1}) = vinv(x)
@inline Base.literal_pow(::typeof(^), x::V, ::Val{0}) where {V <: AbstractStructVec} = vone(V)
@inline Base.literal_pow(::typeof(^), x::AbstractStructVec, ::Val{1}) = x
@inline Base.literal_pow(::typeof(^), x::AbstractStructVec, ::Val{2}) = vmul(x, x)
@inline Base.literal_pow(::typeof(^), x::AbstractStructVec, ::Val{3}) = vmul(x, vmul(x, x))
@inline function Base.literal_pow(::typeof(^), x::AbstractStructVec, ::Val{4})
    x2 = vmul(x, x)
    vmul(x2, x2)
end
@inline function Base.literal_pow(::typeof(^), x::AbstractStructVec, ::Val{5})
    x2 = vmul(x, x)
    vmul(x, vmul(x2, x2))
end
@inline function Base.literal_pow(::typeof(^), x::AbstractStructVec, ::Val{6})
    x2 = vmul(x, x)
    vmul(x2, vmul(x2, x2))
end
@inline function Base.literal_pow(::typeof(^), x::AbstractStructVec, ::Val{7})
    x2 = vmul(x, x)
    x3 = vmul(x, x2)
    x4 = vmul(x2, x2)
    vmul(x3, x4)
end
@inline function Base.literal_pow(::typeof(^), x::AbstractStructVec, ::Val{8})
    x2 = vmul(x, x)
    x4 = vmul(x2, x2)
    vmul(x4, x4)
end
Base.literal_pow(::typeof(^), x::AbstractStructVec, ::Val{P}) where {P} = x ^ P
# @inline literal_power(x, ::Val{P}) = Base.literal_pow(Base.^, x, Val{P}())

@inline relu(x) = max(x,zero(x))

# These don't work yet.
# @static if Base.libllvm_version ≥ v"10.0.0"
#     @generated function vcolumnwiseload(ptr::Ptr{T}, stride::Integer, ::Val{M}, ::Val{N}) where {T,M,N}
#         typ = llvmtype(T)
#         W = M * N
#         vtyp = "<$W x $typ>"
#         ptyp = "i$(8sizeof(Int))"
#         # suffix = 'r' * M * 'c' * N * typ
#         suffix = "v$(W)$(T <: Integer ? 'i' : 'f')$(8sizeof(T))"
#         instr = "@llvm.matrix.columnwise.load.$suffix"
#         decl = "declare $vtyp $instr($vtyp*, i32, i32, i32)"
#         instrs = String[]
#         push!(instrs, "%ptr = inttoptr $ptyp %0 to $vtyp*")
#         push!(instrs, "%res = call $vtyp $instr($vtyp* %ptr, i32 %1, i32 $M, i32 $N)")
#         push!(instrs, "ret $vtyp %res")
#         quote
#             $(Expr(:meta,:inline))
#             llvmcall(
#                 $((decl, join(instrs, "\n"))),
#                 Vec{$W, $T}, Tuple{Ptr{$T}, UInt32},
#                 ptr, stride % UInt32
#             )
#         end
#     end
#     @generated function vcolumnwisestore!(ptr::Ptr{T}, v::Vec{W,T}, stride::Integer, ::Val{M}, ::Val{N}) where {W,T,M,N}
#         typ = llvmtype(T)
#         W = M * N
#         vtyp = "<$W x $typ>"
#         ptyp = "i$(8sizeof(Int))"
#         # suffix = 'r' * M * 'c' * N * typ
#         suffix = "v$(W)$(T <: Integer ? 'i' : 'f')$(8sizeof(T))"
#         instr = "@llvm.matrix.columnwise.store.$suffix"
#         decl = "declare void $instr($vtyp, $typ*, i32, i32, i32)"
#         instrs = String[]
#         push!(instrs, "%ptr = inttoptr $ptyp %0 to $vtyp*")
#         push!(instrs, "%res = call void $instr($vtyp %1, $vtyp* %ptr, i32 %2, i32 $M, i32 $N)")
#         push!(instrs, "ret void")
#         quote
#             $(Expr(:meta,:inline))
#             llvmcall(
#                 $((decl, join(instrs, "\n"))),
#                 Cvoid, Tuple{Ptr{$T}, Vec{$W,$T}, UInt32},
#                 ptr, v, stride % UInt32
#             )
#         end
#     end
#     @generated function vmatmul(vA::Vec{WA,T}, vB::Vec{WB,T}, ::Val{M}, ::Val{K}, ::Val{N}) where {WA, WB, M, K, N, T}
#         typ = llvmtype(
#         WC = M * N
#         vtypA = "<$WA x $typ>"
#         vtypB = "<$WB x $typ>"
#         vtypC = "<$WC x $typ>"
#         sufftyp = "$(T <: Integer ? 'i' : 'f')$(8sizeof(T))"
#         suffix = "v$(WC)$(sufftyp).v$(WA)$(sufftyp).v$(WB)$(sufftyp)"
#         instr = "@llvm.matrix.multiply.$suffix"
#         decl = "declare $vtypC $instr($vtypA, $vtypB, i32, i32, i32)"
#         instrs = String[]
#         push!(instrs, "%res = call $vtypC $instr($vtypA %0, $vtypB %1, i32 $M, i32 $N, i32 $K)")
#         push!(instrs, "ret $vtypC %res")
#         quote
#             $(Expr(:meta, :inline))
#             llvmcall(
#                 $((decl, join(instrs, "\n"))),
#                 Vec{$WC,$T}, Tuple{Vec{$WA,$T}, Vec{$WB,$T}},
#                 vA, vB
#             )
#         end
#     end

# # else
# end


text:
module Molly

using BioStructures
using CellListMap
using ChainRulesCore
import Chemfiles # Imported to avoid clashing names
using Colors
using Combinatorics
using CUDA
using Distances
using Distributions
using EzXML
using ForwardDiff
using KernelDensity
using NearestNeighbors
import NNlib
using ProgressMeter
using Reexport
using Requires
using Unitful
using Zygote

@reexport using AtomsBase
@reexport using StaticArrays
@reexport using Unitful

using Base.Threads
using LinearAlgebra
using SparseArrays

include("types.jl")
include("cutoffs.jl")
include("spatial.jl")
include("forces.jl")
include("simulators.jl")
include("coupling.jl")
include("neighbors.jl")
include("loggers.jl")
include("analysis.jl")
include("chain_rules.jl")
include("zygote.jl")
include("setup.jl")

function __init__()
    @require GLMakie="e9467ef8-e4e7-5192-8a1a-b1aee30e663a" include("makie.jl")
end

end

ext: jl


text:
include("linear_example.jl")
include("second_order_example.jl")


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestNonlinear

using Test
import MathOptInterface
import ForwardDiff
import LinearAlgebra

const MOI = MathOptInterface
const Nonlinear = MOI.Nonlinear

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

_hessian(f, x) = LinearAlgebra.LowerTriangular(ForwardDiff.hessian(f, x))

function test_copy()
    model = Nonlinear.Model()
    @test_throws(
        ErrorException("Copying nonlinear problems not yet implemented"),
        copy(model),
    )
    return
end

function test_parse_unable()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    input = :(f($x))
    @test_throws(
        ErrorException("Unable to parse: $input"),
        Nonlinear.set_objective(model, input),
    )
    return
end

function test_parse_sin()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    Nonlinear.set_objective(model, :(sin($x)))
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.objective_expr(evaluator) == :(sin(x[$x]))
    return
end

function test_parse_sin_squared()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    Nonlinear.set_objective(model, :(sin($x)^2))
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.objective_expr(evaluator) == :(sin(x[$x])^2)
    return
end

function test_parse_ifelse()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    Nonlinear.set_objective(model, :(ifelse($x, 1, 2)))
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.objective_expr(evaluator) == :(ifelse(x[$x], 1, 2))
    return
end

function test_parse_ifelse_inequality_less()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    Nonlinear.set_objective(model, :(ifelse($x < 1, $x - 1, $x + 1)))
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.objective_expr(evaluator) ==
          :(ifelse(x[$x] < 1, x[$x] - 1, x[$x] + 1))
    return
end

function test_parse_ifelse_inequality_greater()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    Nonlinear.set_objective(model, :(ifelse($x > 1, $x - 1, $x + 1)))
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.objective_expr(evaluator) ==
          :(ifelse(x[$x] > 1, x[$x] - 1, x[$x] + 1))
    return
end

function test_parse_ifelse_comparison()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    Nonlinear.set_objective(model, :(ifelse(0 <= $x <= 1, $x - 1, $x + 1)))
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.objective_expr(evaluator) ==
          :(ifelse(0 <= x[$x] <= 1, x[$x] - 1, x[$x] + 1))
    return
end

function test_parse_ifelse_logic_inequality()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    expr = :(ifelse($x < 1.0 || $x > 2.0, $x - 1.0, $x + 1.0))
    Nonlinear.set_objective(model, expr)
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.objective_expr(evaluator) ==
          :(ifelse(x[$x] < 1.0 || x[$x] > 2.0, x[$x] - 1.0, x[$x] + 1.0))
    return
end

function test_parse_splat_prod()
    model = Nonlinear.Model()
    x = MOI.VariableIndex.(1:3)
    Nonlinear.set_objective(model, :(*($x...)))
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.objective_expr(evaluator) ==
          :(x[$(x[1])] * x[$(x[2])] * x[$(x[3])])
    return
end

function test_parse_splat_top_level()
    model = Nonlinear.Model()
    x = MOI.VariableIndex.(1:3)
    @test_throws(
        ErrorException(
            "Unsupported use of the splatting operator. This is only " *
            "supported in the arguments of a function call.",
        ),
        Nonlinear.set_objective(model, :(x...)),
    )
    return
end

function test_parse_splat_expr()
    model = Nonlinear.Model()
    x = MOI.VariableIndex.(1:3)
    @test_throws(
        ErrorException(
            "Unsupported use of the splatting operator. JuMP supports " *
            "splatting only symbols. For example, `x...` is ok, but " *
            "`(x + 1)...`, `[x; y]...` and `g(f(y)...)` are not.",
        ),
        Nonlinear.set_objective(model, :(*((x ./ 2)...))),
    )
    return
end

function test_parse_univariate_prod()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    Nonlinear.set_objective(model, :(*($x)))
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.objective_expr(evaluator) == :(*(x[$x]))
    return
end

function test_parse_string()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    @test_throws(
        ErrorException(
            "Unexpected object abc of type String in nonlinear expression.",
        ),
        Nonlinear.set_objective(model, :($x + "abc")),
    )
    return
end

function test_parse_array()
    model = Nonlinear.Model()
    x = [MOI.VariableIndex(1)]
    c = [1.0]
    @test_throws(
        ErrorException(
            "Unexpected array $(c') in nonlinear expression. Nonlinear " *
            "expressions may contain only scalar expressions.",
        ),
        Nonlinear.set_objective(model, :($(c') * $x)),
    )
    return
end

function test_parse_unsupported_expression()
    model = Nonlinear.Model()
    x = (y = 1,)
    @test_throws(
        ErrorException("Unsupported expression: $(:(x.y))"),
        Nonlinear.set_objective(model, :(x.y)),
    )
    return
end

function test_moi_variable_parse()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    expr = Nonlinear.parse_expression(model, :($x))
    @test expr.nodes == [Nonlinear.Node(Nonlinear.NODE_MOI_VARIABLE, 1, -1)]
    @test isempty(expr.values)
    return
end

function test_expression_parse()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    ex = Nonlinear.add_expression(model, :(sin($x)^2))
    @test model[ex] isa Nonlinear.Expression
    return
end

function test_parameter_parse()
    model = Nonlinear.Model()
    p = Nonlinear.add_parameter(model, 1.2)
    expr = Nonlinear.parse_expression(model, :($p))
    @test expr.nodes == [Nonlinear.Node(Nonlinear.NODE_PARAMETER, 1, -1)]
    @test isempty(expr.values)
    @test model.parameters == [1.2]
    return
end

function test_parameter_set()
    model = Nonlinear.Model()
    p = Nonlinear.add_parameter(model, 1.2)
    @test model.parameters == [1.2]
    @test model[p] == 1.2
    model[p] = 2.1
    @test model.parameters == [2.1]
    @test model[p] == 2.1
    return
end

function test_set_objective()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    input = :($x^2 + 1)
    Nonlinear.set_objective(model, input)
    @test model.objective == Nonlinear.parse_expression(model, input)
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.objective_expr(evaluator) == :(x[$x]^2 + 1)
    return
end

function test_set_objective_subexpression()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    input = :($x^2 + 1)
    expr = Nonlinear.add_expression(model, input)
    Nonlinear.set_objective(model, :($expr^2))
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.objective_expr(evaluator) == :((x[$x]^2 + 1)^2)
    return
end

function test_set_objective_nested_subexpression()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    input = :($x^2 + 1)
    expr = Nonlinear.add_expression(model, input)
    expr_2 = Nonlinear.add_expression(model, :($expr^2))
    Nonlinear.set_objective(model, :($expr_2^2))
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.objective_expr(evaluator) == :(((x[$x]^2 + 1)^2)^2)
    return
end

function test_set_objective_parameter()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    p = Nonlinear.add_parameter(model, 1.2)
    Nonlinear.set_objective(model, :($x^2 + $p))
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.objective_expr(evaluator) == :(x[$x]^2 + 1.2)
    return
end

function test_add_constraint_less_than()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    func = :($x^2 + 1)
    set = MOI.LessThan(1.0)
    c = Nonlinear.add_constraint(model, func, set)
    @test model[c].set == set
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.constraint_expr(evaluator, 1) == :(x[$x]^2 + 1 <= 1.0)
    return
end

function test_add_constraint_delete()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    c1 = Nonlinear.add_constraint(model, :($x^2 + 1), MOI.LessThan(1.0))
    _ = Nonlinear.add_constraint(model, :(sqrt($x)), MOI.LessThan(1.0))
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.constraint_expr(evaluator, 1) == :(x[$x]^2 + 1 <= 1.0)
    @test MOI.constraint_expr(evaluator, 2) == :(sqrt(x[$x]) <= 1.0)
    Nonlinear.delete(model, c1)
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.constraint_expr(evaluator, 1) == :(sqrt(x[$x]) <= 1.0)
    @test_throws BoundsError MOI.constraint_expr(evaluator, 2)
    return
end

function test_add_constraint_greater_than()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    func = :($x^2 + 1)
    set = MOI.GreaterThan(1.0)
    c = Nonlinear.add_constraint(model, func, set)
    @test model[c].set == set
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.constraint_expr(evaluator, 1) == :(x[$x]^2 + 1 >= 1.0)
    return
end

function test_add_constraint_equal_to()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    func, set = :($x^2 + 1), MOI.EqualTo(1.0)
    c = Nonlinear.add_constraint(model, func, set)
    @test model[c].set == set
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.constraint_expr(evaluator, 1) == :(x[$x]^2 + 1 == 1.0)
    return
end

function test_add_constraint_interval()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    func, set = :($x^2 + 1), MOI.Interval(-1.0, 1.0)
    c = Nonlinear.add_constraint(model, func, set)
    @test model[c].set == set
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, [:ExprGraph])
    @test MOI.constraint_expr(evaluator, 1) == :(-1.0 <= x[$x]^2 + 1 <= 1.0)
    return
end

function test_eval_univariate_function()
    r = Nonlinear.OperatorRegistry()
    @test Nonlinear.eval_univariate_function(r, :+, 1.0) == 1.0
    @test Nonlinear.eval_univariate_function(r, :-, 1.0) == -1.0
    @test Nonlinear.eval_univariate_function(r, :abs, -1.1) == 1.1
    @test Nonlinear.eval_univariate_function(r, :abs, 1.1) == 1.1
    return
end

function test_eval_univariate_gradient()
    r = Nonlinear.OperatorRegistry()
    @test Nonlinear.eval_univariate_gradient(r, :+, 1.2) == 1.0
    @test Nonlinear.eval_univariate_gradient(r, :-, 1.2) == -1.0
    @test Nonlinear.eval_univariate_gradient(r, :abs, -1.1) == -1.0
    @test Nonlinear.eval_univariate_gradient(r, :abs, 1.1) == 1.0
    return
end

function test_eval_univariate_hessian()
    r = Nonlinear.OperatorRegistry()
    @test Nonlinear.eval_univariate_hessian(r, :+, 1.2) == 0.0
    @test Nonlinear.eval_univariate_hessian(r, :-, 1.2) == 0.0
    @test Nonlinear.eval_univariate_hessian(r, :abs, -1.1) == 0.0
    @test Nonlinear.eval_univariate_hessian(r, :abs, 1.1) == 0.0
    return
end

function test_eval_univariate_function_registered_method_error()
    r = Nonlinear.Model()
    f(x::Float64) = sin(x)^2
    @test_throws ErrorException Nonlinear.register_operator(r, :f, 1, f)
    return
end

function test_univariate_function_register_twice()
    r = Nonlinear.Model()
    f(x) = x
    Nonlinear.register_operator(r, :f, 1, f)
    @test_throws(
        ErrorException("Operator f is already registered."),
        Nonlinear.register_operator(r, :f, 1, f),
    )
    @test_throws(
        ErrorException("Operator f is already registered."),
        Nonlinear.register_operator(r, :f, 2, f),
    )
    return
end

function test_multivariate_function_register_twice()
    r = Nonlinear.Model()
    f(x, y) = x + y
    Nonlinear.register_operator(r, :f, 2, f)
    @test_throws(
        ErrorException("Operator f is already registered."),
        Nonlinear.register_operator(r, :f, 1, f),
    )
    @test_throws(
        ErrorException("Operator f is already registered."),
        Nonlinear.register_operator(r, :f, 2, f),
    )
    return
end

function test_auto_register()
    r = Nonlinear.OperatorRegistry()
    f(x, y) = x + y
    @test_throws ErrorException Nonlinear.assert_registered(r, :f, 2)
    @test_logs (:warn,) Nonlinear.register_operator_if_needed(r, :f, 2, f)
    Nonlinear.assert_registered(r, :f, 2)
    return
end

function test_register_univariate_function_return_type()
    r = Nonlinear.OperatorRegistry()
    f(x) = x < 1 ? "x" : x
    @test_throws(
        ErrorException(
            "Expected return type of `Float64` from the user-defined " *
            "function :f, but got `String`.",
        ),
        Nonlinear.register_operator(r, :f, 1, f),
    )
    return
end

function test_eval_univariate_function_return_type()
    r = Nonlinear.OperatorRegistry()
    f(x) = x < 1 ? x : "x"
    Nonlinear.register_operator(r, :f, 1, f)
    @test_throws(
        ErrorException(
            "Expected return type of Float64 from a user-defined function, " *
            "but got String.",
        ),
        Nonlinear.eval_univariate_function(r, :f, 1.2),
    )
    return
end

function test_eval_univariate_function_registered()
    r = Nonlinear.OperatorRegistry()
    f(x) = sin(x)^2
    grad_calls = 0
    f′(x) = (grad_calls += 1; 2 * sin(x) * cos(x))
    hess_calls = 0
    f′′(x) = (hess_calls += 1; 2 * (cos(x)^2 - sin(x)^2))
    Nonlinear.register_operator(r, :f, 1, f)
    x = 1.2
    @test Nonlinear.eval_univariate_function(r, :f, x) ≈ f(x)
    @test Nonlinear.eval_univariate_gradient(r, :f, x) ≈ f′(x)
    @test grad_calls == 1
    @test Nonlinear.eval_univariate_hessian(r, :f, x) ≈ f′′(x)
    @test hess_calls == 1
    return
end

function test_eval_univariate_function_registered_grad()
    r = Nonlinear.OperatorRegistry()
    f(x) = sin(x)^2
    grad_calls = 0
    f′(x) = (grad_calls += 1; 2 * sin(x) * cos(x))
    hess_calls = 0
    f′′(x) = (hess_calls += 1; 2 * (cos(x)^2 - sin(x)^2))
    Nonlinear.register_operator(r, :f, 1, f, f′)
    @test grad_calls == 2
    x = 1.2
    @test Nonlinear.eval_univariate_function(r, :f, x) ≈ f(x)
    @test Nonlinear.eval_univariate_gradient(r, :f, x) ≈ f′(x)
    @test grad_calls == 4
    @test Nonlinear.eval_univariate_hessian(r, :f, x) ≈ f′′(x)
    @test hess_calls == 1
    return
end

function test_eval_univariate_function_registered_grad_hess()
    r = Nonlinear.OperatorRegistry()
    f(x) = sin(x)^2
    grad_calls = 0
    f′(x) = (grad_calls += 1; 2 * sin(x) * cos(x))
    hess_calls = 0
    f′′(x) = (hess_calls += 1; 2 * (cos(x)^2 - sin(x)^2))
    Nonlinear.register_operator(r, :f, 1, f, f′, f′′)
    x = 1.2
    @test Nonlinear.eval_univariate_function(r, :f, x) ≈ f(x)
    @test Nonlinear.eval_univariate_gradient(r, :f, x) ≈ f′(x)
    @test grad_calls == 2
    @test Nonlinear.eval_univariate_hessian(r, :f, x) ≈ f′′(x)
    @test hess_calls == 2
    return
end

function test_eval_multivariate_function()
    r = Nonlinear.OperatorRegistry()
    x = [1.1, 2.2]
    @test Nonlinear.eval_multivariate_function(r, :+, x) ≈ 3.3
    @test Nonlinear.eval_multivariate_function(r, :-, x) ≈ -1.1
    @test Nonlinear.eval_multivariate_function(r, :*, x) ≈ 1.1 * 2.2
    @test Nonlinear.eval_multivariate_function(r, :^, x) ≈ 1.1^2.2
    @test Nonlinear.eval_multivariate_function(r, :/, x) ≈ 1.1 / 2.2
    @test Nonlinear.eval_multivariate_function(r, :ifelse, [1; x]) == 1.1
    @test Nonlinear.eval_multivariate_function(r, :ifelse, [0; x]) == 2.2
    return
end

function test_eval_multivariate_gradient()
    r = Nonlinear.OperatorRegistry()
    x = [1.1, 2.2]
    g = zeros(2)
    Nonlinear.eval_multivariate_gradient(r, :+, g, x)
    @test g == [1.0, 1.0]
    Nonlinear.eval_multivariate_gradient(r, :-, g, x)
    @test g == [1.0, -1.0]
    Nonlinear.eval_multivariate_gradient(r, :*, g, x)
    @test g ≈ [2.2, 1.1]
    Nonlinear.eval_multivariate_gradient(r, :^, g, x)
    @test g ≈ [2.2 * 1.1^1.2, 1.1^2.2 * log(1.1)]
    Nonlinear.eval_multivariate_gradient(r, :^, g, [1.1, 1.0])
    @test g ≈ [1.0, 1.1 * log(1.1)]
    Nonlinear.eval_multivariate_gradient(r, :^, g, [1.1, 2.0])
    @test g ≈ [2.0 * 1.1, 1.1^2.0 * log(1.1)]
    Nonlinear.eval_multivariate_gradient(r, :^, g, [-1.1, 2.0])
    @test g[1] ≈ 2.0 * -1.1
    @test isnan(g[2])
    Nonlinear.eval_multivariate_gradient(r, :/, g, x)
    @test g ≈ [1 / 2.2, -1.1 / 2.2^2]
    g = zeros(3)
    Nonlinear.eval_multivariate_gradient(r, :ifelse, g, [1; x])
    @test g ≈ [0.0, 1.0, 0.0]
    Nonlinear.eval_multivariate_gradient(r, :ifelse, g, [0; x])
    @test g ≈ [0.0, 0.0, 1.0]
    return
end

function test_eval_multivariate_gradient_mult()
    r = Nonlinear.OperatorRegistry()
    x = [1.1, 0.0, 2.2]
    g = zeros(3)
    Nonlinear.eval_multivariate_gradient(r, :*, g, x)
    @test g == [0.0, 1.1 * 2.2, 0.0]
    return
end

function test_eval_multivariate_hessian_prod()
    r = Nonlinear.OperatorRegistry()
    # 1-arg *
    x = [1.1]
    H = LinearAlgebra.LowerTriangular(zeros(1, 1))
    @test Nonlinear.eval_multivariate_hessian(r, :*, H, x)
    @test H ≈ _hessian(x -> *(x[1]), x)
    # 2-arg *
    x = [1.1, 2.2]
    H = LinearAlgebra.LowerTriangular(zeros(2, 2))
    @test (@allocated Nonlinear.eval_multivariate_hessian(r, :*, H, x)) == 0
    @test Nonlinear.eval_multivariate_hessian(r, :*, H, x)
    @test H ≈ _hessian(x -> x[1] * x[2], x)
    # 3-arg *
    x = [1.1, 2.2, 3.3]
    H = LinearAlgebra.LowerTriangular(zeros(3, 3))
    @test Nonlinear.eval_multivariate_hessian(r, :*, H, x)
    @test H ≈ _hessian(x -> x[1] * x[2] * x[3], x)
    return
end

function test_eval_multivariate_hessian_exponentiation()
    r = Nonlinear.OperatorRegistry()
    # ^1.0
    x = [1.1, 1.0]
    H = LinearAlgebra.LowerTriangular(zeros(2, 2))
    @test (@allocated Nonlinear.eval_multivariate_hessian(r, :^, H, x)) == 0
    @test Nonlinear.eval_multivariate_hessian(r, :^, H, x)
    @test H ≈ _hessian(x -> x[1]^x[2], x)
    # ^2.0
    x = [1.1, 2.0]
    H = LinearAlgebra.LowerTriangular(zeros(2, 2))
    @test (@allocated Nonlinear.eval_multivariate_hessian(r, :^, H, x)) == 0
    @test Nonlinear.eval_multivariate_hessian(r, :^, H, x)
    @test H ≈ _hessian(x -> x[1]^x[2], x)
    # 2-arg ^
    x = [1.1, 2.2]
    H = LinearAlgebra.LowerTriangular(zeros(2, 2))
    @test (@allocated Nonlinear.eval_multivariate_hessian(r, :^, H, x)) == 0
    @test Nonlinear.eval_multivariate_hessian(r, :^, H, x)
    @test H ≈ _hessian(x -> x[1]^x[2], x)
    return
end

function test_eval_multivariate_hessian_division()
    r = Nonlinear.OperatorRegistry()
    # 2-arg /
    x = [1.1, 2.2]
    H = LinearAlgebra.LowerTriangular(zeros(2, 2))
    @test (@allocated Nonlinear.eval_multivariate_hessian(r, :/, H, x)) == 0
    @test Nonlinear.eval_multivariate_hessian(r, :/, H, x)
    @test H ≈ _hessian(x -> x[1] / x[2], x)
    return
end

function test_eval_multivariate_function_registered()
    r = Nonlinear.OperatorRegistry()
    f(x...) = x[1]^2 + x[1] * x[2] + x[2]^2
    Nonlinear.register_operator(r, :f, 2, f)
    x = [1.1, 2.2]
    @test Nonlinear.eval_multivariate_function(r, :f, x) ≈ f(x...)
    g = zeros(2)
    Nonlinear.eval_multivariate_gradient(r, :f, g, x)
    @test g ≈ [2 * x[1] + x[2], x[1] + 2 * x[2]]
    H = LinearAlgebra.LowerTriangular(zeros(2, 2))
    @test_throws(
        ErrorException,
        Nonlinear.eval_multivariate_hessian(r, :f, H, x),
    )
    return
end

function test_eval_multivariate_function_method_error()
    r = Nonlinear.OperatorRegistry()
    function f(x...)
        if x[1] > 1
            y = zeros(Float64, 1)
            y[1] = x[1]
            return y[1]^2
        end
        return x[1]^2 + x[1] * x[2] + x[2]^2
    end
    Nonlinear.register_operator(r, :f, 2, f)
    x = [1.1, 2.2]
    g = [0.0, 0.0]
    @test_throws(
        ErrorException,
        Nonlinear.eval_multivariate_gradient(r, :f, g, x),
    )
    return
end

function test_eval_multivariate_function_registered_grad()
    r = Nonlinear.OperatorRegistry()
    f(x...) = x[1]^2 + x[1] * x[2] + x[2]^2
    grad_calls = 0
    function ∇f(g, x...)
        grad_calls += 1
        g[1] = 2 * x[1] + x[2]
        g[2] = x[1] + 2 * x[2]
        return
    end
    Nonlinear.register_operator(r, :f, 2, f, ∇f)
    x = [1.1, 2.2]
    @test Nonlinear.eval_multivariate_function(r, :f, x) ≈ f(x...)
    g = zeros(2)
    Nonlinear.eval_multivariate_gradient(r, :f, g, x)
    @test g ≈ [2 * x[1] + x[2], x[1] + 2 * x[2]]
    @test grad_calls == 1
    H = LinearAlgebra.LowerTriangular(zeros(2, 2))
    @test_throws(
        ErrorException("Hessian is not defined for operator f"),
        Nonlinear.eval_multivariate_hessian(r, :f, H, x),
    )
    return
end

function test_eval_multivariate_function_registered_hessian()
    r = Nonlinear.OperatorRegistry()
    f(x...) = x[1]^2 + x[1] * x[2] + x[2]^2
    grad_calls = 0
    function ∇f(g, x...)
        grad_calls += 1
        g[1] = 2 * x[1] + x[2]
        g[2] = x[1] + 2 * x[2]
        return
    end
    hess_calls = 0
    function ∇²f(H, x...)
        hess_calls += 1
        H[1, 1] = 2.0
        H[2, 1] = 1.0
        H[2, 2] = 2.0
        return
    end
    Nonlinear.register_operator(r, :f, 2, f, ∇f, ∇²f)
    x = [1.1, 2.2]
    @test Nonlinear.eval_multivariate_function(r, :f, x) ≈ f(x...)
    g = zeros(2)
    Nonlinear.eval_multivariate_gradient(r, :f, g, x)
    @test g ≈ [2 * x[1] + x[2], x[1] + 2 * x[2]]
    @test grad_calls == 1
    H = LinearAlgebra.LowerTriangular(zeros(2, 2))
    @test Nonlinear.eval_multivariate_hessian(r, :f, H, x) == true
    @test H == [2.0 0.0; 1.0 2.0]
    @test hess_calls == 1
    return
end

function test_eval_logic_function()
    r = Nonlinear.OperatorRegistry()
    for lhs in (true, false), rhs in (true, false)
        @test Nonlinear.eval_logic_function(r, :&&, lhs, rhs) == (lhs && rhs)
        @test Nonlinear.eval_logic_function(r, :||, lhs, rhs) == (lhs || rhs)
        @test_throws(
            AssertionError,
            Nonlinear.eval_logic_function(r, :⊻, lhs, rhs),
        )
    end
    return
end

function test_eval_comprison_function()
    r = Nonlinear.OperatorRegistry()
    for lhs in (true, false), rhs in (true, false)
        @test Nonlinear.eval_comparison_function(r, :<=, lhs, rhs) ==
              (lhs <= rhs)
        @test Nonlinear.eval_comparison_function(r, :>=, lhs, rhs) ==
              (lhs >= rhs)
        @test Nonlinear.eval_comparison_function(r, :(==), lhs, rhs) ==
              (lhs == rhs)
        @test Nonlinear.eval_comparison_function(r, :<, lhs, rhs) == (lhs < rhs)
        @test Nonlinear.eval_comparison_function(r, :>, lhs, rhs) == (lhs > rhs)
        @test_throws(
            AssertionError,
            Nonlinear.eval_comparison_function(r, :⊻, lhs, rhs),
        )
    end
    return
end

function test_features_available()
    model = Nonlinear.Model()
    evaluator = Nonlinear.Evaluator(model)
    @test MOI.features_available(evaluator) == [:ExprGraph]
    return
end

function test_features_available_Default()
    model = Nonlinear.Model()
    evaluator = Nonlinear.Evaluator(
        model,
        Nonlinear.ExprGraphOnly(),
        MOI.VariableIndex[],
    )
    @test MOI.features_available(evaluator) == [:ExprGraph]
    return
end

function test_add_constraint_ordinal_index()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    constraints = [
        Nonlinear.add_constraint(model, :($x), MOI.LessThan(1.0 * i)) for
        i in 1:4
    ]
    for i in 1:4
        @test MOI.is_valid(model, constraints[i])
    end
    evaluator = Nonlinear.Evaluator(model)
    MOI.initialize(evaluator, Symbol[])
    for i in 1:4
        @test Nonlinear.ordinal_index(evaluator, constraints[i]) == i
    end
    Nonlinear.delete(model, constraints[1])
    Nonlinear.delete(model, constraints[3])
    @test !MOI.is_valid(model, constraints[1])
    @test MOI.is_valid(model, constraints[2])
    @test !MOI.is_valid(model, constraints[3])
    @test MOI.is_valid(model, constraints[4])
    MOI.initialize(evaluator, Symbol[])
    @test Nonlinear.ordinal_index(evaluator, constraints[2]) == 1
    @test_throws(
        ErrorException("Invalid constraint index $(constraints[3])"),
        Nonlinear.ordinal_index(evaluator, constraints[3]),
    )
    @test Nonlinear.ordinal_index(evaluator, constraints[4]) == 2
    return
end

function test_show()
    model = Nonlinear.Model()
    evaluator = Nonlinear.Evaluator(model)
    @test occursin(":ExprGraph", sprint(show, evaluator))
    return
end

function test_evaluate_comparison()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    ex = Nonlinear.add_expression(model, :(ifelse($x < 1, -1.0, 1.0)))
    @test Nonlinear.evaluate(Dict(x => 1.1), model, ex) == 1.0
    @test Nonlinear.evaluate(Dict(x => 0.9), model, ex) == -1.0
    return
end

function test_evaluate_logic()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    ex = Nonlinear.add_expression(model, :(ifelse($x < 0 || $x > 1, -1.0, 1.0)))
    @test Nonlinear.evaluate(Dict(x => 1.1), model, ex) == -1.0
    @test Nonlinear.evaluate(Dict(x => 0.9), model, ex) == 1.0
    @test Nonlinear.evaluate(Dict(x => -0.9), model, ex) == -1.0
    return
end

function test_evaluate_subexpressions()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    p = Nonlinear.add_parameter(model, 1.23)
    ex = Nonlinear.add_expression(model, :(*($p, $x, $x)))
    ex = Nonlinear.add_expression(model, :($ex + sqrt($ex)))
    ex_v = 1.23 * 1.1 * 1.1
    @test Nonlinear.evaluate(Dict(x => 1.1), model, ex) ≈ ex_v + sqrt(ex_v)
    model[p] = 3.21
    ex_v = 3.21 * 1.2 * 1.2
    @test Nonlinear.evaluate(Dict(x => 1.2), model, ex) ≈ ex_v + sqrt(ex_v)
    return
end

function test_NLPBlockData()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    evaluator = Nonlinear.Evaluator(model, Nonlinear.ExprGraphOnly(), [x])
    block = MOI.NLPBlockData(evaluator)
    @test block.has_objective == false
    @test length(block.constraint_bounds) == 0
    return
end

end

TestNonlinear.runtests()

include("ReverseAD.jl")

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
"""
    t[idx...]

Returns a `DNDSparse` containing only the elements of `t` where the given indices (`idx`)
match. If `idx` has the same type as the index tuple of the `t`, then this is
considered a scalar indexing (indexing of a single value). In this case the value
itself is looked up and returned.

"""
function Base.getindex(t::DNDSparse{K}, idxs...) where K
    if typeof(idxs) <: astuple(K)
        _getindex_scalar(t, idxs)
    else
        _getindex(t, idxs)
    end
end

function _getindex_scalar(t::DNDSparse{K,V}, idxs) where {K,V}
    # scalar getindex
    rng = map(x->[x], idxs)
    t = collect(fromchunks(delayedmap(x->x[rng...], t.chunks)))
    if isempty(t)
        throw(KeyError(idxs))
    else
        first(t)
    end
end

function _getindex(t::DNDSparse{K,V}, idxs) where {K,V}
    if length(idxs) != ndims(t)
        error("wrong number of indices")
    end
    for idx in idxs
        isa(idx, AbstractVector) && (issorted(idx) || error("indices must be sorted for ranged/vector indexing"))
    end

    # Subset the chunks
    # this is currently a linear search

    mapchunks(t, keeplengths=false) do chunk
        getindex(chunk, idxs...)
    end |> cache_thunks
end

# update a given domain to include a new key
function update_domain(d::IndexSpace{<:NamedTuple}, key::Tuple)
    knt = namedtuple(keys(first(d))...)(key)
    IndexSpace(Interval(min(first(d.interval), knt),
                        max(last(d.interval), knt)),
               Interval(map(min, first(d.interval), knt),
                        map(max, last(d.interval), knt)),
               Nullable{Int}())
end

function update_domain(d::IndexSpace{<:Tuple}, key::Tuple)
    IndexSpace(Interval(min(first(d.interval), key),
                        max(last(d.interval), key)),
               Interval(map(min, first(d.interval), key),
                        map(max, last(d.interval), key)),
               Nullable{Int}())
end

function insert_row!(x::DNDSparse{K,T}, idxs::Tuple, val) where {K,T}
    perm = sortperm(x.domains, by=last)
    cs = convert(Array{Any}, x.chunks[perm])
    ds = x.domains[perm]
    i = searchsortedfirst(astuple.(last.(x.domains)), idxs)
    if i >= length(cs)
        i = length(cs)
    end

    ds[i] = update_domain(ds[i], idxs)
    tmp = cs[i]
    cs[i] = compute(get_context(), delayed(x->(x[idxs...] = val; x))(cs[i]))
    x.domains = ds
    x.chunks = cs
end

function insert_row!(x::DNDSparse{K,T}, idxs::NamedTuple, val) where {K,T}
    insert_row!(s, astuple(idxs), val)
end

function Base.setindex!(x::DNDSparse, val, idxs...)
    insert_row!(x, idxs, val)
end


text:
using GeometricIntegrators
using ChargedParticleDynamics.GuidingCenter4d.SymmetricLoop
using GeometricExamples

const Δt     = 2.
const ntime  = 25000
const nloop  = 2000
const nsave  = 25
const nplot  = 5
const run_id = "poincare_1st_symmetric_vprk_dt2"

include("../guiding_center_4d_settings_vprk.jl")

tableau_list = get_tableau_list_vprk_projection()

pinv = guiding_center_4d_iode_poincare_invariant_1st(Δt, nloop, ntime, nsave)

include("guiding_center_4d_symmetric_poincare_invariant_1st.jl")


text:
module RLHelpers
    import Distributions.sample
    import Base.append!
    using Statistics: mean, std
    using Flux: Zygote.@adjoint

    include("discount.jl")
    include("ReplayMemory.jl")
    include("PGReplayMemory.jl")
    include("PGEpisode.jl")
    include("UDRLMemory.jl")

    macro loop(n, ex)
        :(for i = 1:$(esc(n))
            $(esc(ex))
        end)
    end

    gather(y, a) = y[CartesianIndex.(a, axes(a, 1))]
    gather(y, a::AbstractMatrix) = gather(y, vec(a))

    @adjoint gather(y, a) = gather(y, a),
        ȳ -> begin
            o = zeros(eltype(ȳ), size(y))
            @inbounds for (i, j) in enumerate(a)
                o[j, i] = ȳ[i]
            end
            return (o, nothing)
        end

    function loglikelihood(x, mu, logstd)
        return -5f-1 .* (((x .- mu) ./ (exp.(logstd) .+ 1f-8)) .^ 2f0 .+ 2f0 .* logstd .+ log(2f0*Float32(π)))
    end

    function softupdate!(target::T, model::T, τ=1f-2) where T
        for f in fieldnames(T)
            softupdate!(getfield(target, f), getfield(model, f), τ)
        end
    end

    function softupdate!(dst::A, src::A, τ=T(1f-2)) where {T, A<:AbstractArray{T}}
        dst .= τ .* src .+ (one(T) - τ) .* dst
    end

    export @loop, gather, softupdate!, loglikelihood
    export ReplayMemory, PGReplayMemory, PGEpisode, UDRLMemory, command,
        finish!, sample, discount!, nstepr!, GAEλ!, adv!, batch
end


text:
export conv, plan_conv, conv_psf, plan_conv_psf
export plan_conv_buffer, plan_conv_psf_buffer

"""
    conv(u, v[, dims])

Convolve `u` with `v` over `dims` dimensions with an FFT based method.
Note, that this method introduces wrap-around artifacts without
proper padding/windowing.

# Arguments
* `u` is an array in real space.
* `v` is the array to be convolved in real space as well.
* Per default `ntuple(+, min(N, M)))` means that we perform the convolution
    over all dimensions of that array which has less dimensions.
    If `dims` is an array with integers, we perform convolution
    only over these dimensions. Eg. `dims=[1,3]` would perform the convolution
    over the first and third dimension. Second dimension is not convolved.

If `u` and `v` are both a real valued array we use `rfft` and hence
the output is real as well.
If either `u` or `v` is complex we use `fft` and output is hence complex.

 # Examples
1D with FFT over all dimensions. We choose `v` to be a delta peak.
Therefore convolution should act as identity.
```jldoctest
julia> u = [1 2 3 4 5]
1×5 Array{Int64,2}:
 1  2  3  4  5
julia> v = [0 0 1 0 0]
1×5 Array{Int64,2}:
 0  0  1  0  0

julia> conv(u, v)
1×5 Matrix{Float64}:
 4.0  5.0  1.0  2.0  3.0
```

2D with FFT with different `dims` arguments.
```jldoctest
julia> u = 1im .* [1 2 3; 4 5 6]
2×3 Matrix{Complex{Int64}}:
 0+1im  0+2im  0+3im
 0+4im  0+5im  0+6im

julia> v = [1im 0 0; 1im 0 0]
2×3 Matrix{Complex{Int64}}:
 0+1im  0+0im  0+0im
 0+1im  0+0im  0+0im

julia> conv(u, v)
2×3 Matrix{ComplexF64}:
 -5.0+0.0im  -7.0+0.0im  -9.0+0.0im
 -5.0+0.0im  -7.0+0.0im  -9.0+0.0im
```
"""
function conv(u::AbstractArray{T, N}, v::AbstractArray{D, M}, dims=ntuple(+, min(N, M))) where {T, D, N, M}
    return ifft(fft(u, dims) .* fft(v, dims), dims)
end

function conv(u::AbstractArray{<:Real, N}, v::AbstractArray{<:Real, M}, dims=ntuple(+, min(N, M))) where {N, M}
    return irfft(rfft(u, dims) .* rfft(v, dims), size(u, dims[1]), dims)
end

"""
    conv_psf(u, psf[, dims])

`conv_psf` is a shorthand for `conv(u,ifftshift(psf))`. For examples see `conv`.
"""
function conv_psf(u::AbstractArray{T, N}, psf::AbstractArray{D, M}, dims=ntuple(+, min(N, M))) where {T, D, N, M}
    return conv(u, ifftshift(psf, dims), dims)
end

 # define custom adjoint for conv
 # so far only defined for the derivative regarding the first component
function ChainRulesCore.rrule(::typeof(conv), u::AbstractArray{T, N}, v::AbstractArray{D, M},
                              dims=ntuple(+, min(N, M))) where {T, D, N, M}
    Y = conv(u, v, dims)
    function conv_pullback(barx)
        z = zero(eltype(u))
        return NoTangent(), conv(barx, conj(v), dims), z, z
    end
    return Y, conv_pullback
end

"""
    plan_conv(u, v [, dims]; kwargs...)

Pre-plan an optimized convolution for arrays shaped like `u` and `v` (based on pre-plan FFT)
along the given dimenions `dims`.
`dims = 1:ndims(u)` per default.
The 0 frequency of `u` must be located at the first entry.

We return two arguments:
The first one is `v_ft` (obtained by `fft(v)` or `rfft(v)`).
The second return is the convolution function `pconv`.
`pconv` itself has two arguments. `pconv(u, v_ft=v_ft)` where `u` is the object and `v_ft` the v_ft.
This function achieves faster convolution than `conv(u, u)`.
Depending whether `u` is real or complex we do `fft`s or `rfft`s
Additionally, it is possible to provide `flags=FFTW.MEASURE` as `kwargs`
to change the planning of the FFT.

# Examples
```jldoctest
julia> u = [1 2 3 4 5]
1×5 Matrix{Int64}:
 1  2  3  4  5

julia> v = [1 0 0 0 0]
1×5 Matrix{Int64}:
 1  0  0  0  0

julia> v_ft, pconv = plan_conv(u, v);

julia> pconv(u, v_ft)
1×5 Matrix{Float64}:
 1.0  2.0  3.0  4.0  5.0

julia> pconv(u)
1×5 Matrix{Float64}:
 1.0  2.0  3.0  4.0  5.0
```
"""
function plan_conv(u::AbstractArray{T1, N}, v::AbstractArray{T2, M}, dims=ntuple(+, N);
                   kwargs...) where {T1, T2, N, M}
    eltype_error(T1, T2)
    plan = get_plan(T1)
    # do the preplanning step
    P = let
        # FFTW.MEASURE flag might overwrite input! Hence copy!
        if (:flags in keys(kwargs) &&
            (getindex(kwargs, :flags) == FFTW.MEASURE || getindex(kwargs, :flags) == FFTW.PATIENT))
            P = plan(copy(u), dims; kwargs...)
        else
            P = plan(u, dims; kwargs...)
        end
    end
    P_inv = inv(P)

    v_ft = fft_or_rfft(T1)(v, dims)
    # construct the efficient conv function
    # P and P_inv can be understood like matrices
    # but their computation is fast
    conv(u, v_ft=v_ft) = p_conv_aux(P, P_inv, u, v_ft)
    return v_ft, conv
end

"""
    plan_conv_buffer(u, v [, dims]; kwargs...)

Similar to [`plan_conv`](@ref) but instead uses buffers to prevent memory allocations.
Not AD friendly!

"""
function plan_conv_buffer(u::AbstractArray{T1, N}, v::AbstractArray{T2, M}, dims=ntuple(+, N);
                   kwargs...) where {T1, T2, N, M}
    eltype_error(T1, T2)
    plan = get_plan(T1)
    # do the preplanning step
    P_u = plan(u, dims; kwargs...)
    P_v = plan(v, dims)

    u_buff = P_u * u
    v_ft = P_v * v
    uv_buff = u_buff .* v_ft

    # for fourier space we need a new plan
    P = plan(u .* v, dims; kwargs...)
    P_inv = inv(P)
    out_buff = P_inv * uv_buff

    # construct the efficient conv function
    # P and P_inv can be understood like matrices
    # but their computation is fast
    function conv(u, v_ft=v_ft)
        mul!(u_buff, P_u, u)
        uv_buff .= u_buff .* v_ft
        mul!(out_buff, P_inv, uv_buff)
        return out_buff
    end

    return v_ft, conv
end

"""
    plan_conv_psf_buffer(u, psf [, dims]; kwargs...) where {T, N}

`plan_conv_psf_buffer` is a shorthand for `plan_conv_buffer(u, ifftshift(psf))`. For examples see `plan_conv`.
"""
function plan_conv_psf_buffer(u::AbstractArray{T, N}, psf::AbstractArray{T, M}, dims=ntuple(+, N);
                       kwargs...) where {T, N, M}
    return plan_conv_buffer(u, ifftshift(psf, dims), dims; kwargs...)
end

"""
    plan_conv_psf(u, psf [, dims]; kwargs...) where {T, N}

`plan_conv_psf` is a shorthand for `plan_conv(u, ifftshift(psf))`. For examples see `plan_conv`.
"""
function plan_conv_psf(u::AbstractArray{T, N}, psf::AbstractArray{T, M}, dims=ntuple(+, N);
                       kwargs...) where {T, N, M}
    return plan_conv(u, ifftshift(psf, dims), dims; kwargs...)
end

function p_conv_aux(P, P_inv, u, v_ft)
    tmp = (P_inv.p * ((P * u) .* v_ft))
    tmp .*= P_inv.scale
    return tmp
end

function ChainRulesCore.rrule(::typeof(p_conv_aux), P, P_inv, u, v)
    Y = p_conv_aux(P, P_inv, u, v)
    function conv_pullback(barx)
        z = zero(eltype(u))
        ∇ = p_conv_aux(P, P_inv, barx, conj(v))
        return NoTangent(), z, z, ∇, z
    end
    return Y, conv_pullback
end

"""
    fft_or_rfft(T)

Small helper function to decide whether a real
or a complex valued FFT is appropriate.
"""
function fft_or_rfft(::Type{<:Real})
    return rfft
end

function fft_or_rfft(::Type{T}) where T
    return fft
end

"""
    get_plan(T)

Small helper function to decide whether a real
or a complex valued FFT plan is appropriate.
"""
function get_plan(::Type{<:Real})
    return plan_rfft
end

function get_plan(::Type{T}) where T
    return plan_fft
end


text:
function count_arrangements(limit::Int)
    k1 = floor(sqrt(limit))
    k2 = floor(limit/4) + 1

    count = k1/2 * (k1/2 - 1)

    for n in (k1 + 1):k2
        count += floor((n - ceil(sqrt(n^2 - limit)))/2)
    end

    println(Int(count))
end

count_arrangements(1000000)

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
"""
Base class for implementating the pool of `T` objects.
The pool allows to reduce the stress on GC by collecting the unneeded objects
(`release!(pool, obj)`) and reusing them later (`acquire!(pool)`).
"""
abstract type AbstractPool{T} end

capacity(pool::AbstractPool) = pool.capacity
Base.eltype(::Type{AbstractPool{T}}) where T = T
Base.eltype(pool::AbstractPool) = eltype(typeof(pool))
Base.length(pool::AbstractPool) = length(pool.objs)

acquire!(pool::AbstractPool) = length(pool) > 0 ? pop!(pool.objs) : newelem(pool)

function release!(pool::AbstractPool{T}, obj::T) where T
    if length(pool) >= capacity(pool)
        #@warn "Pool capacity exceeded" FIXME add parameter to enable it selectively
    else
        push!(pool.objs, obj)
    end
end

"""
The default `AbstarctPool` implementation.
"""
struct Pool{T} <: AbstractPool{T}
    objs::Vector{T}
    capacity::Int
end


text:
using Plots, MAT

nprocs = 4

@views function vizme1D_mpi(nprocs)
    T = []
    for ip = 1:nprocs
        file = matopen("T_$(ip-1).mat"); T_loc = read(file, "T"); close(file)
        nx_i = length(T_loc)-2
        i1   = 1 + (ip-1)*nx_i
        if (ip==1)  T = zeros(nprocs*nx_i)  end
        T[i1:i1+nx_i-1] .= T_loc[2:end-1]
    end
    display(plot(T, legend=false, framestyle=:box, linewidth=3, xlims=(1, length(T)), ylims=(0, 1), xlabel="nx", ylabel="heat", title="diffusion"))
    return
end

vizme1D_mpi(nprocs)

max_stars_repo_name: luraess/geo-hpc-course

max_issues_repo_name: luraess/geo-hpc-course

max_forks_repo_name: luraess/geo-hpc-course


text:
"""
    reconstruction_matrix_meshwidth(f, domain, N)

Assemble inverse discrete filtering matrix from a continuous filter `f` width constant width
\$h(x) = \\Delta x / 2\$.
"""
function reconstruction_matrix_meshwidth(f::TopHatFilter, domain::PeriodicIntervalDomain, N)
    x = discretize(domain, N)
    Δx = 1 // N * (domain.right - domain.left)
    h = f.width
    h₀ = h(x[1])

    all(h.(x) .≈ h(x[1])) || error("Filter width must be constant")

    if Δx .≈ 2h₀
        # Three point stencil
        inds = [-1, 0, 1]
        stencil = [-1 // 24, 13 // 12, -1 // 24]
    elseif Δx .≈ h₀
        # Three point stencil
        inds = [-1, 0, 1]
        stencil = [-1 // 6, 8 // 6, -1 // 6]
    else
        error("Filter width must be equal to mesh width")
    end

    # Five point stencil
    # inds = [-2, -1, 0, 1, 2]
    # stencil = [3//640, -29//480, 1067//960, -29//480, 3//640]

    # Construct banded matrix
    diags = [i => fill(s, N - abs(i)) for (i, s) ∈ zip(inds, stencil)]
    R = spdiagm(diags...)

    # Periodic extension of three point stencil
    R[1, end] = stencil[1]
    R[end, 1] = stencil[end]

    # Periodic extension of five point stencil
    # R[1, [end - 1, end]] = stencil[[1, 2]]
    # R[2, end] = stencil[1]
    # R[end - 1, 1] = stencil[end]
    # R[end, [1, 2]] = stencil[[end-1, end]]

    R
end

function reconstruction_matrix_meshwidth(f::TopHatFilter, domain::ClosedIntervalDomain, N)
    x = discretize(domain, N)
    Δx = 1 // N * (domain.right - domain.left)
    h = f.width
    h₀ = h(x[1])

    all(h.(x) .≈ h(x[1])) || error("Filter width must be constant")
    Δx .≈ 2h₀ || error("Filter width must be equal to mesh width")

    # Three point stencil
    inds = [-1, 0, 1]
    stencil = [-1 // 24, 13 // 12, -1 // 24]

    # Five point stencil
    # inds = [-2, -1, 0, 1, 2]
    # stencil = [3 // 640, -29 // 480, 1067 // 960, -29 // 480, 3 // 640]

    # Construct banded matrix
    diags = [i => fill(s, N + 1 - abs(i)) for (i, s) ∈ zip(inds, stencil)]
    R = spdiagm(diags...)

    # Boundary weights for three point stencil
    R[1, [1, 2]] = [4 // 3, -1 // 3]
    R[end, [end, end - 1]] = [4 // 3, -1 // 3]

    # Boundary weights for five point stencil
    # R[1, [1, 2, 3]] = [23 // 15, -41 // 60, 3 // 20]
    # R[2, [1, 2, 3, 4]] = [-8//105, 319//280, -29//420, 1//168]
    # R[end - 1, [end, end - 1, end - 2, end - 3]] = [-8//105, 319//280, -29//420, 1//168]
    # R[end, [end, end - 1, end - 2]] = [23 // 15, -41 // 60, 3 // 20]

    R
end