text:
module TestFunctions

using MathOptInterface: VectorOfVariables
using Test
using MathOptInterface
import MutableArithmetics

const MOI = MathOptInterface
const MA = MutableArithmetics

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

const w = MOI.VariableIndex(0)
const x = MOI.VariableIndex(1)
const y = MOI.VariableIndex(2)
const z = MOI.VariableIndex(3)

# Number-like but not subtype of `Number`
struct NonNumber
    value::Int
end
Base.:*(a::NonNumber, b::NonNumber) = NonNumber(a.value * b.value)
Base.:+(a::NonNumber, b::NonNumber) = NonNumber(a.value + b.value)
Base.zero(::Type{NonNumber}) = NonNumber(0)
function Base.isapprox(a::NonNumber, b::NonNumber; kws...)
    return isapprox(a.value, b.value; kws...)
end

function test_NonNumber()
    two = NonNumber(2)
    three = NonNumber(3)
    six = NonNumber(6)
    three_x = MOI.Utilities.operate(*, NonNumber, three, x)
    six_x = MOI.Utilities.operate(*, NonNumber, six, x)
    @test six_x ≈ two * three_x
    @test six_x ≈ three_x * two
    return
end

function test_Vectorization_vectorize()
    g = MOI.VectorAffineFunction(
        MOI.VectorAffineTerm.([3, 1], MOI.ScalarAffineTerm.([5, 2], [y, x])),
        [3, 1, 4],
    )
    g1 = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(2, x)], 3)
    g2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Int}[], 1)
    g3 = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(5, y)], 4)
    @test g ≈ MOI.Utilities.vectorize([g1, g2, g3])
    vov = MOI.Utilities.vectorize(MOI.VariableIndex[])
    @test MOI.output_dimension(vov) == 0
    @test vov isa MOI.VectorOfVariables
    aff = MOI.Utilities.vectorize(MOI.ScalarAffineFunction{Int}[])
    @test MOI.output_dimension(aff) == 0
    @test aff isa MOI.VectorAffineFunction{Int}
    quad = MOI.Utilities.vectorize(MOI.ScalarQuadraticFunction{Int}[])
    @test MOI.output_dimension(quad) == 0
    @test quad isa MOI.VectorQuadraticFunction{Int}
    return
end

function test_Vectorization_operate_vcat()
    g = MOI.VectorAffineFunction(
        MOI.VectorAffineTerm.([3, 1], MOI.ScalarAffineTerm.([5, 2], [y, x])),
        [3, 1, 4],
    )
    v = MOI.VectorOfVariables([y, w])
    for T in [Int, Float64]
        @test MOI.VectorOfVariables == MOI.Utilities.promote_operation(
            vcat,
            T,
            typeof(w),
            typeof(v),
            typeof(x),
        )
        vov = MOI.Utilities.operate(vcat, T, w, v, x)
        @test vov.variables == [w, y, w, x]
        @test MOI.VectorOfVariables == MOI.Utilities.promote_operation(
            vcat,
            T,
            typeof(v),
            typeof(w),
            typeof(x),
        )
        vov = MOI.Utilities.operate(vcat, T, v, w, x)
        @test vov.variables == [y, w, w, x]
        @test MOI.VectorOfVariables == MOI.Utilities.promote_operation(
            vcat,
            T,
            typeof(w),
            typeof(x),
            typeof(v),
        )
        vov = MOI.Utilities.operate(vcat, T, w, x, v)
        @test vov.variables == [w, x, y, w]
    end
    f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2, 4], [x, z]), 5)
    g = MOI.VectorAffineFunction(
        MOI.VectorAffineTerm.([3, 1], MOI.ScalarAffineTerm.([5, 2], [y, x])),
        [3, 1, 4],
    )
    @test MOI.Utilities.promote_operation(
        vcat,
        Int,
        typeof(w),
        typeof(f),
        typeof(v),
        Int,
        typeof(g),
        typeof(x),
        Int,
    ) == MOI.VectorAffineFunction{Int}
    F = MOI.Utilities.operate(vcat, Int, w, f, v, 3, g, x, -4)
    expected_terms =
        MOI.VectorAffineTerm.(
            [1, 2, 2, 3, 4, 8, 6, 9],
            MOI.ScalarAffineTerm.(
                [1, 2, 4, 1, 1, 5, 2, 1],
                [w, x, z, y, w, y, x, x],
            ),
        )
    expected_constants = [0, 5, 0, 0, 3, 3, 1, 4, 0, -4]
    F = MOI.Utilities.operate(vcat, Int, w, f, v, 3, g, x, -4)
    @test F.terms == expected_terms
    @test F.constants == expected_constants
    return
end

function test_promote_operation_Quadratic()
    @test MOI.Utilities.promote_operation(
        vcat,
        Int,
        MOI.VectorQuadraticFunction{Int},
        MOI.VariableIndex,
        MOI.ScalarQuadraticFunction{Int},
        MOI.VectorOfVariables,
        Int,
        MOI.ScalarQuadraticFunction{Int},
        MOI.VectorQuadraticFunction{Int},
        MOI.VariableIndex,
        Int,
    ) == MOI.VectorQuadraticFunction{Int}
    return
end

function test_MultirowChange_construction()
    chg1 = MOI.MultirowChange(w, [(Int32(2), 2.0), (Int32(1), 3.0)])
    chg2 = MOI.MultirowChange(w, [(Int64(2), 2.0), (Int64(1), 3.0)])
    @test chg1.variable == chg2.variable
    @test chg1.new_coefficients == chg2.new_coefficients
    return
end

function test_VectorAffineTerm_VectorQuadraticTerm_construction()
    scalaraffine = MOI.ScalarAffineTerm(2.0, z)
    @test MOI.VectorAffineTerm(Int32(3), scalaraffine) ===
          MOI.VectorAffineTerm(Int64(3), scalaraffine)
    @test MOI.VectorAffineTerm{Float64}(Int32(3), scalaraffine) ===
          MOI.VectorAffineTerm(Int64(3), scalaraffine)
    scalarquad = MOI.ScalarQuadraticTerm(2.0, y, z)
    @test MOI.VectorQuadraticTerm(Int32(3), scalarquad) ===
          MOI.VectorQuadraticTerm(Int64(3), scalarquad)
    @test MOI.VectorQuadraticTerm{Float64}(Int32(3), scalarquad) ===
          MOI.VectorQuadraticTerm(Int64(3), scalarquad)
    return
end

function test_eval_variables()
    # We do tests twice to make sure the function is not modified
    vals = Dict(w => 0, x => 3, y => 1, z => 5)
    @test MOI.output_dimension(z) == 1
    @test MOI.Utilities.eval_variables(vi -> vals[vi], z) ≈ 5
    @test MOI.Utilities.eval_variables(vi -> vals[vi], z) ≈ 5
    fvv = MOI.VectorOfVariables([x, z, y])
    @test MOI.output_dimension(fvv) == 3
    @test MOI.Utilities.eval_variables(vi -> vals[vi], fvv) ≈ [3, 5, 1]
    @test MOI.Utilities.eval_variables(vi -> vals[vi], fvv) ≈ [3, 5, 1]
    fsa = MOI.ScalarAffineFunction(
        [
            MOI.ScalarAffineTerm(1.0, x),
            MOI.ScalarAffineTerm(3.0, z),
            MOI.ScalarAffineTerm(2.0, y),
        ],
        2.0,
    )
    @test MOI.output_dimension(fsa) == 1
    @test MOI.Utilities.eval_variables(vi -> vals[vi], fsa) ≈ 22
    @test MOI.Utilities.eval_variables(vi -> vals[vi], fsa) ≈ 22
    fva = MOI.VectorAffineFunction(
        MOI.VectorAffineTerm.(
            [2, 1, 2],
            MOI.ScalarAffineTerm.([1.0, 3.0, 2.0], [x, z, y]),
        ),
        [-3.0, 2.0],
    )
    @test MOI.output_dimension(fva) == 2
    @test MOI.Utilities.eval_variables(vi -> vals[vi], fva) ≈ [12, 7]
    @test MOI.Utilities.eval_variables(vi -> vals[vi], fva) ≈ [12, 7]
    fsq = MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm.(1.0, [x, w, w], [z, z, y]),
        MOI.ScalarAffineTerm.(1.0, [x, y]),
        -3.0,
    )
    @test MOI.output_dimension(fsq) == 1
    @test MOI.Utilities.eval_variables(vi -> vals[vi], fsq) ≈ 16
    @test MOI.Utilities.eval_variables(vi -> vals[vi], fsq) ≈ 16
    fvq = MOI.VectorQuadraticFunction(
        MOI.VectorQuadraticTerm.(
            [1, 2, 2],
            MOI.ScalarQuadraticTerm.(1.0, [x, w, w], [z, z, y]),
        ),
        MOI.VectorAffineTerm.([2, 1], MOI.ScalarAffineTerm.(1.0, [x, y])),
        [-3.0, -2.0],
    )
    @test MOI.output_dimension(fvq) == 2
    @test MOI.Utilities.eval_variables(vi -> vals[vi], fvq) ≈ [13, 1]
    @test MOI.Utilities.eval_variables(vi -> vals[vi], fvq) ≈ [13, 1]
    return
end

function test_substitute_variables()
    # We do tests twice to make sure the function is not modified
    subs = Dict(w => 1.0y + 1.0z, x => 2.0y + 1.0, y => 1.0y, z => -1.0w)
    vals = Dict(w => 0.0, x => 3.0, y => 1.0, z => 5.0)
    subs_vals = Dict(w => 6.0, x => 3.0, y => 1.0, z => 0.0)
    fsa = x + 3.0z + 2.0y + 2.0
    subs_sa = -3.0w + 4.0y + 3.0
    @test MOI.Utilities.eval_variables(vi -> subs_vals[vi], fsa) ==
          MOI.Utilities.eval_variables(vi -> vals[vi], subs_sa)
    @test MOI.Utilities.substitute_variables(vi -> subs[vi], fsa) ≈ subs_sa
    @test MOI.Utilities.substitute_variables(vi -> subs[vi], fsa) ≈ subs_sa
    fva = MOI.Utilities.operate(vcat, Float64, 3.0z - 3.0, x + 2.0y + 2.0)
    subs_va = MOI.Utilities.operate(vcat, Float64, -3.0w - 3.0, 4.0y + 3.0)
    @test MOI.Utilities.eval_variables(vi -> subs_vals[vi], fva) ==
          MOI.Utilities.eval_variables(vi -> vals[vi], subs_va)
    @test MOI.Utilities.substitute_variables(vi -> subs[vi], fva) ≈ subs_va
    @test MOI.Utilities.substitute_variables(vi -> subs[vi], fva) ≈ subs_va
    fsq = 1.0x + 1.0y + 1.0x * z + 1.0w * z + 1.0w * y + 2.0w * w - 3.0
    #     2y + 1 +     y - 2yw - w    + (y + z) * (-w) + (y + z) * y  + 2 (y + z) * (y + z) - 3
    #     2y     +     y - 2yw - w    - yw - zw        + y^2 + zy     + 2y^2 + 2z^2 + 4yz   - 2
    #     3y - w + 3y^2 + 2z^2 - 3yw - zw + 5yz - 2
    subs_sq =
        3.0y - 1.0w + 3.0y * y + 2.0z * z - 3.0y * w - 1.0z * w + 5.0y * z - 2.0
    @test MOI.Utilities.eval_variables(vi -> subs_vals[vi], fsq) ==
          MOI.Utilities.eval_variables(vi -> vals[vi], subs_sq)
    @test MOI.Utilities.substitute_variables(vi -> subs[vi], fsq) ≈ subs_sq
    @test MOI.Utilities.substitute_variables(vi -> subs[vi], fsq) ≈ subs_sq
    fvq = MOI.Utilities.operate(
        vcat,
        Float64,
        1.0y + 1.0x * z - 3.0,
        1.0x + 1.0w * z + 1.0w * y - 2.0,
    )
    subs_vq = MOI.Utilities.operate(
        vcat,
        Float64,
        1.0y - w - 2.0y * w - 3.0,
        2.0y + 1.0y * y - 1.0w * y - 1.0w * z + 1.0y * z - 1.0,
    )
    @test MOI.Utilities.eval_variables(vi -> subs_vals[vi], fvq) ==
          MOI.Utilities.eval_variables(vi -> vals[vi], subs_vq)
    @test MOI.Utilities.substitute_variables(vi -> subs[vi], fvq) ≈ subs_vq
    @test MOI.Utilities.substitute_variables(vi -> subs[vi], fvq) ≈ subs_vq

    complex_aff = 2.0im * y
    # Test that variables can be substituted for `MOI.ScalarAffineFunction{S}`
    # in a `MOI.ScalarAffineFunction{T}` where `S != T`.
    @test MOI.Utilities.substitute_variables(vi -> 1.5x, complex_aff) ≈
          3.0im * x
    float_aff = 2.0 * y
    @test MOI.Utilities.substitute_variables(vi -> 3x, float_aff) ≈ 6.0 * x
    complex_quad = 1.5im * y * y
    @test MOI.Utilities.substitute_variables(vi -> 2x, complex_quad) ≈
          6.0im * x * x
    int_quad = 3 * y * x
    @test MOI.Utilities.substitute_variables(vi -> true * y, int_quad) ≈
          3 * y * y
    return
end

