Generalise Zeros +/-

Project.toml

@@ -1,6 +1,6 @@
 name = "FillArrays"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "0.13.8"
+version = "0.13.9"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/FillArrays.jl

@@ -6,7 +6,7 @@ import Base: size, getindex, setindex!, IndexStyle, checkbounds, convert,
     +, -, *, /, \, diff, sum, cumsum, maximum, minimum, sort, sort!,
     any, all, axes, isone, iterate, unique, allunique, permutedims, inv,
     copy, vec, setindex!, count, ==, reshape, _throw_dmrs, map, zero,
-    show, view, in, mapreduce, one, reverse
+    show, view, in, mapreduce, one, reverse, promote_op
 
 import LinearAlgebra: rank, svdvals!, tril, triu, tril!, triu!, diag, transpose, adjoint, fill!,
     dot, norm2, norm1, normInf, normMinusInf, normp, lmul!, rmul!, diagzero, AdjointAbsVec, TransposeAbsVec,
@@ -208,35 +208,6 @@ sort(a::AbstractFill; kwds...) = a
 sort!(a::AbstractFill; kwds...) = a
 svdvals!(a::AbstractFillMatrix) = [getindex_value(a)*sqrt(prod(size(a))); Zeros(min(size(a)...)-1)]
 
-+(a::AbstractFill) = a
--(a::AbstractFill) = Fill(-getindex_value(a), size(a))
-
-# Fill +/- Fill
-function +(a::AbstractFill{T, N}, b::AbstractFill{V, N}) where {T, V, N}
-    axes(a) ≠ axes(b) && throw(DimensionMismatch("dimensions must match."))
-    return Fill(getindex_value(a) + getindex_value(b), axes(a))
-end
--(a::AbstractFill, b::AbstractFill) = a + (-b)
-
-function +(a::FillVector{T}, b::AbstractRange) where {T}
-    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
-    Tout = promote_type(T, eltype(b))
-    return a.value .+ b
-end
-+(a::AbstractRange, b::AbstractFill) = b + a
-# LinearAlgebra defines `+(a::UniformScaling, b::AbstractMatrix) = b + a`,
-# so the implementation of `+(a::AbstractFill{<:Any,2}, b::UniformScaling)` is sufficient
-function +(a::AbstractFillMatrix, b::UniformScaling)
-    n = LinearAlgebra.checksquare(a)
-    return a + Diagonal(Fill(b.λ, n))
-end
-
--(a::AbstractFill, b::AbstractRange) = a + (-b)
--(a::AbstractRange, b::AbstractFill) = a + (-b)
-# LinearAlgebra defines `-(a::AbstractMatrix, b::UniformScaling) = a + (-b)`,
-# so the implementation of `-(a::UniformScaling, b::AbstractFill{<:Any,2})` is sufficient
--(a::UniformScaling, b::AbstractFillMatrix) = a + (-b)
-
 function fill_reshape(parent, dims::Integer...)
     n = length(parent)
     prod(dims) == n || _throw_dmrs(n, "size", dims)

src/fillalgebra.jl

@@ -211,61 +211,48 @@ function dot(u::AbstractVector{T}, D::Diagonal{U,<:Zeros}, v::AbstractVector{V})
     zero(promote_type(T,U,V))
 end
 
-+(a::Zeros) = a
--(a::Zeros) = a
-
-# Zeros +/- Zeros
-function +(a::Zeros{T}, b::Zeros{V}) where {T, V}
-    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
-    return Zeros{promote_type(T,V)}(size(a)...)
+# Addition and Subtraction
+function +(a::Zeros{T}, b::Zeros{V}) where {T, V} # for disambiguity
+    promote_shape(a,b)
+    return elconvert(promote_op(+,T,V),a)
+end
+for TYPE in (:AbstractArray, :AbstractFill) # AbstractFill for disambiguity
+    @eval function +(a::$TYPE{T}, b::Zeros{V}) where {T, V}
+        promote_shape(a,b)
+        return elconvert(promote_op(+,T,V),a)
+    end
+    @eval +(a::Zeros, b::$TYPE) = b + a
 end
--(a::Zeros, b::Zeros) = -(a + b)
--(a::Ones, b::Ones) = Zeros(a)+Zeros(b)
 
