Make Zero/One subtypes of Real

Project.toml

@@ -1,7 +1,7 @@
 name = "VectorInterface"
 uuid = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 authors = ["Jutho Haegeman <jutho.haegeman@ugent.be> and contributors"]
-version = "0.6.1"
+version = "0.7.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

docs/src/man/interface.md

@@ -154,10 +154,14 @@ One
 Zero
 ```
 
-`One` and `Zero` are singleton subtypes of `Number` used to represent hard-coded constant coefficients in linear combinations.
+`One` and `Zero` are singleton subtypes of `Real` used to represent hard-coded constant coefficients in linear combinations.
 They allow methods like [`add`](@ref) to dispatch on a unit coefficient at compile time, avoiding unnecessary multiplications.
 They are the default values for the `α` and `β` coefficients in [`add`](@ref), [`add!`](@ref), and [`add!!`](@ref).
 
+Both behave as ordinary real numbers: arithmetic and promotion with any `Number`, `conj`, `real`, `imag`, `abs`, `abs2`, `sign`, the `is*` predicates (`iszero`, `isone`, `isreal`, `isinteger`, `isfinite`, …), comparison against other reals, hashing consistent with `0` and `1`, and `isapprox`.
+Wherever an operation can, it returns `One()` or `Zero()` rather than a plain `1` or `0`, so that `α === One()` fast paths keep firing downstream.
+Deliberately unsupported operations are those with no meaning for a hard-coded coefficient, such as `widen`, `big` and `eps`; dividing by `Zero()` throws a `DivideError` rather than producing `Inf`.
+
 ## Supported types
 
 Out of the box, `VectorInterface.jl` provides implementations for:

src/onezero.jl

@@ -4,7 +4,7 @@
 Singleton type for representing a hard-coded constant 1 in vector addition / linear
 combinations.
 """
-struct One <: Number end
+struct One <: Real end
 
 """
     struct Zero end
@@ -16,7 +16,9 @@ combinations.
     Now I am become Zero, the destroyer of NaN. - Vishnu
 ```
 """
-struct Zero <: Number end
+struct Zero <: Real end
+
+const ZeroOne = Union{Zero, One}
 
 # Base Arithmetic
 # ---------------
@@ -59,10 +61,10 @@ Base.conj(::Zero) = Zero()
 
 Base.one(::Type{One}) = One()
 Base.one(::Type{Zero}) = One()
-Base.one(::Union{Zero, One}) = One()
+Base.one(::ZeroOne) = One()
 Base.zero(::Type{One}) = Zero()
 Base.zero(::Type{Zero}) = Zero()
-Base.zero(::Union{Zero, One}) = Zero()
+Base.zero(::ZeroOne) = Zero()
 
 Base.:(==)(::One, ::One) = true
 Base.:(==)(::Zero, ::Zero) = true
@@ -78,7 +80,84 @@ Base.promote_rule(::Type{Zero}, ::Type{T}) where {T <: Number} = T
 # disambiguate:
 Base.promote_rule(::Type{Bool}, ::Type{One}) = Bool
 Base.promote_rule(::Type{Bool}, ::Type{Zero}) = Bool
-Base.convert(::Type{T}, ::One) where {T <: Number} = one(T)
+
+# NOTE: deliberately no `convert(::Type{T}, ::Zero) where {T <: Number}` methods here.
+# Since `Zero`/`One` are `<: Real`, such methods supersede `convert(::Type{T}, x::Number)`
+# for real argument types and invalidate a large amount of precompiled Base code (~200
+# invalidations, mostly under `convert(::Type{Int64}, ::Real)`). Base's own
+# `convert(::Type{T}, x::Number) = T(x)::T` routes through the constructors below instead,
+# which is equivalent and invalidation-free. Do not reintroduce them.
 (T::Type{<:Number})(::One) = one(T)
-Base.convert(::Type{T}, ::Zero) where {T <: Number} = zero(T)
 (T::Type{<:Number})(::Zero) = zero(T)
+
+# Disambiguation
+# --------------
+# Being `<: Real` makes Base's Real-vs-Complex methods applicable, which collide with the
+# broad `::Number` methods above. The bodies match those methods; these only fix dispatch.
+for C in (:Complex, :(Complex{Bool}))
+    @eval begin
+        Base.:(+)(::Zero, x::$C) = x
+        Base.:(+)(x::$C, ::Zero) = x
+        Base.:(-)(::Zero, x::$C) = -x
+        Base.:(-)(x::$C, ::Zero) = x
+        Base.:(*)(::Zero, x::$C) = zero(x)
+        Base.:(*)(x::$C, ::Zero) = zero(x)
+        Base.:(*)(::One, x::$C) = x
+        Base.:(*)(x::$C, ::One) = x
+        Base.:(/)(::$C, ::Zero) = throw(DivideError())
+        Base.:(/)(x::$C, ::One) = x
+    end
+end
+Base.:(/)(::Zero, x::S) where {S <: Complex} = iszero(x) ? throw(DivideError()) : zero(x)
+Base.:(/)(::One, x::S) where {S <: Complex} = inv(x)
+
+# against `Bool(::Real)`, `Complex(::Real)` and `Complex{T}(::Real)`
+Base.Bool(::One) = true
+Base.Bool(::Zero) = false
+Base.Complex(::One) = Complex(true)
+Base.Complex(::Zero) = Complex(false)
+Base.Complex{T}(::One) where {T <: Real} = Complex{T}(one(T))
+Base.Complex{T}(::Zero) where {T <: Real} = Complex{T}(zero(T))
+
+# Utility
+# -------
+# `real`, `imag`, `abs`, `abs2`, `signbit`, `isreal`, `isnan` and mixed-type comparisons are
+# inherited from the `Real` fallbacks. What follows is what `Real` does not provide.
+
+Base.isinteger(::ZeroOne) = true
+Base.isfinite(::ZeroOne) = true
+Base.isinf(::ZeroOne) = false
+
+Base.iszero(::Zero) = true
+Base.iszero(::One) = false
+Base.isone(::Zero) = false
+Base.isone(::One) = true
+
+# `isequal(Zero(), 0)` holds through promotion, so the hashes have to agree as well
+Base.hash(::Zero, h::UInt) = hash(0, h)
+Base.hash(::One, h::UInt) = hash(1, h)
+
+# same-type comparisons: `<(x::T, y::T) where {T <: Real}` is a Base no-op error, and
+# `promote(Zero(), Zero())` is a fixed point, so these cannot be inherited
+Base.isless(::Zero, ::Zero) = false
+Base.isless(::One, ::One) = false
+Base.isless(::Zero, ::One) = true
+Base.isless(::One, ::Zero) = false
+Base.:(<)(::Zero, ::Zero) = false
+Base.:(<)(::One, ::One) = false
+Base.:(<)(::Zero, ::One) = true
+Base.:(<)(::One, ::Zero) = false
+Base.:(<=)(::Zero, ::Zero) = true
+Base.:(<=)(::One, ::One) = true
+Base.:(<=)(::Zero, ::One) = true
+Base.:(<=)(::One, ::Zero) = false
+
+# `sign(x::Real)` compares `x` against itself and hits the same no-op error
+Base.sign(::Zero) = Zero()
+Base.sign(::One) = One()
+
+# identity-preserving: the inherited versions promote through `Bool`
+Base.min(::Zero, ::One) = Zero()
+Base.min(::One, ::Zero) = Zero()
+Base.max(::Zero, ::One) = One()
+Base.max(::One, ::Zero) = One()

test/onezero.jl

@@ -1,5 +1,6 @@
 using Test
 using VectorInterface
+using LinearAlgebra
 
 const Z = Zero()
 const I = One()
@@ -19,6 +20,115 @@ const typelist = (
     @test isone(I) == true
     @test isone(Z) == false
     @test iszero(I) == false
+
+    @test isreal(Z)
+    @test isreal(I)
+    @test isinteger(Z)
+    @test isinteger(I)
+    @test isfinite(Z)
+    @test isfinite(I)
+    @test !isinf(Z)
+    @test !isinf(I)
+    @test !isnan(Z)
+    @test !isnan(I)
+    @test !signbit(Z)
+    @test !signbit(I)
+end
+
+@testset "traits" begin
+    @test @inferred(real(Z)) === Z
+    @test @inferred(real(I)) === I
+    @test @inferred(imag(Z)) === Z
+    @test @inferred(imag(I)) === Z
+    @test real(Zero) === Zero
+    @test real(One) === One
+
+    @test @inferred(abs(Z)) === Z
+    @test @inferred(abs(I)) === I
+    @test @inferred(abs2(Z)) === Z
+    @test @inferred(abs2(I)) === I
+    @test @inferred(sign(Z)) === Z
+    @test @inferred(sign(I)) === I
+
+    # inherited from `Real`/`Number`, pinned here so a Base change shows up as a failure
+    @test @inferred(conj(Z)) === Z
+    @test @inferred(conj(I)) === I
+    @test @inferred(adjoint(Z)) === Z
+    @test @inferred(adjoint(I)) === I
+    @test @inferred(transpose(Z)) === Z
+    @test @inferred(transpose(I)) === I
+    @test @inferred(oneunit(Z)) === I
+    @test @inferred(oneunit(I)) === I
+    @test Z^2 === Z  # not `@inferred`: `Zero()^n` is `Union{Zero, One}`, since `Zero()^0 == One()`
+    @test @inferred(I^2) === I
+    @test Z^0 === I
+    @test @inferred(float(Z)) === 0.0
+    @test @inferred(float(I)) === 1.0
+end
+
+@testset "ordering" begin
+    @test !(Z < Z)
+    @test !(I < I)
+    @test Z < I
+    @test !(I < Z)
+    @test Z <= Z
+    @test I <= I
+    @test Z <= I
+    @test !(I <= Z)
+
+    @test isless(Z, I)
+    @test !isless(I, Z)
+    @test !isless(Z, Z)
+    @test cmp(Z, I) == -1
+    @test cmp(I, Z) == 1
+
+    @test @inferred(min(Z, I)) === Z
+    @test @inferred(min(I, Z)) === Z
+    @test @inferred(max(Z, I)) === I
+    @test @inferred(max(I, Z)) === I
+
+    for T in typelist
+        T <: Real || continue
+        @test Z < one(T)
+        @test !(Z < zero(T))
+        @test Z <= zero(T)
+        @test one(T) > Z
+        @test zero(T) < I
+        @test isless(Z, one(T))
+        @test isless(zero(T), I)
+    end
+end
+
+@testset "hashing" begin
+    @test hash(Z) == hash(0)
+    @test hash(I) == hash(1)
+    @test isequal(Z, 0)
+    @test isequal(I, 1)
+    @test !isequal(Z, I)
+    # the reason the hashes have to agree in the first place
+    d = Dict(0 => :a, 1 => :b)
+    @test d[Z] === :a
+    @test d[I] === :b
+    @test length(Set([Z, 0, false])) == 1
+end
+
+@testset "isapprox" begin
+    @test Z ≈ Z
+    @test I ≈ I
+    @test !(Z ≈ I)
+
+    for T in typelist
+        @test Z ≈ zero(T)
+        @test I ≈ one(T)
+        @test !(Z ≈ one(T))
+        @test !(I ≈ zero(T))
+        @test zero(T) ≈ Z
+        @test one(T) ≈ I
+        if T <: AbstractFloat
+            @test isapprox(Z, T(1.0e-3); atol = 1.0e-1)
+            @test !isapprox(Z, T(1.0e-3); atol = 1.0e-6)
+        end
+    end
 end
 
 @testset "arithmetic" begin
@@ -91,3 +201,51 @@ end
         @test @inferred(convert(T, Z)) == zero(T)
     end
 end
+
+@testset "complex disambiguation" begin
+    # `<: Real` makes Base's Real-vs-Complex methods applicable; these cover the methods
+    # added to resolve the resulting ambiguities.
+    for x in (im, ComplexF32(1.0f0, 2.0f0), ComplexF64(1.0, 2.0), Complex{Int}(3, -4))
+        @test @inferred(Z + x) == x
+        @test @inferred(x + Z) == x
+        @test @inferred(Z - x) == -x
+        @test @inferred(x - Z) == x
+        @test @inferred(Z * x) == zero(x)
+        @test @inferred(x * Z) == zero(x)
+        @test @inferred(I * x) == x
+        @test @inferred(x * I) == x
+        @test @inferred(x / I) == x
+        @test @inferred(I / x) == inv(x)
+        @test @inferred(Z / x) == zero(x)
+        @test_throws DivideError x / Z
+    end
+
+    @test @inferred(Bool(Z)) === false
+    @test @inferred(Bool(I)) === true
+    @test @inferred(Complex(Z)) == 0
+    @test @inferred(Complex(I)) == 1
+    @test @inferred(Complex{Float64}(Z)) === ComplexF64(0.0)
+    @test @inferred(Complex{Float64}(I)) === ComplexF64(1.0)
+end
+
+@testset "LinearAlgebra scalars" begin
+    @test norm(Z) === 0.0
+    @test norm(I) === 1.0
+    @test @inferred(dot(Z, I)) === Z
+    @test @inferred(dot(I, I)) === I
+
+    # 5-arg `mul!` is what motivated widening the utility surface in the first place
+    for T in (Float64, ComplexF64)
+        A, B = rand(T, 3, 3), rand(T, 3, 3)
+        @test mul!(zeros(T, 3, 3), A, B, I, Z) ≈ A * B
+        C = rand(T, 3, 3)
+        @test mul!(copy(C), A, B, I, I) ≈ C + A * B
+
+        x, y = rand(T, 3), rand(T, 3)
+        @test mul!(zeros(T, 3), A, x, I, Z) ≈ A * x
+        @test axpy!(I, x, copy(y)) ≈ y + x
+        @test rmul!(copy(x), I) ≈ x
+        @test lmul!(I, copy(x)) ≈ x
+        @test rmul!(copy(x), Z) ≈ zero(x)
+    end
+end