function test_map_indices()
    fsq = MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm.(1.0, [x, w, w], [z, z, y]),
        MOI.ScalarAffineTerm.(1.0, [x, y]),
        -3.0,
    )
    index_map = Dict(x => y, y => z, w => w, z => x)
    gsq = MOI.Utilities.map_indices(index_map, fsq)
    sats = MOI.ScalarAffineTerm.(1.0, [y, z])
    sqts = MOI.ScalarQuadraticTerm.(1.0, [y, w, w], [x, x, z])
    @test gsq.affine_terms == sats
    @test gsq.quadratic_terms == sqts
    @test gsq.constant == -3.0
    fvq = MOI.VectorQuadraticFunction(
        MOI.VectorQuadraticTerm.(
            [1, 2, 2],
            MOI.ScalarQuadraticTerm.(1.0, [x, w, w], [z, z, y]),
        ),
        MOI.VectorAffineTerm.([2, 1], MOI.ScalarAffineTerm.(1.0, [x, y])),
        [-3.0, -2.0],
    )
    gvq = MOI.Utilities.map_indices(index_map, fvq)
    @test gvq.affine_terms == MOI.VectorAffineTerm.([2, 1], sats)
    @test gvq.quadratic_terms == MOI.VectorQuadraticTerm.([1, 2, 2], sqts)
    @test MOI.Utilities.constant_vector(gvq) == [-3.0, -2.0]
    @test MOI.Utilities.map_indices(index_map, :s) == :s
    @test MOI.Utilities.map_indices(index_map, [:a, :b]) == [:a, :b]
    @test MOI.Utilities.map_indices(index_map, "s") == "s"
    @test MOI.Utilities.map_indices(index_map, ["a", "b"]) == ["a", "b"]
    return
end

function test_Conversion_VectorOfVariables_VectorAffineFunction()
    f = MOI.VectorAffineFunction{Int}(MOI.VectorOfVariables([z, x, y]))
    @test f isa MOI.VectorAffineFunction{Int}
    @test f.terms ==
          MOI.VectorAffineTerm.(
        1:3,
        MOI.ScalarAffineTerm.(ones(Int, 3), [z, x, y]),
    )
    @test all(iszero.(MOI.Utilities.constant_vector(f)))
    f = MOI.VectorAffineFunction{Float64}(MOI.VectorOfVariables([x, w]))
    @test f isa MOI.VectorAffineFunction{Float64}
    @test f.terms ==
          MOI.VectorAffineTerm.(1:2, MOI.ScalarAffineTerm.(1.0, [x, w]))
    @test all(iszero.(MOI.Utilities.constant_vector(f)))
    return
end

function test_iteration_and_indexing_on_VectorOfVariables()
    f = MOI.VectorOfVariables([z, w, x, y])
    it = MOI.Utilities.eachscalar(f)
    @test length(it) == 4
    @test eltype(it) == MOI.VariableIndex
    @test collect(it) == [z, w, x, y]
    @test it[2] == w
    @test it[end] == y
    return
end

function test_indexing_on_VectorAffineFunction()
    f = MOI.VectorAffineFunction(
        MOI.VectorAffineTerm.(
            [2, 1, 3, 2, 2, 1, 3, 1, 2],
            MOI.ScalarAffineTerm.(
                [1, 7, 2, 9, 3, 1, 6, 4, 1],
                [x, y, z, z, y, z, x, x, y],
            ),
        ),
        [2, 7, 5],
    )
    it = MOI.Utilities.eachscalar(f)
    @test length(it) == 3
    @test eltype(it) == MOI.ScalarAffineFunction{Int}
    g = it[2]
    @test g isa MOI.ScalarAffineFunction
    @test g.terms == MOI.ScalarAffineTerm.([1, 9, 3, 1], [x, z, y, y])
    @test g.constant == 7
    g = it[1]
    @test g isa MOI.ScalarAffineFunction
    @test g.terms == MOI.ScalarAffineTerm.([7, 1, 4], [y, z, x])
    @test g.constant == 2
    g = it[end]
    @test g isa MOI.ScalarAffineFunction
    @test g.terms == MOI.ScalarAffineTerm.([2, 6], [z, x])
    @test g.constant == 5
    h = it[[3, 1]]
    @test h isa MOI.VectorAffineFunction
    @test sort(h.terms, by = t -> t.output_index) ==
          MOI.VectorAffineTerm.(
        [1, 1, 2, 2, 2],
        MOI.ScalarAffineTerm.([2, 6, 7, 1, 4], [z, x, y, z, x]),
    )
    @test MOI.Utilities.constant_vector(h) == [5, 2]
    F = MOI.Utilities.operate(vcat, Int, it[[1, 2]], it[3])
    @test F isa MOI.VectorAffineFunction{Int}
    @test sort(F.terms, by = t -> t.output_index) ==
          MOI.VectorAffineTerm.(
        [1, 1, 1, 2, 2, 2, 2, 3, 3],
        MOI.ScalarAffineTerm.(
            [7, 1, 4, 1, 9, 3, 1, 2, 6],
            [y, z, x, x, z, y, y, z, x],
        ),
    )
    @test MOI.Utilities.constant_vector(F) == MOI.Utilities.constant_vector(f)
    return
end

function test_indexing_on_VectorQuadraticFunction()
    f = MOI.VectorQuadraticFunction(
        MOI.VectorQuadraticTerm.(
            [2, 3, 1, 2],
            MOI.ScalarQuadraticTerm.([1, 6, 4, 3], [z, x, x, y], [y, z, z, y]),
        ),
        MOI.VectorAffineTerm.(
            [2, 1, 3, 2, 2],
            MOI.ScalarAffineTerm.([1, 7, 2, 9, 3], [x, y, z, z, y]),
        ),
        [2, 7, 5],
    )
    it = MOI.Utilities.eachscalar(f)
    @test length(it) == 3
    @test eltype(it) == MOI.ScalarQuadraticFunction{Int}
    g = it[2]
    @test g isa MOI.ScalarQuadraticFunction
    @test g.affine_terms == MOI.ScalarAffineTerm.([1, 9, 3], [x, z, y])
    @test g.quadratic_terms == MOI.ScalarQuadraticTerm.([1, 3], [z, y], [y, y])
    @test g.constant == 7
    g = it[end]
    @test g isa MOI.ScalarQuadraticFunction
    @test g.affine_terms == MOI.ScalarAffineTerm.([2], [z])
    @test g.quadratic_terms == MOI.ScalarQuadraticTerm.([6], [x], [z])
    @test g.constant == 5
    return
end

function test_Scalar_Variable_isone()
    f = MOI.VariableIndex(0)
    g = MOI.VariableIndex(1)
    @test !isone(f)
    @test !isone(g)
    return
end

function test_Scalar_Variable_iszero()
    f = MOI.VariableIndex(0)
    g = MOI.VariableIndex(1)
    @test !iszero(f)
    @test !iszero(g)
    @test f + 1 ≈ 1 + f
    @test (f + 1.0) - 1.0 ≈ (2.0f) / 2.0
    @test (f - 1.0) + 1.0 ≈ (2.0f) / 2.0
    @test (1.0 + f) - 1.0 ≈ (f * 2.0) / 2.0
    @test 1.0 - (1.0 - f) ≈ (f / 2.0) * 2.0
    return
end

function test_Scalar_Variable_complex()
    f = MOI.VariableIndex(0)
    @test real(f) === f
    @test MA.promote_operation(real, typeof(f)) == typeof(f)
    for T in [Int, Int32, Float64]
        z = MOI.Utilities.operate(imag, T, f)
        @test z isa MOI.ScalarAffineFunction{T}
        @test iszero(z)
        @test MOI.Utilities.promote_operation(imag, T, typeof(f)) == typeof(z)
    end
    @test conj(f) === f
    @test MA.promote_operation(conj, typeof(f)) == typeof(f)
    return
end

function test_Scalar_Affine_zero()
    f = @inferred MOI.Utilities.zero(MOI.ScalarAffineFunction{Float64})
    @test iszero(f)
    @test MOI.Utilities.isapprox_zero(f, 1e-16)
    return
end

function test_Scalar_Affine_complex()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    r = 3x + 4y
    c = 2x + 5y
    f = (1 + 0im) * r + c * im
    @test real(f) ≈ r
    @test MOI.Utilities.operate(imag, Int, f) ≈ c
    @test imag(f) ≈ c
    @test conj(f) ≈ (1 + 0im) * r - c * im
    @test MA.promote_operation(real, typeof(f)) == typeof(r)
    @test MA.promote_operation(imag, typeof(f)) == typeof(c)
    @test MA.promote_operation(conj, typeof(f)) == typeof(f)
    return
end

function test_Scalar_Affine_promote_operation()
    @test MOI.Utilities.promote_operation(-, Int, MOI.VariableIndex) ==
          MOI.ScalarAffineFunction{Int}
    @test MOI.Utilities.promote_operation(
        -,
        Int,
        MOI.ScalarAffineFunction{Int},
    ) == MOI.ScalarAffineFunction{Int}
    @test MOI.Utilities.promote_operation(
        +,
        Float64,
        MOI.VariableIndex,
        MOI.VariableIndex,
    ) == MOI.ScalarAffineFunction{Float64}
    @test MOI.Utilities.promote_operation(
        +,
        Float64,
        MOI.ScalarAffineFunction{Float64},
        Float64,
    ) == MOI.ScalarAffineFunction{Float64}
    @test MOI.Utilities.promote_operation(
        +,
        Int,
        MOI.ScalarAffineFunction{Int},
        MOI.ScalarAffineFunction{Int},
    ) == MOI.ScalarAffineFunction{Int}
    return
end

function test_Scalar_Affine_comparison()
    @test MOI.Utilities.operate(+, Float64, x, z) + 1.0 ≈
          MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.([1, 1e-7, 1], [x, y, z]),
        1.0,
    ) atol = 1e-6
    f1 = MOI.ScalarAffineFunction(
        [MOI.ScalarAffineTerm(1.0, x), MOI.ScalarAffineTerm(1e-7, y)],
        1.0,
    )
    f2 = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 1.0)
    @test f1 ≈ f2 atol = 1e-6
    fdiff = f1 - f2
    @test !iszero(fdiff)
    @test iszero(f1 - f1)
    @test iszero(f2 - f2)
    MOI.Utilities.canonicalize!(fdiff)
    @test !MOI.Utilities.isapprox_zero(fdiff, 1e-8)
    @test MOI.Utilities.isapprox_zero(fdiff, 1e-6)
    return
end

function test_Scalar_Affine_canonical()
    f = MOI.Utilities.canonical(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([2, 1, 3, -2, -3], [y, x, z, x, z]),
            5,
        ),
    )
    @test MOI.output_dimension(f) == 1
    @test f.terms == MOI.ScalarAffineTerm.([-1, 2], [x, y])
    @test f.constant == 5
    f = MOI.Utilities.canonical(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1, 3, 1, 2, -3, 2], [w, y, w, x, x, z]),
            2,
        ) + MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([-1, -2, -2, 3, 2], [y, z, w, x, y]),
            3,
        ),
    )
    @test f.terms == MOI.ScalarAffineTerm.([2, 4], [x, y])
    @test f.constant == 5
end

function test_Scalar_Affine_convert()
    f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 0.5], [x, y]), 0.5)
    @test_throws InexactError convert(MOI.VariableIndex, f)
    @test_throws InexactError MOI.Utilities.convert_approx(MOI.VariableIndex, f)
    @test MOI.Utilities.convert_approx(MOI.VariableIndex, f, tol = 0.5) == x
    @test convert(typeof(f), f) === f
    quad_f = MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm{Float64}[],
        f.terms,
        f.constant,
    )
    @test convert(MOI.ScalarQuadraticFunction{Float64}, f) ≈ quad_f
    for g in [
        convert(MOI.ScalarAffineFunction{Float64}, x),
        convert(MOI.ScalarAffineFunction{Float64}, 1x),
    ]
        @test g isa MOI.ScalarAffineFunction{Float64}
        @test convert(MOI.VariableIndex, g) == x
        @test MOI.Utilities.convert_approx(MOI.VariableIndex, g) == x
    end
    return
end

function test_Scalar_Affine_operate_with_Float_coefficient()
    f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 4.0], [x, y]), 5.0)
    @test f ≈ 2.0f / 2.0
    return
end

function test_Scalar_Affine_operate_with_Int_coefficient()
    f = MOI.Utilities.canonical(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1, 3, 1, 2, -3, 2], [w, y, w, x, x, z]),
            2,
        ) + MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([-1, -2, -2, 3, 2], [y, z, w, x, y]),
            3,
        ),
    )
    @test f === +f
    @test f ≈
          x + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1, 4], [x, y]), 5)
    @test f ≈ f * 1
    @test f ≈
          MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1, 2], [x, y]), 2) *
          2 + 1
    @test f ≈
          x -
          MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-1, -4], [x, y]), -5)
    @test f ≈
          MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([3, 4], [x, y]), 5) - x
    return
end

function test_Scalar_Affine_modification()
    f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2, 4], [x, y]), 0)
    f = MOI.Utilities.modify_function(f, MOI.ScalarConstantChange(6))
    @test f.constant == 6
    g = deepcopy(f)
    @test g ≈ f
    f = MOI.Utilities.modify_function(f, MOI.ScalarCoefficientChange(y, 3))
    @test !(g ≈ f)
    @test g.terms == MOI.ScalarAffineTerm.([2, 4], [x, y])
    @test f.terms == MOI.ScalarAffineTerm.([2, 3], [x, y])
    f = MOI.Utilities.modify_function(f, MOI.ScalarCoefficientChange(x, 0))
    @test f.terms == MOI.ScalarAffineTerm.([3], [y])
    return
end

function test_Scalar_Quadratic_zero()
    f = @inferred MOI.Utilities.zero(MOI.ScalarQuadraticFunction{Float64})
    @test MOI.Utilities.isapprox_zero(f, 1e-16)
    return
end