-# Zeros +/- Fill and Fill +/- Zeros
-function +(a::AbstractFill{T}, b::Zeros{V}) where {T, V}
-    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
-    return AbstractFill{promote_type(T, V)}(a)
-end
-+(a::Zeros, b::AbstractFill) = b + a
--(a::AbstractFill, b::Zeros) = a + b
--(a::Zeros, b::AbstractFill) = a + (-b)
-
-# Zeros +/- Array and Array +/- Zeros
-function +(a::Zeros{T, N}, b::AbstractArray{V, N}) where {T, V, N}
-    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
-    return AbstractArray{promote_type(T,V),N}(b)
-end
-function +(a::Array{T, N}, b::Zeros{V, N}) where {T, V, N}
-    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
-    return AbstractArray{promote_type(T,V),N}(a)
+# for VERSION other than 1.6, could use ZerosMatrix only
+function +(a::AbstractFillMatrix{T}, b::UniformScaling) where {T}
+    n = checksquare(a)
+    return a + Diagonal(Fill(zero(T) + b.λ, n))
 end
 
-function -(a::Zeros{T, N}, b::AbstractArray{V, N}) where {T, V, N}
-    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
-    return -b + a
-end
--(a::Array{T, N}, b::Zeros{V, N}) where {T, V, N} = a + b
+# LinearAlgebra defines `-(a::AbstractMatrix, b::UniformScaling) = a + (-b)`,
+# so the implementation of `-(a::UniformScaling, b::AbstractFill{<:Any,2})` is sufficient
+-(a::UniformScaling, b::AbstractFill) = -b + a # @test I-Zeros(3,3) === Diagonal(Ones(3))
 
+-(a::Ones, b::Ones) = Zeros(a) + Zeros(b)
 
-+(a::AbstractRange, b::Zeros) = b + a
+# necessary for AbstractRange, Diagonal, etc
++(a::AbstractFill, b::AbstractFill) = fill_add(a, b)
++(a::AbstractFill, b::AbstractArray) = fill_add(b, a)
++(a::AbstractArray, b::AbstractFill) = fill_add(a, b)
+-(a::AbstractFill, b::AbstractFill) = a + (-b)
+-(a::AbstractFill, b::AbstractArray) = a + (-b)
+-(a::AbstractArray, b::AbstractFill) = a + (-b)
 
-function +(a::ZerosVector{T}, b::AbstractRange) where {T}
-    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
-    Tout = promote_type(T, eltype(b))
-    return Tout(first(b)):Tout(step(b)):Tout(last(b))
-end
-function +(a::ZerosVector{T}, b::UnitRange) where {T<:Integer}
-    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
-    Tout = promote_type(T, eltype(b))
-    return AbstractUnitRange{Tout}(b)
+@inline function fill_add(a, b::AbstractFill)
+    promote_shape(a, b)
+    a .+ getindex_value(b)
 end
 
-function -(a::ZerosVector, b::AbstractRange)
-    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
-    return -b + a
-end
--(a::AbstractRange, b::ZerosVector) = a + b
+# following needed since as of Julia v1.8 convert(AbstractArray{T}, ::AbstractRange) might return a Vector
+@inline elconvert(::Type{T}, A::AbstractRange) where T = T(first(A)):T(step(A)):T(last(A))
+@inline elconvert(::Type{T}, A::AbstractUnitRange) where T<:Integer = AbstractUnitRange{T}(A)
+@inline elconvert(::Type{T}, A::AbstractArray) where T = AbstractArray{T}(A)
 
 ####
 # norm

test/runtests.jl

@@ -146,9 +146,9 @@ include("infinitearrays.jl")
         y = x + x
         @test y isa Fill{Int,1}
         @test y[1] == 2
-        @test x + Zeros{Bool}(5) ≡ x
-        @test x - Zeros{Bool}(5) ≡ x
-        @test Zeros{Bool}(5) + x ≡ x
+        @test x + Zeros{Bool}(5) ≡ Ones{Int}(5)
+        @test x - Zeros{Bool}(5) ≡ Ones{Int}(5)
+        @test Zeros{Bool}(5) + x ≡ Ones{Int}(5)
         @test -x ≡ Fill(-1,5)
     end