function test_Scalar_Quadratic_promote_operation()
    @test MOI.Utilities.promote_operation(
        -,
        Int,
        MOI.ScalarQuadraticFunction{Int},
    ) == MOI.ScalarQuadraticFunction{Int}
    @test MOI.Utilities.promote_operation(
        +,
        Int,
        MOI.ScalarQuadraticFunction{Int},
        MOI.ScalarQuadraticFunction{Int},
    ) == MOI.ScalarQuadraticFunction{Int}
    @test MOI.Utilities.promote_operation(
        +,
        Int,
        MOI.ScalarQuadraticFunction{Int},
        MOI.ScalarAffineFunction{Int},
    ) == MOI.ScalarQuadraticFunction{Int}
    @test MOI.Utilities.promote_operation(
        +,
        Int,
        MOI.ScalarAffineFunction{Int},
        MOI.ScalarQuadraticFunction{Int},
    ) == MOI.ScalarQuadraticFunction{Int}
    @test MOI.Utilities.promote_operation(
        *,
        Int,
        MOI.VariableIndex,
        MOI.VariableIndex,
    ) == MOI.ScalarQuadraticFunction{Int}
    @test MOI.Utilities.promote_operation(
        *,
        Float64,
        MOI.VariableIndex,
        MOI.ScalarAffineFunction{Float64},
    ) == MOI.ScalarQuadraticFunction{Float64}
    @test MOI.Utilities.promote_operation(
        *,
        Int,
        MOI.ScalarAffineFunction{Int},
        MOI.VariableIndex,
    ) == MOI.ScalarQuadraticFunction{Int}
    @test MOI.Utilities.promote_operation(
        *,
        Float64,
        MOI.ScalarAffineFunction{Float64},
        MOI.ScalarAffineFunction{Float64},
    ) == MOI.ScalarQuadraticFunction{Float64}
    @test MOI.Utilities.promote_operation(
        /,
        Float64,
        MOI.ScalarQuadraticFunction{Float64},
        Float64,
    ) == MOI.ScalarQuadraticFunction{Float64}
    return
end

function test_Scalar_Quadratic_Comparison_With_iszero()
    f = 7 + 3x + 1x * x + 2y * y + 3x * y
    MOI.Utilities.canonicalize!(f)
    @test !iszero(f)
    @test iszero(0 * f)
    @test iszero(f - f)
    return
end

function test_Scalar_Quadratic_Comparison_With_tolerance()
    f = 7 + 3x + 1x * x + 2y * y + 3x * y
    MOI.Utilities.canonicalize!(f)
    @test !MOI.Utilities.isapprox_zero(f, 1e-8)
    # Test isapprox_zero with zero terms
    @test MOI.Utilities.isapprox_zero(0 * f, 1e-8)
    g = 1.0x * y - (1 + 1e-6) * y * x
    MOI.Utilities.canonicalize!(g)
    @test MOI.Utilities.isapprox_zero(g, 1e-5)
    @test !MOI.Utilities.isapprox_zero(g, 1e-7)
    return
end

function test_Scalar_Quadratic_convert()
    f = 7 + 3x + 1x * x + 2y * y + 3x * y
    MOI.Utilities.canonicalize!(f)
    @test_throws InexactError convert(MOI.VariableIndex, f)
    @test_throws InexactError convert(MOI.ScalarAffineFunction{Int}, f)
    g = convert(MOI.ScalarQuadraticFunction{Float64}, x)
    @test convert(MOI.VariableIndex, g) == x
    return
end

function test_Scalar_Quadratic_Power_Affine()
    f = 7 + 3x + 1x * x + 2y * y + 3x * y
    aff = 1x + 2 + y
    @test isone(@inferred aff^0)
    @test convert(typeof(f), aff) ≈ @inferred aff^1
    @test aff * aff ≈ @inferred aff^2
    err = ArgumentError("Cannot take $(typeof(aff)) to the power 3.")
    @test_throws err aff^3
    return
end

function test_Scalar_Quadratic_Power_Quadratic()
    f = 7 + 3x + 1x * x + 2y * y + 3x * y
    @test isone(@inferred f^0)
    @test f ≈ @inferred f^1
    err = ArgumentError("Cannot take $(typeof(f)) to the power 2.")
    @test_throws err f^2
    return
end

function test_Scalar_Quadratic_no_zero_affine_term()
    # Test that no affine 0y term is created when multiplying 1x by y
    for xy in [1x * y, x * 1y]
        @test isempty(xy.affine_terms)
        @test length(xy.quadratic_terms) == 1
        @test xy.quadratic_terms[1] == MOI.ScalarQuadraticTerm(1, x, y) ||
              xy.quadratic_terms[1] == MOI.ScalarQuadraticTerm(1, y, x)
    end
    for xx in [1x * x, x * 1x]
        @test isempty(xx.affine_terms)
        @test length(xx.quadratic_terms) == 1
        @test xx.quadratic_terms[1] == MOI.ScalarQuadraticTerm(2, x, x)
    end
    return
end

function test_Scalar_Quadratic_operate!()
    q = 1.0x + 1.0y + (1.0x) * z + (1.0w) * z
    @test q ≈ 1.0x + 1.0y + (1.0w) * z + (1.0x) * z
    # This calls
    aff = 1.0x + 1.0y
    # which tries to mutate `aff`, gets a quadratic expression
    # and mutate it with the remaining term
    @test MOI.Utilities.operate!(+, Float64, aff, (1.0x) * z, (1.0w) * z) ≈ q
    return
end

function test_Scalar_Quadratic_operate()
    f = 7 + 3x + 1x * x + 2y * y + 3x * y
    @test f ≈ 7 + (x + 2y) * (1x + y) + 3x
    @test f ≈ -(-7 - 3x) + (x + 2y) * (1x + y)
    @test f ≈ -((x + 2y) * (MOI.Utilities.operate(-, Int, x) - y)) + 3x + 7
    @test f ≈ 7 + MOI.Utilities.operate(*, Int, x, x) + 3x * (y + 1) + 2y * y
    @test f ≈ (x + 2) * (x + 1) + (y + 1) * (2y + 3x) + (5 - 3x - 2y)
    @test f ≈ begin
        MOI.ScalarQuadraticFunction(
            MOI.ScalarQuadraticTerm.([2], [x], [x]),
            [MOI.ScalarAffineTerm(3, x)],
            4,
        ) + MOI.ScalarQuadraticFunction(
            MOI.ScalarQuadraticTerm.([4, 3], [y, x], [y, y]),
            MOI.ScalarAffineTerm{Int}[],
            3,
        )
    end
    @test f ≈ begin
        MOI.ScalarQuadraticFunction(
            MOI.ScalarQuadraticTerm.([2], [x], [x]),
            [MOI.ScalarAffineTerm(3, x)],
            10,
        ) - MOI.ScalarQuadraticFunction(
            MOI.ScalarQuadraticTerm.([-4, -3], [y, x], [y, y]),
            MOI.ScalarAffineTerm{Int}[],
            3,
        )
    end
    @test f ≈ begin
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(3, x)], 5) +
        MOI.ScalarQuadraticFunction(
            MOI.ScalarQuadraticTerm.([2, 4, 3], [x, y, x], [x, y, y]),
            MOI.ScalarAffineTerm{Int}[],
            2,
        )
    end
    @test f ≈ begin
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(3, x)], 5) -
        MOI.ScalarQuadraticFunction(
            MOI.ScalarQuadraticTerm.([-2, -4, -3], [x, y, x], [x, y, y]),
            MOI.ScalarAffineTerm{Int}[],
            -2,
        )
    end
    @test f ≈ begin
        MOI.ScalarQuadraticFunction(
            MOI.ScalarQuadraticTerm.([2, 4, 3], [x, y, x], [x, y, y]),
            MOI.ScalarAffineTerm{Int}[],
            2,
        ) + MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(3, x)], 5)
    end
    @test f ≈ begin
        MOI.ScalarQuadraticFunction(
            MOI.ScalarQuadraticTerm.([2, 4, 3], [x, y, x], [x, y, y]),
            MOI.ScalarAffineTerm{Int}[],
            12,
        ) - MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-3, x)], 5)
    end
    @test f ≈ begin
        MOI.ScalarQuadraticFunction(
            MOI.ScalarQuadraticTerm.([2, 4, 3], [x, y, x], [x, y, y]),
            MOI.ScalarAffineTerm.([2], [x]),
            7,
        ) + x
    end
    @test f ≈
          MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm.([2, 4, 3], [x, y, x], [x, y, y]),
        MOI.ScalarAffineTerm.([3], [x]),
        10,
    ) - 3
    @test f ≈
          2.0 * MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm.([2.0, 4.0, 3.0], [x, y, x], [x, y, y]),
        MOI.ScalarAffineTerm.([3.0], [x]),
        7.0,
    ) / 2.0
    return
end

function test_Scalar_Quadratic_modification()
    f = 7 + 3x + 1x * x + 2y * y + 3x * y
    f = MOI.Utilities.modify_function(f, MOI.ScalarConstantChange(9))
    @test f.constant == 9
    f = MOI.Utilities.modify_function(f, MOI.ScalarCoefficientChange(y, 0))
    @test f.affine_terms == MOI.ScalarAffineTerm.([3], [x])
    g = deepcopy(f)
    @test f ≈ g
    f = MOI.Utilities.modify_function(f, MOI.ScalarCoefficientChange(y, 2))
    @test !(f ≈ g)
    @test g.affine_terms == MOI.ScalarAffineTerm.([3], [x])
    @test f.affine_terms == MOI.ScalarAffineTerm.([3, 2], [x, y])
    return
end

function test_Vector_Variable()
    f = MOI.VectorOfVariables([MOI.VariableIndex(1), MOI.VariableIndex(2)])
    @test real(f) === f
    @test MA.promote_operation(real, typeof(f)) == typeof(f)
    for T in [Int, Int32, Float64]
        z = MOI.Utilities.operate(imag, T, f)
        @test z isa MOI.VectorAffineFunction{T}
        @test isempty(z.terms)
        @test all(iszero, z.constants)
        @test MOI.output_dimension(z) == MOI.output_dimension(f)
        @test MOI.Utilities.promote_operation(imag, T, typeof(f)) == typeof(z)
    end
    @test conj(f) === f
    @test MA.promote_operation(conj, typeof(f)) == typeof(f)
end

function test_Vector_Affine_complex()
    x = MOI.VectorOfVariables([MOI.VariableIndex(1)])
    y = MOI.VectorOfVariables([MOI.VariableIndex(2)])
    r = 3x + 4y
    c = 2x + 5y
    f = (1 + 0im) * r + c * im
    @test real(f) ≈ r
    @test MOI.Utilities.operate(imag, Int, f) ≈ c
    @test imag(f) ≈ c
    @test conj(f) ≈ (1 + 0im) * r - c * im
    @test MA.promote_operation(real, typeof(f)) == typeof(r)
    @test MA.promote_operation(imag, typeof(f)) == typeof(c)
    @test MA.promote_operation(conj, typeof(f)) == typeof(f)
    return
end

function test_Vector_Affine()
    f = MOI.Utilities.canonical(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [2, 1, 2, 1, 1, 2, 2, 2, 2, 1, 1, 2, 1, 2],
                MOI.ScalarAffineTerm.(
                    [3, 2, 3, -3, -1, -2, 3, -2, 1, 3, 5, -2, 0, -1],
                    [x, x, z, y, y, x, y, z, x, y, y, x, x, z],
                ),
            ),
            [5, 7],
        ),
    )
    @test MOI.output_dimension(f) == 2
    @test f.terms ==
          MOI.VectorAffineTerm.(
        [1, 1, 2],
        MOI.ScalarAffineTerm.([2, 4, 3], [x, y, y]),
    )
    @test MOI.Utilities.constant_vector(f) == [5, 7]
    f = MOI.Utilities.modify_function(f, MOI.VectorConstantChange([6, 8]))
    @test MOI.Utilities.constant_vector(f) == [6, 8]
    g = deepcopy(f)
    @test f ≈ g
    f = MOI.Utilities.modify_function(f, MOI.MultirowChange(y, [(2, 9)]))
    @test !(f ≈ g)
    @test f.terms ==
          MOI.VectorAffineTerm.(
        [1, 1, 2],
        MOI.ScalarAffineTerm.([2, 4, 9], [x, y, y]),
    )
    @test g.terms ==
          MOI.VectorAffineTerm.(
        [1, 1, 2],
        MOI.ScalarAffineTerm.([2, 4, 3], [x, y, y]),
    )
    f = MOI.Utilities.modify_function(f, MOI.MultirowChange(y, [(1, 0)]))
    @test f.terms ==
          MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([2, 9], [x, y]))
    return
end

function test_Vector_Quadratic()
    f = MOI.VectorQuadraticFunction(
        MOI.VectorQuadraticTerm.(
            [1, 1, 2],
            MOI.ScalarQuadraticTerm.([1, 2, 3], [x, y, x], [x, y, y]),
        ),
        MOI.VectorAffineTerm.(
            [1, 2, 2],
            MOI.ScalarAffineTerm.([3, 1, 2], [x, x, y]),
        ),
        [7, 3, 4],
    )
    @test MOI.output_dimension(f) == 3
    f = MOI.Utilities.modify_function(f, MOI.VectorConstantChange([10, 11, 12]))
    @test MOI.Utilities.constant_vector(f) == [10, 11, 12]
    f = MOI.Utilities.modify_function(
        f,
        MOI.MultirowChange(y, [(2, 0), (1, 1)]),
    )
    @test f.affine_terms ==
          MOI.VectorAffineTerm.(
        [1, 2, 1],
        MOI.ScalarAffineTerm.([3, 1, 1], [x, x, y]),
    )
    g = deepcopy(f)
    f = MOI.Utilities.modify_function(
        f,
        MOI.MultirowChange(x, [(1, 0), (3, 4)]),
    )
    @test f.affine_terms ==
          MOI.VectorAffineTerm.(
        [2, 1, 3],
        MOI.ScalarAffineTerm.([1, 1, 4], [x, y, x]),
    )
    @test g.affine_terms ==
          MOI.VectorAffineTerm.(
        [1, 2, 1],
        MOI.ScalarAffineTerm.([3, 1, 1], [x, x, y]),
    )
    return
end

function _isapprox_ordered(
    f1::T,
    f2::T,
) where {T<:Union{MOI.ScalarAffineFunction,MOI.VectorAffineFunction}}
    return (
        (MOI.term_indices.(f1.terms) == MOI.term_indices.(f2.terms)) &&
        (MOI.constant(f1) ≈ MOI.constant(f2)) &&
        (MOI.coefficient.(f1.terms) ≈ MOI.coefficient.(f2.terms))
    )
end

function _isapprox_ordered(
    f1::T,
    f2::T,
) where {T<:Union{MOI.ScalarQuadraticFunction,MOI.VectorQuadraticFunction}}
    return (
        (
            MOI.term_indices.(f1.affine_terms) ==
            MOI.term_indices.(f2.affine_terms)
        ) &&
        (
            MOI.term_indices.(f1.quadratic_terms) ==
            MOI.term_indices.(f2.quadratic_terms)
        ) &&
        (MOI.constant(f1) ≈ MOI.constant(f2)) &&
        (
            MOI.coefficient.(f1.affine_terms) ≈
            MOI.coefficient.(f2.affine_terms)
        ) &&
        (
            MOI.coefficient.(f1.quadratic_terms) ≈
            MOI.coefficient.(f2.quadratic_terms)
        )
    )
end

function _test_canonicalization(
    f::T,
    expected::T,
) where {
    T<:Union{
        MOI.ScalarAffineFunction,
        MOI.VectorAffineFunction,
        MOI.ScalarQuadraticFunction,
        MOI.VectorQuadraticFunction,
    },
}
    @test MOI.Utilities.is_canonical(expected)
    g = @inferred(MOI.Utilities.canonical(f))
    @test _isapprox_ordered(g, expected)
    @test MOI.Utilities.is_canonical(g)
    @test MOI.Utilities.is_canonical(expected)
    @test _isapprox_ordered(MOI.Utilities.canonical(g), g)
    @test MOI.Utilities.canonical(g) !== g
    @test @allocated(MOI.Utilities.canonicalize!(f)) == 0
    return
end

function test_canonical_ScalarAffine()
    @test MOI.Utilities.is_canonical(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex.(Int[])),
            1.0,
        ),
    )
    @test !MOI.Utilities.is_canonical(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([0.0], MOI.VariableIndex.([1])),
            1.0,
        ),
    )
    @test !MOI.Utilities.is_canonical(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 2.0], MOI.VariableIndex.([2, 1])),
            2.0,
        ),
    )
    @test !MOI.Utilities.is_canonical(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 0.0], MOI.VariableIndex.([1, 2])),
            2.0,
        ),
    )

    _test_canonicalization(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex[]),
            1.5,
        ),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex[]),
            1.5,
        ),
    )
    _test_canonicalization(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([0.0], [MOI.VariableIndex(1)]),
            -2.0,
        ),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex[]),
            -2.0,
        ),
    )
    _test_canonicalization(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([0.0], [MOI.VariableIndex(2)]),
            -2.0,
        ),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex[]),
            -2.0,
        ),
    )
    _test_canonicalization(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([3.0], [MOI.VariableIndex(2)]),
            1.0,
        ),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([3.0], [MOI.VariableIndex(2)]),
            1.0,
        ),
    )
    _test_canonicalization(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([3.0, 4.0], MOI.VariableIndex.([2, 1])),
            0.0,
        ),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([4.0, 3.0], MOI.VariableIndex.([1, 2])),
            0.0,
        ),
    )
    _test_canonicalization(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 0.1], MOI.VariableIndex.([1, 1])),
            5.0,
        ),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.1], MOI.VariableIndex.([1])),
            5.0,
        ),
    )
    _test_canonicalization(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.(
                [1.0, 0.5, 2.0],
                MOI.VariableIndex.([1, 3, 1]),
            ),
            5.0,
        ),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([3.0, 0.5], MOI.VariableIndex.([1, 3])),
            5.0,
        ),
    )
    _test_canonicalization(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.(
                [1.0, 0.5, 2.0],
                MOI.VariableIndex.([1, 3, 1]),
            ),
            5.0,
        ),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([3.0, 0.5], MOI.VariableIndex.([1, 3])),
            5.0,
        ),
    )
    _test_canonicalization(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.(
                [1.0, 2.0, 3.0, 4.0, 5.0],
                MOI.VariableIndex.([1, 3, 1, 3, 2]),
            ),
            5.0,
        ),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.(
                [4.0, 5.0, 6.0],
                MOI.VariableIndex.([1, 2, 3]),
            ),
            5.0,
        ),
    )
    _test_canonicalization(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, -1.0], MOI.VariableIndex.([7, 7])),
            5.0,
        ),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex.(Int[])),
            5.0,
        ),
    )
    _test_canonicalization(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.(
                [1.0, -1.0, 0.5],
                MOI.VariableIndex.([7, 7, 2]),
            ),
            5.0,
        ),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([0.5], MOI.VariableIndex.([2])),
            5.0,
        ),
    )
    _test_canonicalization(
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.(
                [1.0, -0.5, 0.5, -0.5],
                MOI.VariableIndex.([7, 7, 2, 7]),
            ),
            5.0,
        ),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([0.5], MOI.VariableIndex.([2])),
            5.0,
        ),
    )
    return
end

function test_canonical_VectorAffine()
    @test MOI.Utilities.is_canonical(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                Int[],
                MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex.(Int[])),
            ),
            Float64[],
        ),
    )
    @test !MOI.Utilities.is_canonical(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1],
                MOI.ScalarAffineTerm.([0.0], MOI.VariableIndex.([1])),
            ),
            [1.0],
        ),
    )
    @test !MOI.Utilities.is_canonical(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 2],
                MOI.ScalarAffineTerm.([0.0, 1.0], MOI.VariableIndex.([1, 2])),
            ),
            [1.0],
        ),
    )
    @test !MOI.Utilities.is_canonical(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [2, 1],
                MOI.ScalarAffineTerm.([1.0, 1.0], MOI.VariableIndex.([1, 2])),
            ),
            [1.0],
        ),
    )
    @test !MOI.Utilities.is_canonical(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 1],
                MOI.ScalarAffineTerm.([1.0, -1.0], MOI.VariableIndex.([1, 1])),
            ),
            [1.0],
        ),
    )
    @test !MOI.Utilities.is_canonical(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 1],
                MOI.ScalarAffineTerm.([1.0, 1.0], MOI.VariableIndex.([1, 1])),
            ),
            [1.0],
        ),
    )

    _test_canonicalization(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                Int[],
                MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex.(Int[])),
            ),
            Float64[],
        ),
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                Int[],
                MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex.(Int[])),
            ),
            Float64[],
        ),
    )
    _test_canonicalization(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [3],
                MOI.ScalarAffineTerm.([0.0], MOI.VariableIndex.([5])),
            ),
            [1.0],
        ),
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                Int[],
                MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex.(Int[])),
            ),
            [1.0],
        ),
    )
    _test_canonicalization(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [3, 2],
                MOI.ScalarAffineTerm.([1.0, 2.0], MOI.VariableIndex.([5, 6])),
            ),
            [1.0, 2.0, 3.0],
        ),
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [2, 3],
                MOI.ScalarAffineTerm.([2.0, 1.0], MOI.VariableIndex.([6, 5])),
            ),
            [1.0, 2.0, 3.0],
        ),
    )
    _test_canonicalization(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [3, 3],
                MOI.ScalarAffineTerm.([1.0, 2.0], MOI.VariableIndex.([6, 5])),
            ),
            [1.0, 2.0, 3.0],
        ),
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [3, 3],
                MOI.ScalarAffineTerm.([2.0, 1.0], MOI.VariableIndex.([5, 6])),
            ),
            [1.0, 2.0, 3.0],
        ),
    )
    _test_canonicalization(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [3, 3],
                MOI.ScalarAffineTerm.([1.0, 2.0], MOI.VariableIndex.([1, 1])),
            ),
            [1.0, 2.0, 3.0],
        ),
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [3],
                MOI.ScalarAffineTerm.([3.0], MOI.VariableIndex.([1])),
            ),
            [1.0, 2.0, 3.0],
        ),
    )
    _test_canonicalization(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 1, 1],
                MOI.ScalarAffineTerm.(
                    [1.0, 2.0, 3.0],
                    MOI.VariableIndex.([1, 2, 1]),
                ),
            ),
            [4.0, 5.0, 6.0],
        ),
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 1],
                MOI.ScalarAffineTerm.([4.0, 2.0], MOI.VariableIndex.([1, 2])),
            ),
            [4.0, 5.0, 6.0],
        ),
    )
    _test_canonicalization(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 1, 1],
                MOI.ScalarAffineTerm.(
                    [1.0, 2.0, 3.0],
                    MOI.VariableIndex.([2, 1, 1]),
                ),
            ),
            [4.0, 5.0, 6.0],
        ),
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 1],
                MOI.ScalarAffineTerm.([5.0, 1.0], MOI.VariableIndex.([1, 2])),
            ),
            [4.0, 5.0, 6.0],
        ),
    )
    _test_canonicalization(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [2, 3, 3],
                MOI.ScalarAffineTerm.(
                    [1.0, 2.0, 3.0],
                    MOI.VariableIndex.([1, 1, 1]),
                ),
            ),
            [4.0, 5.0, 6.0],
        ),
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [2, 3],
                MOI.ScalarAffineTerm.([1.0, 5.0], MOI.VariableIndex.([1, 1])),
            ),
            [4.0, 5.0, 6.0],
        ),
    )
    _test_canonicalization(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [2, 3, 3],
                MOI.ScalarAffineTerm.(
                    [1.0, -3.0, 3.0],
                    MOI.VariableIndex.([1, 1, 1]),
                ),
            ),
            [4.0, 5.0, 6.0],
        ),
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [2],
                MOI.ScalarAffineTerm.([1.0], MOI.VariableIndex.([1])),
            ),
            [4.0, 5.0, 6.0],
        ),
    )
    _test_canonicalization(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [2, 3, 3, 3],
                MOI.ScalarAffineTerm.(
                    [1.0, 3.0, -1.0, -2.0],
                    MOI.VariableIndex.([1, 1, 1, 1]),
                ),
            ),
            [4.0, 5.0, 6.0],
        ),
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [2],
                MOI.ScalarAffineTerm.([1.0], MOI.VariableIndex.([1])),
            ),
            [4.0, 5.0, 6.0],
        ),
    )
    return
end

function test_canonical_ScalarQuadratic()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    @test MOI.Utilities.is_canonical(
        convert(MOI.ScalarQuadraticFunction{Float64}, 1.0),
    )
    @test !MOI.Utilities.is_canonical(1.0x * y + 2.0x * x + 2.0)
    @test !MOI.Utilities.is_canonical(1.0x * x + 0.0x * y + 2.0)

    _test_canonicalization(
        convert(MOI.ScalarQuadraticFunction{Float64}, 1.5),
        convert(MOI.ScalarQuadraticFunction{Float64}, 1.5),
    )
    _test_canonicalization(
        0.0x * y - 2.0,
        convert(MOI.ScalarQuadraticFunction{Float64}, -2.0),
    )
    _test_canonicalization(3.0x * y + 4.0x * x + 0.0, 4.0x * x + 3.0x * y + 0.0)
    _test_canonicalization(1.0x * y + 0.1x * y + 5.0, 1.1x * y + 5.0)
    return
end

function test_canonical_VectorQuadratic()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    @test MOI.Utilities.is_canonical(
        MOI.Utilities.operate(
            vcat,
            Float64,
            convert(MOI.ScalarQuadraticFunction{Float64}, 1.0),
        ),
    )
    @test !MOI.Utilities.is_canonical(
        MOI.Utilities.operate(vcat, Float64, 1.0x * y + 2.0x * x + 2.0),
    )
    @test !MOI.Utilities.is_canonical(
        MOI.Utilities.operate(vcat, Float64, 1.0x * x + 0.0x * y + 2.0),
    )

    _test_canonicalization(
        MOI.Utilities.operate(
            vcat,
            Float64,
            1.5,
            0.0x * y - 2.0,
            3.0x * y + 4.0x * x + 0.0,
            1.0x * y + 0.1x * y + 5.0,
        ),
        MOI.Utilities.operate(
            vcat,
            Float64,
            1.5,
            -2.0,
            4.0x * x + 3.0x * y + 0.0,
            1.1x * y + 5.0,
        ),
    )
    return
end

function test_vector_terms()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    at = MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(10.0, x))
    qt = MOI.VectorQuadraticTerm(3, MOI.ScalarQuadraticTerm(6.0, x, y))
    @test MOI.Utilities.operate_term(*, 3.0, at) ==
          MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(30.0, x))
    @test MOI.Utilities.operate_term(*, at, at) ==
          MOI.VectorQuadraticTerm(3, MOI.ScalarQuadraticTerm(200.0, x, x))
    @test MOI.Utilities.operate_term(*, 3.0, qt) ==
          MOI.VectorQuadraticTerm(3, MOI.ScalarQuadraticTerm(18.0, x, y))
    @test MOI.Utilities.operate_term(/, at, 2.0) ==
          MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(5.0, x))
    @test MOI.Utilities.operate_term(/, qt, 3.0) ==
          MOI.VectorQuadraticTerm(3, MOI.ScalarQuadraticTerm(2.0, x, y))
end

function test_vector_promote()
    T = Float64
    for t1 in [MOI.VectorAffineFunction{T}, MOI.VectorOfVariables, Vector{T}]
        for t2 in [MOI.VectorAffineFunction{T}, MOI.VectorOfVariables]
            @test MOI.Utilities.promote_operation(+, T, t1, t2) ==
                  MOI.VectorAffineFunction{T}
            @test MOI.Utilities.promote_operation(-, T, t1, t2) ==
                  MOI.VectorAffineFunction{T}
        end
    end
    for t1 in [
        MOI.VectorQuadraticFunction{T},
        MOI.VectorAffineFunction{T},
        MOI.VectorOfVariables,
        Vector{T},
    ]
        for t2 in [MOI.VectorQuadraticFunction{T}]
            @test MOI.Utilities.promote_operation(+, T, t1, t2) ==
                  MOI.VectorQuadraticFunction{T}
            @test MOI.Utilities.promote_operation(-, T, t1, t2) ==
                  MOI.VectorQuadraticFunction{T}
        end
    end
    for t in [MOI.VectorOfVariables, MOI.VectorAffineFunction{T}]
        @test MOI.Utilities.promote_operation(-, T, t) ==
              MOI.VectorAffineFunction{T}
    end
    t = MOI.VectorQuadraticFunction{T}
    @test MOI.Utilities.promote_operation(-, T, t) ==
          MOI.VectorQuadraticFunction{T}
end

function test_vector_plus_and_neg()
    w = MOI.VariableIndex(0)
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    α = [1, 2, 3]
    v = MOI.VectorOfVariables([y, w, y])
    g = MOI.VectorAffineFunction(
        MOI.VectorAffineTerm.([3, 1], MOI.ScalarAffineTerm.([5, 2], [y, x])),
        [3, 1, 4],
    )
    f = MOI.VectorQuadraticFunction(
        MOI.VectorQuadraticTerm.(
            [1, 1, 2],
            MOI.ScalarQuadraticTerm.([1, 2, 3], [x, y, x], [x, y, y]),
        ),
        MOI.VectorAffineTerm.(
            [1, 2, 2],
            MOI.ScalarAffineTerm.([3, 1, 2], [x, x, y]),
        ),
        [7, 3, 4],
    )
    v_plus_g = MOI.VectorAffineFunction(
        MOI.VectorAffineTerm.(
            [3, 1, 1, 2, 3],
            MOI.ScalarAffineTerm.([5, 2, 1, 1, 1], [y, x, y, w, y]),
        ),
        [3, 1, 4],
    )
    g_plus_α = MOI.VectorAffineFunction(
        MOI.VectorAffineTerm.([3, 1], MOI.ScalarAffineTerm.([5, 2], [y, x])),
        [4, 3, 7],
    )
    α_minus_v = MOI.VectorAffineFunction(
        MOI.VectorAffineTerm.(
            [1, 2, 3],
            MOI.ScalarAffineTerm.([-1, -1, -1], [y, w, y]),
        ),
        [1, 2, 3],
    )
    v_minus_v_plus_v = MOI.VectorAffineFunction(
        MOI.VectorAffineTerm.(
            [1, 2, 3, 1, 2, 3, 1, 2, 3],
            MOI.ScalarAffineTerm.(
                [1, 1, 1, -1, -1, -1, 1, 1, 1],
                [y, w, y, y, w, y, y, w, y],
            ),
        ),
        [0, 0, 0],
    )
    f_plus_α = MOI.VectorQuadraticFunction(
        MOI.VectorQuadraticTerm.(
            [1, 1, 2],
            MOI.ScalarQuadraticTerm.([1, 2, 3], [x, y, x], [x, y, y]),
        ),
        MOI.VectorAffineTerm.(
            [1, 2, 2],
            MOI.ScalarAffineTerm.([3, 1, 2], [x, x, y]),
        ),
        [8, 5, 7],
    )
    f_minus_g = MOI.VectorQuadraticFunction(
        MOI.VectorQuadraticTerm.(
            [1, 1, 2],
            MOI.ScalarQuadraticTerm.([1, 2, 3], [x, y, x], [x, y, y]),
        ),
        MOI.VectorAffineTerm.(
            [1, 2, 2, 3, 1],
            MOI.ScalarAffineTerm.([3, 1, 2, -5, -2], [x, x, y, y, x]),
        ),
        [4, 2, 0],
    )
    @test v + g ≈ v_plus_g
    @test g + α ≈ g_plus_α
    @test α + g ≈ g_plus_α
    @test α - v ≈ α_minus_v
    @test MOI.Utilities.operate(
        +,
        Int,
        MOI.Utilities.operate(-, Int, v, v),
        v,
    ) ≈ v_minus_v_plus_v
    @test f + α ≈ f_plus_α
    @test f - g ≈ f_minus_g
    @test f - f + f - g ≈ f_minus_g
    @test v + f + α - v ≈ f_plus_α
    @test v - f - α - v ≈ -f_plus_α
    @test MOI.Utilities.operate!(-, Int, v, f) - v ≈ -f
    @test (g + v + g + v + f) - (v + g + v + g) ≈ f
    @test v - α ≈ -α_minus_v
    @test g - f ≈ -f_minus_g
    return
end

function test_vector_mult_and_div()
    w = MOI.VariableIndex(0)
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    v = MOI.VectorOfVariables([y, w, y])
    v2 = MOI.Utilities.operate(vcat, Float64, 2.0y, 2.0w, 2.0y)
    f = MOI.Utilities.operate(
        vcat,
        Float64,
        2.0x + 3.0y + 1.0,
        3.0x + 2.0y + 3.0,
    )
    f2 = MOI.Utilities.operate(
        vcat,
        Float64,
        4.0x + 6.0y + 2.0,
        6.0x + 4.0y + 6.0,
    )
    g = MOI.Utilities.operate(
        vcat,
        Float64,
        7.0x * y + 2.0x + 3.0y + 1.0,
        6.0x * x + 5.0y * y + 3.0x + 2.0y + 3.0,
    )
    g2 = MOI.Utilities.operate(
        vcat,
        Float64,
        14.0x * y + 4.0x + 6.0y + 2.0,
        12.0x * x + 10.0y * y + 6.0x + 4.0y + 6.0,
    )
    for (a, a2) in [(v, v2), (f, f2), (g, g2)]
        @test a * 2.0 ≈ a2
        @test 2.0 * a ≈ a2
        @test a / 0.5 ≈ a2
    end
    return
end

function test_modifycoefficient_duplicates()
    model = MOI.Utilities.Model{Float64}()
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction(
        [MOI.ScalarAffineTerm(1.0, x), MOI.ScalarAffineTerm(1.0, x)],
        0.0,
    )
    new_f =
        MOI.Utilities.modify_function(f, MOI.ScalarCoefficientChange(x, 3.0))
    @test new_f ≈ MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(3.0, x)], 0.0)
    return
end

function test_modifycoefficients_duplicates()
    model = MOI.Utilities.Model{Float64}()
    x = MOI.add_variable(model)
    f = MOI.VectorAffineFunction(
        [
            MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x)),
            MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x)),
            MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x)),
            MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(1.0, x)),
        ],
        [0.0, 0.0, 0.0],
    )
    new_f = MOI.Utilities.modify_function(
        f,
        MOI.MultirowChange(x, [(1, 3.0), (2, 0.5)]),
    )
    @test new_f ≈ MOI.VectorAffineFunction(
        [
            MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(3.0, x)),
            MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(0.5, x)),
            MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(1.0, x)),
        ],
        [0.0, 0.0, 0.0],
    )
    return
end

function test_zero_with_output_dimension()
    f = MOI.Utilities.zero_with_output_dimension(Vector{Int}, 2)
    @test f == [0, 0]
    @test f isa Vector{Int}

    g = MOI.Utilities.zero_with_output_dimension(Vector{Float64}, 2)
    @test g == [0.0, 0.0]
    @test g isa Vector{Float64}

    h = MOI.Utilities.zero_with_output_dimension(
        MOI.VectorAffineFunction{Int},
        2,
    )
    @test MOI.output_dimension(h) == 2
    @test isapprox(
        h,
        MOI.VectorAffineFunction{Int}(MOI.VectorAffineTerm{Int}[], [0, 0]),
    )
    @test h isa MOI.VectorAffineFunction{Int}

    i = MOI.Utilities.zero_with_output_dimension(
        MOI.VectorQuadraticFunction{Int},
        2,
    )
    @test MOI.output_dimension(i) == 2
    @test isapprox(
        i,
        MOI.VectorQuadraticFunction{Int}(
            MOI.VectorQuadraticTerm{Int}[],
            MOI.VectorAffineTerm{Int}[],
            [0, 0],
        ),
    )
    @test i isa MOI.VectorQuadraticFunction{Int}
    return
end

"""
    test_SingleVariable_operators()

Test the three Base.:(op) fallbacks for pure SingleVariable operations.
"""
function test_SingleVariable_operators()
    x = MOI.VariableIndex(1)
    @test_throws ErrorException x + x
    @test_throws ErrorException x - x
    @test_throws ErrorException x * x
    return
end

end  # module

TestFunctions.runtests()


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

const IEEEFloat = Union{Float16, Float32, Float64}

## floating point traits ##

"""
    Inf16

Positive infinity of type [`Float16`](@ref).
"""
const Inf16 = bitcast(Float16, 0x7c00)
"""
    NaN16

A not-a-number value of type [`Float16`](@ref).
"""
const NaN16 = bitcast(Float16, 0x7e00)
"""
    Inf32

Positive infinity of type [`Float32`](@ref).
"""
const Inf32 = bitcast(Float32, 0x7f800000)
"""
    NaN32

A not-a-number value of type [`Float32`](@ref).
"""
const NaN32 = bitcast(Float32, 0x7fc00000)
const Inf64 = bitcast(Float64, 0x7ff0000000000000)
const NaN64 = bitcast(Float64, 0x7ff8000000000000)

const Inf = Inf64
"""
    Inf, Inf64

Positive infinity of type [`Float64`](@ref).
"""
Inf, Inf64

const NaN = NaN64
"""
    NaN, NaN64

A not-a-number value of type [`Float64`](@ref).
"""
NaN, NaN64

# bit patterns
reinterpret(::Type{Unsigned}, x::Float64) = reinterpret(UInt64, x)
reinterpret(::Type{Unsigned}, x::Float32) = reinterpret(UInt32, x)
reinterpret(::Type{Unsigned}, x::Float16) = reinterpret(UInt16, x)
reinterpret(::Type{Signed}, x::Float64) = reinterpret(Int64, x)
reinterpret(::Type{Signed}, x::Float32) = reinterpret(Int32, x)
reinterpret(::Type{Signed}, x::Float16) = reinterpret(Int16, x)

sign_mask(::Type{Float64}) =        0x8000_0000_0000_0000
exponent_mask(::Type{Float64}) =    0x7ff0_0000_0000_0000
exponent_one(::Type{Float64}) =     0x3ff0_0000_0000_0000
exponent_half(::Type{Float64}) =    0x3fe0_0000_0000_0000
significand_mask(::Type{Float64}) = 0x000f_ffff_ffff_ffff

sign_mask(::Type{Float32}) =        0x8000_0000
exponent_mask(::Type{Float32}) =    0x7f80_0000
exponent_one(::Type{Float32}) =     0x3f80_0000
exponent_half(::Type{Float32}) =    0x3f00_0000
significand_mask(::Type{Float32}) = 0x007f_ffff

sign_mask(::Type{Float16}) =        0x8000
exponent_mask(::Type{Float16}) =    0x7c00
exponent_one(::Type{Float16}) =     0x3c00
exponent_half(::Type{Float16}) =    0x3800
significand_mask(::Type{Float16}) = 0x03ff

for T in (Float16, Float32, Float64)
    @eval significand_bits(::Type{$T}) = $(trailing_ones(significand_mask(T)))
    @eval exponent_bits(::Type{$T}) = $(sizeof(T)*8 - significand_bits(T) - 1)
    @eval exponent_bias(::Type{$T}) = $(Int(exponent_one(T) >> significand_bits(T)))
    # maximum float exponent
    @eval exponent_max(::Type{$T}) = $(Int(exponent_mask(T) >> significand_bits(T)) - exponent_bias(T))
    # maximum float exponent without bias
    @eval exponent_raw_max(::Type{$T}) = $(Int(exponent_mask(T) >> significand_bits(T)))
end

## conversions to floating-point ##

# TODO: deprecate in 2.0
Float16(x::Integer) = convert(Float16, convert(Float32, x)::Float32)

for t1 in (Float16, Float32, Float64)
    for st in (Int8, Int16, Int32, Int64)
        @eval begin
            (::Type{$t1})(x::($st)) = sitofp($t1, x)
            promote_rule(::Type{$t1}, ::Type{$st}) = $t1
        end
    end
    for ut in (Bool, UInt8, UInt16, UInt32, UInt64)
        @eval begin
            (::Type{$t1})(x::($ut)) = uitofp($t1, x)
            promote_rule(::Type{$t1}, ::Type{$ut}) = $t1
        end
    end
end

Bool(x::Real) = x==0 ? false : x==1 ? true : throw(InexactError(:Bool, Bool, x))

promote_rule(::Type{Float64}, ::Type{UInt128}) = Float64
promote_rule(::Type{Float64}, ::Type{Int128}) = Float64
promote_rule(::Type{Float32}, ::Type{UInt128}) = Float32
promote_rule(::Type{Float32}, ::Type{Int128}) = Float32
promote_rule(::Type{Float16}, ::Type{UInt128}) = Float16
promote_rule(::Type{Float16}, ::Type{Int128}) = Float16

function Float64(x::UInt128)
    x == 0 && return 0.0
    n = 128-leading_zeros(x) # ndigits0z(x,2)
    if n <= 53
        y = ((x % UInt64) << (53-n)) & 0x000f_ffff_ffff_ffff
    else
        y = ((x >> (n-54)) % UInt64) & 0x001f_ffff_ffff_ffff # keep 1 extra bit
        y = (y+1)>>1 # round, ties up (extra leading bit in case of next exponent)
        y &= ~UInt64(trailing_zeros(x) == (n-54)) # fix last bit to round to even
    end
    d = ((n+1022) % UInt64) << 52
    reinterpret(Float64, d + y)
end

function Float64(x::Int128)
    x == 0 && return 0.0
    s = ((x >>> 64) % UInt64) & 0x8000_0000_0000_0000 # sign bit
    x = abs(x) % UInt128
    n = 128-leading_zeros(x) # ndigits0z(x,2)
    if n <= 53
        y = ((x % UInt64) << (53-n)) & 0x000f_ffff_ffff_ffff
    else
        y = ((x >> (n-54)) % UInt64) & 0x001f_ffff_ffff_ffff # keep 1 extra bit
        y = (y+1)>>1 # round, ties up (extra leading bit in case of next exponent)
        y &= ~UInt64(trailing_zeros(x) == (n-54)) # fix last bit to round to even
    end
    d = ((n+1022) % UInt64) << 52
    reinterpret(Float64, s | d + y)
end

function Float32(x::UInt128)
    x == 0 && return 0f0
    n = 128-leading_zeros(x) # ndigits0z(x,2)
    if n <= 24
        y = ((x % UInt32) << (24-n)) & 0x007f_ffff
    else
        y = ((x >> (n-25)) % UInt32) & 0x00ff_ffff # keep 1 extra bit
        y = (y+one(UInt32))>>1 # round, ties up (extra leading bit in case of next exponent)
        y &= ~UInt32(trailing_zeros(x) == (n-25)) # fix last bit to round to even
    end
    d = ((n+126) % UInt32) << 23
    reinterpret(Float32, d + y)
end

function Float32(x::Int128)
    x == 0 && return 0f0
    s = ((x >>> 96) % UInt32) & 0x8000_0000 # sign bit
    x = abs(x) % UInt128
    n = 128-leading_zeros(x) # ndigits0z(x,2)
    if n <= 24
        y = ((x % UInt32) << (24-n)) & 0x007f_ffff
    else
        y = ((x >> (n-25)) % UInt32) & 0x00ff_ffff # keep 1 extra bit
        y = (y+one(UInt32))>>1 # round, ties up (extra leading bit in case of next exponent)
        y &= ~UInt32(trailing_zeros(x) == (n-25)) # fix last bit to round to even
    end
    d = ((n+126) % UInt32) << 23
    reinterpret(Float32, s | d + y)
end

# TODO: optimize
Float16(x::UInt128) = convert(Float16, Float32(x))
Float16(x::Int128)  = convert(Float16, Float32(x))

Float16(x::Float32) = fptrunc(Float16, x)
Float16(x::Float64) = fptrunc(Float16, x)
Float32(x::Float64) = fptrunc(Float32, x)

Float32(x::Float16) = fpext(Float32, x)
Float64(x::Float32) = fpext(Float64, x)
Float64(x::Float16) = fpext(Float64, x)

AbstractFloat(x::Bool)    = Float64(x)
AbstractFloat(x::Int8)    = Float64(x)
AbstractFloat(x::Int16)   = Float64(x)
AbstractFloat(x::Int32)   = Float64(x)
AbstractFloat(x::Int64)   = Float64(x) # LOSSY
AbstractFloat(x::Int128)  = Float64(x) # LOSSY
AbstractFloat(x::UInt8)   = Float64(x)
AbstractFloat(x::UInt16)  = Float64(x)
AbstractFloat(x::UInt32)  = Float64(x)
AbstractFloat(x::UInt64)  = Float64(x) # LOSSY
AbstractFloat(x::UInt128) = Float64(x) # LOSSY

Bool(x::Float16) = x==0 ? false : x==1 ? true : throw(InexactError(:Bool, Bool, x))

"""
    float(x)

Convert a number or array to a floating point data type.
"""
float(x) = AbstractFloat(x)

"""
    float(T::Type)

Return an appropriate type to represent a value of type `T` as a floating point value.
Equivalent to `typeof(float(zero(T)))`.

# Examples
```jldoctest
julia> float(Complex{Int})
ComplexF64 (alias for Complex{Float64})

julia> float(Int)
Float64
```
"""
float(::Type{T}) where {T<:Number} = typeof(float(zero(T)))
float(::Type{T}) where {T<:AbstractFloat} = T

"""
    unsafe_trunc(T, x)

Return the nearest integral value of type `T` whose absolute value is
less than or equal to `x`. If the value is not representable by `T`, an arbitrary value will
be returned.
"""
function unsafe_trunc end

for Ti in (Int8, Int16, Int32, Int64)
    @eval begin
        unsafe_trunc(::Type{$Ti}, x::Float16) = fptosi($Ti, x)
        unsafe_trunc(::Type{$Ti}, x::Float32) = fptosi($Ti, x)
        unsafe_trunc(::Type{$Ti}, x::Float64) = fptosi($Ti, x)
    end
end
for Ti in (UInt8, UInt16, UInt32, UInt64)
    @eval begin
        unsafe_trunc(::Type{$Ti}, x::Float16) = fptoui($Ti, x)
        unsafe_trunc(::Type{$Ti}, x::Float32) = fptoui($Ti, x)
        unsafe_trunc(::Type{$Ti}, x::Float64) = fptoui($Ti, x)
    end
end

function unsafe_trunc(::Type{UInt128}, x::Float64)
    xu = reinterpret(UInt64,x)
    k = Int(xu >> 52) & 0x07ff - 1075
    xu = (xu & 0x000f_ffff_ffff_ffff) | 0x0010_0000_0000_0000
    if k <= 0
        UInt128(xu >> -k)
    else
        UInt128(xu) << k
    end
end
function unsafe_trunc(::Type{Int128}, x::Float64)
    copysign(unsafe_trunc(UInt128,x) % Int128, x)
end

function unsafe_trunc(::Type{UInt128}, x::Float32)
    xu = reinterpret(UInt32,x)
    k = Int(xu >> 23) & 0x00ff - 150
    xu = (xu & 0x007f_ffff) | 0x0080_0000
    if k <= 0
        UInt128(xu >> -k)
    else
        UInt128(xu) << k
    end
end
function unsafe_trunc(::Type{Int128}, x::Float32)
    copysign(unsafe_trunc(UInt128,x) % Int128, x)
end

unsafe_trunc(::Type{UInt128}, x::Float16) = unsafe_trunc(UInt128, Float32(x))
unsafe_trunc(::Type{Int128}, x::Float16) = unsafe_trunc(Int128, Float32(x))

# matches convert methods
# also determines floor, ceil, round
trunc(::Type{Signed}, x::Float16) = trunc(Int,x)
trunc(::Type{Signed}, x::Float32) = trunc(Int,x)
trunc(::Type{Signed}, x::Float64) = trunc(Int,x)
trunc(::Type{Unsigned}, x::Float16) = trunc(UInt,x)
trunc(::Type{Unsigned}, x::Float32) = trunc(UInt,x)
trunc(::Type{Unsigned}, x::Float64) = trunc(UInt,x)
trunc(::Type{Integer}, x::Float16) = trunc(Int,x)
trunc(::Type{Integer}, x::Float32) = trunc(Int,x)
trunc(::Type{Integer}, x::Float64) = trunc(Int,x)

# fallbacks
floor(::Type{T}, x::AbstractFloat) where {T<:Integer} = trunc(T,round(x, RoundDown))
ceil(::Type{T}, x::AbstractFloat) where {T<:Integer} = trunc(T,round(x, RoundUp))
round(::Type{T}, x::AbstractFloat) where {T<:Integer} = trunc(T,round(x, RoundNearest))

round(x::Float64, r::RoundingMode{:ToZero})  = trunc_llvm(x)
round(x::Float32, r::RoundingMode{:ToZero})  = trunc_llvm(x)
round(x::Float16, r::RoundingMode{:ToZero})  = trunc_llvm(x)
round(x::Float64, r::RoundingMode{:Down})    = floor_llvm(x)
round(x::Float32, r::RoundingMode{:Down})    = floor_llvm(x)
round(x::Float16, r::RoundingMode{:Down})    = floor_llvm(x)
round(x::Float64, r::RoundingMode{:Up})      = ceil_llvm(x)
round(x::Float32, r::RoundingMode{:Up})      = ceil_llvm(x)
round(x::Float16, r::RoundingMode{:Up})      = ceil_llvm(x)
round(x::Float64, r::RoundingMode{:Nearest}) = rint_llvm(x)
round(x::Float32, r::RoundingMode{:Nearest}) = rint_llvm(x)
round(x::Float16, r::RoundingMode{:Nearest}) = rint_llvm(x)

## floating point promotions ##
promote_rule(::Type{Float32}, ::Type{Float16}) = Float32
promote_rule(::Type{Float64}, ::Type{Float16}) = Float64
promote_rule(::Type{Float64}, ::Type{Float32}) = Float64

widen(::Type{Float16}) = Float32
widen(::Type{Float32}) = Float64

## floating point arithmetic ##
-(x::Float64) = neg_float(x)
-(x::Float32) = neg_float(x)
-(x::Float16) = neg_float(x)

+(x::Float16, y::Float16) = add_float(x, y)
+(x::Float32, y::Float32) = add_float(x, y)
+(x::Float64, y::Float64) = add_float(x, y)
-(x::Float16, y::Float16) = sub_float(x, y)
-(x::Float32, y::Float32) = sub_float(x, y)
-(x::Float64, y::Float64) = sub_float(x, y)
*(x::Float16, y::Float16) = mul_float(x, y)
*(x::Float32, y::Float32) = mul_float(x, y)
*(x::Float64, y::Float64) = mul_float(x, y)
/(x::Float16, y::Float16) = div_float(x, y)
/(x::Float32, y::Float32) = div_float(x, y)
/(x::Float64, y::Float64) = div_float(x, y)

muladd(x::Float16, y::Float16, z::Float16) = muladd_float(x, y, z)
muladd(x::Float32, y::Float32, z::Float32) = muladd_float(x, y, z)
muladd(x::Float64, y::Float64, z::Float64) = muladd_float(x, y, z)

# TODO: faster floating point div?
# TODO: faster floating point fld?
# TODO: faster floating point mod?

rem(x::Float16, y::Float16) = rem_float(x, y)
rem(x::Float32, y::Float32) = rem_float(x, y)
rem(x::Float64, y::Float64) = rem_float(x, y)

cld(x::T, y::T) where {T<:AbstractFloat} = -fld(-x,y)

function mod(x::T, y::T) where T<:AbstractFloat
    r = rem(x,y)
    if r == 0
        copysign(r,y)
    elseif (r > 0) ⊻ (y > 0)
        r+y
    else
        r
    end
end

## floating point comparisons ##
==(x::Float16, y::Float16) = eq_float(x, y)
==(x::Float32, y::Float32) = eq_float(x, y)
==(x::Float64, y::Float64) = eq_float(x, y)
!=(x::Float16, y::Float16) = ne_float(x, y)
!=(x::Float32, y::Float32) = ne_float(x, y)
!=(x::Float64, y::Float64) = ne_float(x, y)
<( x::Float16, y::Float16) = lt_float(x, y)
<( x::Float32, y::Float32) = lt_float(x, y)
<( x::Float64, y::Float64) = lt_float(x, y)
<=(x::Float16, y::Float16) = le_float(x, y)
<=(x::Float32, y::Float32) = le_float(x, y)
<=(x::Float64, y::Float64) = le_float(x, y)

isequal(x::Float16, y::Float16) = fpiseq(x, y)
isequal(x::Float32, y::Float32) = fpiseq(x, y)
isequal(x::Float64, y::Float64) = fpiseq(x, y)
isless( x::Float16, y::Float16) = fpislt(x, y)
isless( x::Float32, y::Float32) = fpislt(x, y)
isless( x::Float64, y::Float64) = fpislt(x, y)

# Exact Float (Tf) vs Integer (Ti) comparisons
# Assumes:
# - typemax(Ti) == 2^n-1
# - typemax(Ti) can't be exactly represented by Tf:
#   => Tf(typemax(Ti)) == 2^n or Inf
# - typemin(Ti) can be exactly represented by Tf
#
# 1. convert y::Ti to float fy::Tf
# 2. perform Tf comparison x vs fy
# 3. if x == fy, check if (1) resulted in rounding:
#  a. convert fy back to Ti and compare with original y
#  b. unsafe_convert undefined behaviour if fy == Tf(typemax(Ti))
#     (but consequently x == fy > y)
for Ti in (Int64,UInt64,Int128,UInt128)
    for Tf in (Float16,Float32,Float64)
        @eval begin
            function ==(x::$Tf, y::$Ti)
                fy = ($Tf)(y)
                (x == fy) & (fy != $(Tf(typemax(Ti)))) & (y == unsafe_trunc($Ti,fy))
            end
            ==(y::$Ti, x::$Tf) = x==y

            function <(x::$Ti, y::$Tf)
                fx = ($Tf)(x)
                (fx < y) | ((fx == y) & ((fx == $(Tf(typemax(Ti)))) | (x < unsafe_trunc($Ti,fx)) ))
            end
            function <=(x::$Ti, y::$Tf)
                fx = ($Tf)(x)
                (fx < y) | ((fx == y) & ((fx == $(Tf(typemax(Ti)))) | (x <= unsafe_trunc($Ti,fx)) ))
            end

            function <(x::$Tf, y::$Ti)
                fy = ($Tf)(y)
                (x < fy) | ((x == fy) & (fy < $(Tf(typemax(Ti)))) & (unsafe_trunc($Ti,fy) < y))
            end
            function <=(x::$Tf, y::$Ti)
                fy = ($Tf)(y)
                (x < fy) | ((x == fy) & (fy < $(Tf(typemax(Ti)))) & (unsafe_trunc($Ti,fy) <= y))
            end
        end
    end
end
for op in (:(==), :<, :<=)
    @eval begin
        ($op)(x::Float16, y::Union{Int128,UInt128,Int64,UInt64}) = ($op)(Float64(x), Float64(y))
        ($op)(x::Union{Int128,UInt128,Int64,UInt64}, y::Float16) = ($op)(Float64(x), Float64(y))

        ($op)(x::Union{Float16,Float32}, y::Union{Int32,UInt32}) = ($op)(Float64(x), Float64(y))
        ($op)(x::Union{Int32,UInt32}, y::Union{Float16,Float32}) = ($op)(Float64(x), Float64(y))

        ($op)(x::Float16, y::Union{Int16,UInt16}) = ($op)(Float32(x), Float32(y))
        ($op)(x::Union{Int16,UInt16}, y::Float16) = ($op)(Float32(x), Float32(y))
    end
end

abs(x::Float16) = abs_float(x)
abs(x::Float32) = abs_float(x)
abs(x::Float64) = abs_float(x)

"""
    isnan(f) -> Bool

Test whether a number value is a NaN, an indeterminate value which is neither an infinity
nor a finite number ("not a number").
"""
isnan(x::AbstractFloat) = (x != x)::Bool
isnan(x::Number) = false

isfinite(x::AbstractFloat) = x - x == 0
isfinite(x::Real) = decompose(x)[3] != 0
isfinite(x::Integer) = true

"""
    isinf(f) -> Bool

Test whether a number is infinite.
"""
isinf(x::Real) = !isnan(x) & !isfinite(x)

const hx_NaN = hash_uint64(reinterpret(UInt64, NaN))
let Tf = Float64, Tu = UInt64, Ti = Int64
    @eval function hash(x::$Tf, h::UInt)
        # see comments on trunc and hash(Real, UInt)
        if $(Tf(typemin(Ti))) <= x < $(Tf(typemax(Ti)))
            xi = fptosi($Ti, x)
            if isequal(xi, x)
                return hash(xi, h)
            end
        elseif $(Tf(typemin(Tu))) <= x < $(Tf(typemax(Tu)))
            xu = fptoui($Tu, x)
            if isequal(xu, x)
                return hash(xu, h)
            end
        elseif isnan(x)
            return hx_NaN ⊻ h # NaN does not have a stable bit pattern
        end
        return hash_uint64(bitcast(UInt64, x)) - 3h
    end
end

hash(x::Float32, h::UInt) = hash(Float64(x), h)
hash(x::Float16, h::UInt) = hash(Float64(x), h)

## generic hashing for rational values ##

function hash(x::Real, h::UInt)
    # decompose x as num*2^pow/den
    num, pow, den = decompose(x)

    # handle special values
    num == 0 && den == 0 && return hash(NaN, h)
    num == 0 && return hash(ifelse(den > 0, 0.0, -0.0), h)
    den == 0 && return hash(ifelse(num > 0, Inf, -Inf), h)

    # normalize decomposition
    if den < 0
        num = -num
        den = -den
    end
    z = trailing_zeros(num)
    if z != 0
        num >>= z
        pow += z
    end
    z = trailing_zeros(den)
    if z != 0
        den >>= z
        pow -= z
    end

    # handle values representable as Int64, UInt64, Float64
    if den == 1
        left = ndigits0z(num,2) + pow
        right = trailing_zeros(num) + pow
        if -1074 <= right
            if 0 <= right && left <= 64
                left <= 63                     && return hash(Int64(num) << Int(pow), h)
                signbit(num) == signbit(den)   && return hash(UInt64(num) << Int(pow), h)
            end # typemin(Int64) handled by Float64 case
            left <= 1024 && left - right <= 53 && return hash(ldexp(Float64(num),pow), h)
        end
    end

    # handle generic rational values
    h = hash_integer(den, h)
    h = hash_integer(pow, h)
    h = hash_integer(num, h)
    return h
end

#=
`decompose(x)`: non-canonical decomposition of rational values as `num*2^pow/den`.

The decompose function is the point where rational-valued numeric types that support
hashing hook into the hashing protocol. `decompose(x)` should return three integer
values `num, pow, den`, such that the value of `x` is mathematically equal to

    num*2^pow/den

The decomposition need not be canonical in the sense that it just needs to be *some*
way to express `x` in this form, not any particular way – with the restriction that
`num` and `den` may not share any odd common factors. They may, however, have powers
of two in common – the generic hashing code will normalize those as necessary.

Special values:

 - `x` is zero: `num` should be zero and `den` should have the same sign as `x`
 - `x` is infinite: `den` should be zero and `num` should have the same sign as `x`
 - `x` is not a number: `num` and `den` should both be zero
=#

decompose(x::Integer) = x, 0, 1

function decompose(x::Float16)::NTuple{3,Int}
    isnan(x) && return 0, 0, 0
    isinf(x) && return ifelse(x < 0, -1, 1), 0, 0
    n = reinterpret(UInt16, x)
    s = (n & 0x03ff) % Int16
    e = ((n & 0x7c00) >> 10) % Int
    s |= Int16(e != 0) << 10
    d = ifelse(signbit(x), -1, 1)
    s, e - 25 + (e == 0), d
end

function decompose(x::Float32)::NTuple{3,Int}
    isnan(x) && return 0, 0, 0
    isinf(x) && return ifelse(x < 0, -1, 1), 0, 0
    n = reinterpret(UInt32, x)
    s = (n & 0x007fffff) % Int32
    e = ((n & 0x7f800000) >> 23) % Int
    s |= Int32(e != 0) << 23
    d = ifelse(signbit(x), -1, 1)
    s, e - 150 + (e == 0), d
end

function decompose(x::Float64)::Tuple{Int64, Int, Int}
    isnan(x) && return 0, 0, 0
    isinf(x) && return ifelse(x < 0, -1, 1), 0, 0
    n = reinterpret(UInt64, x)
    s = (n & 0x000fffffffffffff) % Int64
    e = ((n & 0x7ff0000000000000) >> 52) % Int
    s |= Int64(e != 0) << 52
    d = ifelse(signbit(x), -1, 1)
    s, e - 1075 + (e == 0), d
end

"""
    precision(num::AbstractFloat)

Get the precision of a floating point number, as defined by the effective number of bits in
the significand.
"""
function precision end

precision(::Type{Float16}) = 11
precision(::Type{Float32}) = 24
precision(::Type{Float64}) = 53
precision(::T) where {T<:AbstractFloat} = precision(T)

"""
    uabs(x::Integer)

Return the absolute value of `x`, possibly returning a different type should the
operation be susceptible to overflow. This typically arises when `x` is a two's complement
signed integer, so that `abs(typemin(x)) == typemin(x) < 0`, in which case the result of
`uabs(x)` will be an unsigned integer of the same size.
"""
uabs(x::Integer) = abs(x)
uabs(x::BitSigned) = unsigned(abs(x))

"""
    nextfloat(x::AbstractFloat, n::Integer)

The result of `n` iterative applications of `nextfloat` to `x` if `n >= 0`, or `-n`
applications of `prevfloat` if `n < 0`.
"""
function nextfloat(f::IEEEFloat, d::Integer)
    F = typeof(f)
    fumax = reinterpret(Unsigned, F(Inf))
    U = typeof(fumax)

    isnan(f) && return f
    fi = reinterpret(Signed, f)
    fneg = fi < 0
    fu = unsigned(fi & typemax(fi))

    dneg = d < 0
    da = uabs(d)
    if da > typemax(U)
        fneg = dneg
        fu = fumax
    else
        du = da % U
        if fneg ⊻ dneg
            if du > fu
                fu = min(fumax, du - fu)
                fneg = !fneg
            else
                fu = fu - du
            end
        else
            if fumax - fu < du
                fu = fumax
            else
                fu = fu + du
            end
        end
    end
    if fneg
        fu |= sign_mask(F)
    end
    reinterpret(F, fu)
end

"""
    nextfloat(x::AbstractFloat)

Return the smallest floating point number `y` of the same type as `x` such `x < y`. If no
such `y` exists (e.g. if `x` is `Inf` or `NaN`), then return `x`.
"""
nextfloat(x::AbstractFloat) = nextfloat(x,1)

"""
    prevfloat(x::AbstractFloat, n::Integer)

The result of `n` iterative applications of `prevfloat` to `x` if `n >= 0`, or `-n`
applications of `nextfloat` if `n < 0`.
"""
prevfloat(x::AbstractFloat, d::Integer) = nextfloat(x, -d)

"""
    prevfloat(x::AbstractFloat)

Return the largest floating point number `y` of the same type as `x` such `y < x`. If no
such `y` exists (e.g. if `x` is `-Inf` or `NaN`), then return `x`.
"""
prevfloat(x::AbstractFloat) = nextfloat(x,-1)

for Ti in (Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UInt128)
    for Tf in (Float16, Float32, Float64)
        if Ti <: Unsigned || sizeof(Ti) < sizeof(Tf)
            # Here `Tf(typemin(Ti))-1` is exact, so we can compare the lower-bound
            # directly. `Tf(typemax(Ti))+1` is either always exactly representable, or
            # rounded to `Inf` (e.g. when `Ti==UInt128 && Tf==Float32`).
            @eval begin
                function trunc(::Type{$Ti},x::$Tf)
                    if $(Tf(typemin(Ti))-one(Tf)) < x < $(Tf(typemax(Ti))+one(Tf))
                        return unsafe_trunc($Ti,x)
                    else
                        throw(InexactError(:trunc, $Ti, x))
                    end
                end
                function (::Type{$Ti})(x::$Tf)
                    if ($(Tf(typemin(Ti))) <= x <= $(Tf(typemax(Ti)))) && (round(x, RoundToZero) == x)
                        return unsafe_trunc($Ti,x)
                    else
                        throw(InexactError($(Expr(:quote,Ti.name.name)), $Ti, x))
                    end
                end
            end
        else
            # Here `eps(Tf(typemin(Ti))) > 1`, so the only value which can be truncated to
            # `Tf(typemin(Ti)` is itself. Similarly, `Tf(typemax(Ti))` is inexact and will
            # be rounded up. This assumes that `Tf(typemin(Ti)) > -Inf`, which is true for
            # these types, but not for `Float16` or larger integer types.
            @eval begin
                function trunc(::Type{$Ti},x::$Tf)
                    if $(Tf(typemin(Ti))) <= x < $(Tf(typemax(Ti)))
                        return unsafe_trunc($Ti,x)
                    else
                        throw(InexactError(:trunc, $Ti, x))
                    end
                end
                function (::Type{$Ti})(x::$Tf)
                    if ($(Tf(typemin(Ti))) <= x < $(Tf(typemax(Ti)))) && (round(x, RoundToZero) == x)
                        return unsafe_trunc($Ti,x)
                    else
                        throw(InexactError($(Expr(:quote,Ti.name.name)), $Ti, x))
                    end
                end
            end
        end
    end
end

"""
    issubnormal(f) -> Bool

Test whether a floating point number is subnormal.
"""
function issubnormal(x::T) where {T<:IEEEFloat}
    y = reinterpret(Unsigned, x)
    (y & exponent_mask(T) == 0) & (y & significand_mask(T) != 0)
end

ispow2(x::AbstractFloat) = !iszero(x) && frexp(x)[1] == 0.5

@eval begin
    typemin(::Type{Float16}) = $(bitcast(Float16, 0xfc00))
    typemax(::Type{Float16}) = $(Inf16)
    typemin(::Type{Float32}) = $(-Inf32)
    typemax(::Type{Float32}) = $(Inf32)
    typemin(::Type{Float64}) = $(-Inf64)
    typemax(::Type{Float64}) = $(Inf64)
    typemin(x::T) where {T<:Real} = typemin(T)
    typemax(x::T) where {T<:Real} = typemax(T)

    floatmin(::Type{Float16}) = $(bitcast(Float16, 0x0400))
    floatmin(::Type{Float32}) = $(bitcast(Float32, 0x00800000))
    floatmin(::Type{Float64}) = $(bitcast(Float64, 0x0010000000000000))
    floatmax(::Type{Float16}) = $(bitcast(Float16, 0x7bff))
    floatmax(::Type{Float32}) = $(bitcast(Float32, 0x7f7fffff))
    floatmax(::Type{Float64}) = $(bitcast(Float64, 0x7fefffffffffffff))

    eps(x::AbstractFloat) = isfinite(x) ? abs(x) >= floatmin(x) ? ldexp(eps(typeof(x)), exponent(x)) : nextfloat(zero(x)) : oftype(x, NaN)
    eps(::Type{Float16}) = $(bitcast(Float16, 0x1400))
    eps(::Type{Float32}) = $(bitcast(Float32, 0x34000000))
    eps(::Type{Float64}) = $(bitcast(Float64, 0x3cb0000000000000))
    eps() = eps(Float64)
end

"""
    floatmin(T = Float64)

Return the smallest positive normal number representable by the floating-point
type `T`.

# Examples
```jldoctest
julia> floatmin(Float16)
Float16(6.104e-5)

julia> floatmin(Float32)
1.1754944f-38

julia> floatmin()
2.2250738585072014e-308
```
"""
floatmin(x::T) where {T<:AbstractFloat} = floatmin(T)

"""
    floatmax(T = Float64)

Return the largest finite number representable by the floating-point type `T`.

# Examples
```jldoctest
julia> floatmax(Float16)
Float16(6.55e4)

julia> floatmax(Float32)
3.4028235f38

julia> floatmax()
1.7976931348623157e308
```
"""
floatmax(x::T) where {T<:AbstractFloat} = floatmax(T)

floatmin() = floatmin(Float64)
floatmax() = floatmax(Float64)

"""
    eps(::Type{T}) where T<:AbstractFloat
    eps()

Return the *machine epsilon* of the floating point type `T` (`T = Float64` by
default). This is defined as the gap between 1 and the next largest value representable by
`typeof(one(T))`, and is equivalent to `eps(one(T))`.  (Since `eps(T)` is a
bound on the *relative error* of `T`, it is a "dimensionless" quantity like [`one`](@ref).)

# Examples
```jldoctest
julia> eps()
2.220446049250313e-16

julia> eps(Float32)
1.1920929f-7

julia> 1.0 + eps()
1.0000000000000002

julia> 1.0 + eps()/2
1.0
```
"""
eps(::Type{<:AbstractFloat})

"""
    eps(x::AbstractFloat)

Return the *unit in last place* (ulp) of `x`. This is the distance between consecutive
representable floating point values at `x`. In most cases, if the distance on either side
of `x` is different, then the larger of the two is taken, that is

    eps(x) == max(x-prevfloat(x), nextfloat(x)-x)

The exceptions to this rule are the smallest and largest finite values
(e.g. `nextfloat(-Inf)` and `prevfloat(Inf)` for [`Float64`](@ref)), which round to the
smaller of the values.

The rationale for this behavior is that `eps` bounds the floating point rounding
error. Under the default `RoundNearest` rounding mode, if ``y`` is a real number and ``x``
is the nearest floating point number to ``y``, then

```math
|y-x| \\leq \\operatorname{eps}(x)/2.
```

# Examples
```jldoctest
julia> eps(1.0)
2.220446049250313e-16

julia> eps(prevfloat(2.0))
2.220446049250313e-16

julia> eps(2.0)
4.440892098500626e-16

julia> x = prevfloat(Inf)      # largest finite Float64
1.7976931348623157e308

julia> x + eps(x)/2            # rounds up
Inf

julia> x + prevfloat(eps(x)/2) # rounds down
1.7976931348623157e308
```
"""
eps(::AbstractFloat)

## byte order swaps for arbitrary-endianness serialization/deserialization ##
bswap(x::IEEEFloat) = bswap_int(x)

# integer size of float
uinttype(::Type{Float64}) = UInt64
uinttype(::Type{Float32}) = UInt32
uinttype(::Type{Float16}) = UInt16

## Array operations on floating point numbers ##

float(A::AbstractArray{<:AbstractFloat}) = A

function float(A::AbstractArray{T}) where T
    if !isconcretetype(T)
        error("`float` not defined on abstractly-typed arrays; please convert to a more specific type")
    end
    convert(AbstractArray{typeof(float(zero(T)))}, A)
end

float(r::StepRange) = float(r.start):float(r.step):float(last(r))
float(r::UnitRange) = float(r.start):float(last(r))
float(r::StepRangeLen{T}) where {T} =
    StepRangeLen{typeof(float(T(r.ref)))}(float(r.ref), float(r.step), length(r), r.offset)
function float(r::LinRange)
    LinRange(float(r.start), float(r.stop), length(r))
end

ext: jl

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
import Pkg
Pkg.add("LowRankApprox")
using LowRankApprox

m=128
n=64
M = matrixlib(:fourier, rand(m), rand(n))
A=real(M)
F = psvdfact(A)
@time U, S, V = psvd(A)

using LinearAlgebra
@time U1, S1, V1 = svd(A)

max_stars_repo_name: jamesjcai/ScTenifold

max_issues_repo_name: jamesjcai/ScTenifold

max_forks_repo_name: jamesjcai/ScTenifold


text:
using RobotZoo
using RobotDynamics
using StaticArrays
using BenchmarkTools
using FiniteDiff
using LinearAlgebra
using Test
using Random
Random.seed!(1)
const RD = RobotDynamics

## Double Integrator
D = 2
linmodel = RobotZoo.LinearModels.DoubleIntegrator(D)
model = RobotZoo.DoubleIntegrator(D)
x,u = rand(model)
@test RD.dynamics(model, x, u) ≈ RD.dynamics(linmodel, x, u)

h = 0.1
dmodel_exp = RobotZoo.DiscreteLinearModel(linmodel, h)
dmodel_lin = RobotZoo.LinearModels.DiscreteDoubleIntegrator(h, D)
@test RD.discrete_dynamics(dmodel_exp, x, u, 0.0, h) ≈
    RD.discrete_dynamics(dmodel_lin, x, u, 0.0, h)

## FlexibleSatellite
model = RobotZoo.LinearModels.FlexibleSatellite()
@test RD.dims(model) == (12,3,12)


text:
```
gascum(args...)  -->Plot

Plot the Gas Cum of a well in a single Plot

The next table show the list of variables allowed:

|PropertyName|Args|Default|Input|Description
|---|---|---|---|---|
|Date|Mandatory|--|Array{Date,1}| Array of dates|
|bsw|Mandatory|--|Array{Number,1}| Array of Gascum|
|Dlim|Optional|Drange=false|Dlim=[Date(y,m,d) Date(y,m,d)]|xlim Dates range|
|mrange|Optional|mrange=false|mrange=number| Time interval in months to set xticks|
|dtf|Optional|dft="u yy"|dtf=DateFormat|Date format of xticks|
```
@userplot GasCum

@recipe function f(h::GasCum; Dlim=false, mrange=6, dft="u yy")
                  legend --> false
                  x, y = h.args

    if Dlim==false
        r = x[1]:Dates.Month(mrange):x[end]
        d = Dates.format.(r,dft)
    else
        r = Dlim[1]:Dates.Month(mrange):Dlim[2]
        d = Dates.format.(r,dft)
    end

        @series begin
        seriestype := :path
        linestyle := :dash
        linewidth := 1
        seriescolor := RGB(0.6, 0.0, 0.0)
        ylabel := "Gas Cum [MMscf]"
        xticks := (r,d)

        if Dlim!=false
        xlim := (Dates.value(Dlim[1]),Dates.value(Dlim[2]))
        end
        x, y
        end
end


text:
# PNML id, id registry and related function.

"""
Holds a set of pnml id symbols and a lock to allow safe reentrancy.

$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct IDRegistry
    ids::Set{Symbol}
    lk::ReentrantLock
end

IDRegistry() = IDRegistry(Set{Symbol}(), ReentrantLock())

function Base.show(io::IO, reg::IDRegistry)
    print(io, typeof(reg), " ", length(reg.ids), " ids: ", reg.ids)
end

"""
$(TYPEDSIGNATURES)

Duplicated pnml node id has been detected.
Default `action` is to issue a warning.
"""
function duplicate_id_action(id::Symbol; action=:warn)
    action === :warn && @warn "ID '$id' already registered"
    action === :error && error("ID '$id' already registered")
    return nothing
end

"""
$(TYPEDSIGNATURES)

Register `id` symbol and return the symbol.
"""
register_id!(reg::IDRegistry, s::AbstractString) = register_id!(reg, Symbol(s))
function register_id!(reg::IDRegistry, id::Symbol)
    lock(reg.lk) do
        id ∈ reg.ids && duplicate_id_action(id)
        push!(reg.ids, id)
    end
    id
end

"""
$(TYPEDSIGNATURES)

Return `true` if `s` is registered in `reg`.
"""
isregistered(reg::IDRegistry, s::AbstractString) = isregistered(reg, Symbol(s))
function isregistered(reg::IDRegistry, id::Symbol)
    lock(reg.lk) do
        id ∈ reg.ids
    end
end

"""
$(TYPEDSIGNATURES)

Empty the set of id symbols. Use case is unit tests.
In normal use it should never be needed.
"""
function reset_registry!(reg::IDRegistry)
    lock(reg.lk) do
        empty!(reg.ids)
    end
end

function Base.isempty(reg::IDRegistry)
    lock(reg.lk) do
        isempty(reg.ids)
    end
end


text:
##===----------------------------------------------------------------------===##
##                                   ACCELERATE.JL                            ##
##                                                                            ##
##===----------------------------------------------------------------------===##
##                                                                            ##
##  This file provides access to Apple's Accelerate on OS X systems           ##
##                                                                            ##
##===----------------------------------------------------------------------===##
module Accelerate
import Vectorize: functions, addfunction

const libacc = "/System/Library/Frameworks/Accelerate.framework/Accelerate"

## This defines the mapping between Julia function names and Accelerate vecLib names.
## Each tuple is of the form (:julia_name, :accelerate_name).accelerate names should
## be stored without the "vv" prefix. The accelerate names are stored without type
## suffix; it is assumed that there is no suffix for Float64, and "f" suffix for Float32.
##
## This assumes arguments of the form (output_vector, input_vector, length)
##
const veclibfunctions =
    ((:ceil, :ceil), (:floor, :floor), (:sqrt, :sqrt), (:invsqrt, :rsqrt), (:rec, :rec),
     (:exp, :exp), (:exp2, :exp2), (:expm1, :expm1), (:log, :log), (:log1p, :log1p),
     (:log2, :log2), (:log10, :log10), (:sin, :sin), (:sinpi, :sinpi), (:cos, :cos),
     (:cospi, :cospi), (:tan, :tan), (:tanpi, :tanpi), (:asin, :asin), (:acos, :acos),
     (:atan, :atan), (:sinh, :sinh), (:cosh, :cosh), (:tanh, :tanh), (:asinh, :asinh),
     (:acosh, :acosh), (:atanh, :atanh), (:trunc,:int), (:round,:nint), (:exponent,:logb),
     (:abs,:fabs))

## For vecLibFunctions
for (T, suff) in ((Float64, ""), (Float32, "f"))

    for (f, fa) in veclibfunctions
        name = string(f)
        f! = Symbol("$(f)!")
        addfunction(functions, (f, (T,)), "Vectorize.Accelerate.$f")
        addfunction(functions, (f!, (T,T)), "Vectorize.Accelerate.$(f!)")
        @eval begin
            @doc """
            `$($f)(X::Vector{$($T)})`
            Implements element-wise **$($name)** over a **Vector{$($T)}**. Allocates
            memory to store result. *Returns:* **Vector{$($T)}**
            """ ->
            function ($f)(X::Vector{$T})
                out = similar(X)
                return ($f!)(out, X)
            end
            @doc """
            `$($f!)(result::Vector{$($T)}, X::Vector{$($T)})`
            Implements element-wise **$($name)** over a **Vector{$($T)}** and overwrites
            the result vector with computed value. *Returns:* **Vector{$($T)}** `result`
            """ ->
            function ($f!)(out::Vector{$T}, X::Vector{$T})
                ccall(($(string("vv", fa,suff)),libacc),Void,
                      (Ptr{$T},Ptr{$T},Ptr{Cint}),
                      out,X,&length(X))
                return out
            end
        end
    end
end

## This defines the mapping between Julia function names and Accelerate vDSP names.
## Each tuple is of the form (:julia_name, :accelerate_name); accelerate names should
## be stored without the "vDSP" prefix. The accelerate names are stored without type
## suffix; it is assumed that there is no suffix for Float32, and "D" suffix for Float64.
##
## This assumes arguments of the form (input_vec_a, stride_a, input_vec_b, stride_b,
##                                                   output_vec, stride_out, length)
##
const vDSPfunctions = ((:add, :vadd, "addition"), (:sub, :vsub, "subtraction"),
                       (:div, :vdiv, "division"), (:mul, :vmul, "multiplication"),
                       (:max, :vmax, "maximum"), (:min, :vmin, "minimum"))

for (T, suff) in ((Float64, "D"), (Float32, ""))

    for (f, fa, name) in vDSPfunctions
        f! = Symbol("$(f)!")
        addfunction(functions, (f, (T,T)), "Vectorize.Accelerate.$f")
        addfunction(functions, (f!, (T,T,T)), "Vectorize.Accelerate.$f!")
        @eval begin
            @doc """
            `$($f)(X::Vector{$($T)}, Y::Vector{$($T)})`
            Implements element-wise **$($name)** over two **Vector{$($T)}**. Allocates
            memory to store result. *Returns:* **Vector{$($T)}**
            """ ->
            function ($f)(X::Vector{$T}, Y::Vector{$T})
                out = Array($T, length(X))
                return ($f!)(out, Y, X)
            end
        end
        @eval begin
            @doc """
            `$($f!)(result::Vector{$($T)}, X::Vector{$($T)}, Y::Vector{$($T)})`
            Implements element-wise **$($name)** over two **Vector{$($T)}** and overwrites
            the result vector with computed value. *Returns:* **Vector{$($T)}** `result`
            """ ->
            function ($f!)(out::Vector{$T}, X::Vector{$T}, Y::Vector{$T})
                ccall(($(string("vDSP_", fa,suff)),libacc), Void,
                      (Ptr{$T},Int64, Ptr{$T}, Int64, Ptr{$T}, Int64, Int64),
                      X, 1, Y, 1, out, 1, length(out))
                return out
            end
        end
    end
end

## This defines the mapping between Julia function names and Accelerate vDSP names that
## return a scalar value.
## Each tuple is of the form (:julia_name, :accelerate_name); accelerate names should
## be stored without the "vDSP" prefix. The accelerate names are stored without type
## suffix; it is assumed that there is no suffix for Float32, and "D" suffix for Float64.
##
## This assumes arguments of the form (input_vec_a, stride_a,  output_scalar, length)
const vDSPscalar = ((:sum,  :sve, "sum"), (:summag, :svemg, "sum-of-magnitudes"),
                    (:sumsqr, :svesq, "sum-of-squares"), (:maximum, :maxv, "maximum value"),
                    (:minimum, :minv, "minimum value"),  (:mean, :meanv, "mean value"))

for (T, suff) in ((Float64, "D"), (Float32, ""))

    for (f, fa, name) in vDSPscalar
        addfunction(functions, (f, (T,)), "Vectorize.Accelerate.$f")
        @eval begin
            @doc """
            `$($f)(X::Vector{$($T)}, Y::Vector{$($T)})`
            Computes the **$($name)** of a **Vector{$($T)}**.
            *Returns:* **$($T)**
            """ ->
            function ($f)(X::Vector{$T})
                out = Ref{$T}(0.0)
                ccall(($(string("vDSP_", fa,suff)),libacc),Void,
                      (Ptr{$T},Cint, Ref{$T}, Cint),
                      X, 1, out, length(X))
                return out[]
            end
        end
    end

    for (f, fa, name) in ((:findmax, :maxvi, "maximum-value with index"),
                          (:findmin, :minvi, "minimum-value with index"))
        @eval begin
            @doc """
            `$($f)(X::Vector{$($T)})`
            Computes the **$($name)** element-wise over a **Vector{$($T)}**.
            *Returns:* **$(($T, UInt64))**
            """ ->
            function ($f)(X::Vector{$T})
                index = Ref{Int}(0)
                val = Ref{$T}(0.0)
                ccall(($(string("vDSP_", fa, suff), libacc)),  Void,
                      (Ptr{$T}, Int64,  Ref{$T}, Ref{Int}, UInt64),
                      X, 1, val, index, length(X))
                return (val[], index[]+1)
            end
        end
    end
end

end # module


text:
module ModuleReplicable

using Test
using Hyperspecialize

global A = Set{Type}([])

struct Weeble <: Real
  x::Int
end

f(::Real) = false

using Qux
import Qux.h

Qux.h(::Weeble, ::Real) = true

@testset "Module Replicable" begin

  # First, a test for module local widening

  @concretize TypicalTag []
  @replicable f(::@hyperspecialize(TypicalTag)) = true

  @test !f(1)
  @test !f(1.0)
  @test !Qux.f(1)
  @test !Qux.f(1.0)

  @widen TypicalTag Int
  @widen Qux.TypicalTag Float64

  @test f(1)
  @test !f(1.0)
  @test !Qux.f(1)
  @test Qux.f(1.0)

  # Do all methods corresponding to a type get widened?

  @test !Qux.g(1)
  @test Qux.g(1.0)

  # Can we resolve ambiguities?

  @test_throws MethodError h(Weeble(1), Wobble(2))

  @replicable h(::Weeble, ::@hyperspecialize(Real)) = true

  @test h(Weeble(1), Wobble(2))

end

end #module