diff --git a/Project.toml b/Project.toml
index 0d75506b..7438ca20 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "FillArrays"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "0.13.11"
+version = "1.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -9,7 +9,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-Aqua = "0.5"
+Aqua = "0.5, 0.6"
 julia = "1.6"
 
 [extras]
diff --git a/README.md b/README.md
index 997b3fec..86d7137e 100644
--- a/README.md
+++ b/README.md
@@ -14,12 +14,8 @@ as well as identity matrices.  This package exports the following types:
 
 
 The primary purpose of this package is to present a unified way of constructing
-matrices. For example, to construct a 5-by-5 `CLArray` of all zeros, one would use
-```julia
-julia> CLArray(Zeros(5,5))
-```
-Because `Zeros` is lazy, this can be accomplished on the GPU with no memory transfer.
-Similarly, to construct a 5-by-5 `BandedMatrix` of all zeros with bandwidths `(1,2)`, one would use  
+matrices. 
+For example, to construct a 5-by-5 `BandedMatrix` of all zeros with bandwidths `(1,2)`, one would use  
 ```julia
 julia> BandedMatrix(Zeros(5,5), (1, 2))
 ```
diff --git a/src/FillArrays.jl b/src/FillArrays.jl
index bf090490..e82141f8 100644
--- a/src/FillArrays.jl
+++ b/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, promote_op
+    show, view, in, mapreduce, one, reverse, promote_op, promote_rule
 
 import LinearAlgebra: rank, svdvals!, tril, triu, tril!, triu!, diag, transpose, adjoint, fill!,
     dot, norm2, norm1, normInf, normMinusInf, normp, lmul!, rmul!, diagzero, AdjointAbsVec, TransposeAbsVec,
@@ -18,7 +18,7 @@ import Base.Broadcast: broadcasted, DefaultArrayStyle, broadcast_shape
 import Statistics: mean, std, var, cov, cor
 
 
-export Zeros, Ones, Fill, Eye, Trues, Falses
+export Zeros, Ones, Fill, Eye, Trues, Falses, OneElement
 
 import Base: oneto
 
@@ -34,6 +34,7 @@ const AbstractFillVecOrMat{T} = Union{AbstractFillVector{T},AbstractFillMatrix{T
 
 ==(a::AbstractFill, b::AbstractFill) = axes(a) == axes(b) && getindex_value(a) == getindex_value(b)
 
+
 @inline function _fill_getindex(F::AbstractFill, kj::Integer...)
     @boundscheck checkbounds(F, kj...)
     getindex_value(F)
@@ -147,15 +148,47 @@ Fill{T,0}(x::T, ::Tuple{}) where T = Fill{T,0,Tuple{}}(x, ()) # ambiguity fix
 @inline axes(F::Fill) = F.axes
 @inline size(F::Fill) = map(length, F.axes)
 
+"""
+    getindex_value(F::AbstractFill)
+
+Return the value that `F` is filled with.
+
+# Examples
+
+```jldoctest
+julia> f = Ones(3);
+
+julia> FillArrays.getindex_value(f)
+1.0
+
+julia> g = Fill(2, 10);
+
+julia> FillArrays.getindex_value(g)
+2
+```
+"""
+getindex_value
+
 @inline getindex_value(F::Fill) = F.value
 
 AbstractArray{T}(F::Fill{V,N}) where {T,V,N} = Fill{T}(convert(T, F.value)::T, F.axes)
 AbstractArray{T,N}(F::Fill{V,N}) where {T,V,N} = Fill{T}(convert(T, F.value)::T, F.axes)
 AbstractFill{T}(F::AbstractFill) where T = AbstractArray{T}(F)
+AbstractFill{T,N}(F::AbstractFill) where {T,N} = AbstractArray{T,N}(F)
+AbstractFill{T,N,Ax}(F::AbstractFill{<:Any,N,Ax}) where {T,N,Ax} = AbstractArray{T,N}(F)
+
+convert(::Type{AbstractFill{T}}, F::AbstractFill) where T = convert(AbstractArray{T}, F)
+convert(::Type{AbstractFill{T,N}}, F::AbstractFill) where {T,N} = convert(AbstractArray{T,N}, F)
+convert(::Type{AbstractFill{T,N,Ax}}, F::AbstractFill{<:Any,N,Ax}) where {T,N,Ax} = convert(AbstractArray{T,N}, F)
 
 copy(F::Fill) = Fill(F.value, F.axes)
 
-""" Throws an error if `arr` does not contain one and only one unique value. """
+"""
+    unique_value(arr::AbstractArray)
+
+Return `only(unique(arr))` without intermediate allocations.
+Throws an error if `arr` does not contain one and only one unique value.
+"""
 function unique_value(arr::AbstractArray)
     if isempty(arr) error("Cannot convert empty array to Fill") end
     val = first(arr)
@@ -262,6 +295,7 @@ for (Typ, funcs, func) in ((:Zeros, :zeros, :zero), (:Ones, :ones, :one))
         @inline $Typ{T,N}(A::AbstractArray{V,N}) where{T,V,N} = $Typ{T,N}(size(A))
         @inline $Typ{T}(A::AbstractArray) where{T} = $Typ{T}(size(A))
         @inline $Typ(A::AbstractArray) = $Typ{eltype(A)}(A)
+        @inline $Typ(::Type{T}, m...) where T = $Typ{T}(m...)
 
         @inline axes(Z::$Typ) = Z.axes
         @inline size(Z::$Typ) = length.(Z.axes)
@@ -273,6 +307,14 @@ for (Typ, funcs, func) in ((:Zeros, :zeros, :zero), (:Ones, :ones, :one))
         copy(F::$Typ) = F
 
         getindex(F::$Typ{T,0}) where T = getindex_value(F)
+
+        promote_rule(::Type{$Typ{T, N, Axes}}, ::Type{$Typ{V, N, Axes}}) where {T,V,N,Axes} = $Typ{promote_type(T,V),N,Axes}
+        function convert(::Type{$Typ{T,N,Axes}}, A::$Typ{V,N,Axes}) where {T,V,N,Axes}
+            convert(T, getindex_value(A)) # checks that the types are convertible
+            $Typ{T,N,Axes}(axes(A))
+        end
+        convert(::Type{$Typ{T,N}}, A::$Typ{V,N,Axes}) where {T,V,N,Axes} = convert($Typ{T,N,Axes}, A)
+        convert(::Type{$Typ{T}}, A::$Typ{V,N,Axes}) where {T,V,N,Axes} = convert($Typ{T,N,Axes}, A)
     end
 end
 
@@ -284,6 +326,8 @@ for TYPE in (:Fill, :AbstractFill, :Ones, :Zeros), STYPE in (:AbstractArray, :Ab
     end
 end
 
+promote_rule(::Type{<:AbstractFill{T, N, Axes}}, ::Type{<:AbstractFill{V, N, Axes}}) where {T,V,N,Axes} = Fill{promote_type(T,V),N,Axes}
+
 """
     fillsimilar(a::AbstractFill, axes)
 
@@ -426,7 +470,8 @@ end
 
 
 ## Array
-Base.Array{T,N}(F::AbstractFill{V,N}) where {T,V,N} = fill(convert(T, getindex_value(F)), size(F))
+Base.Array{T,N}(F::AbstractFill{V,N}) where {T,V,N} =
+    convert(Array{T,N}, fill(convert(T, getindex_value(F)), size(F)))
 
 # These are in case `zeros` or `ones` are ever faster than `fill`
 for (Typ, funcs, func) in ((:Zeros, :zeros, :zero), (:Ones, :ones, :one))
@@ -437,7 +482,7 @@ end
 
 # temporary patch. should be a PR(#48895) to LinearAlgebra
 Diagonal{T}(A::AbstractFillMatrix) where T = Diagonal{T}(diag(A))
-function convert(::Type{T}, A::AbstractFillMatrix) where T<:Diagonal 
+function convert(::Type{T}, A::AbstractFillMatrix) where T<:Diagonal
     checksquare(A)
     isdiag(A) ? T(A) : throw(InexactError(:convert, T, A))
 end
@@ -496,14 +541,14 @@ sum(x::AbstractFill) = getindex_value(x)*length(x)
 sum(f, x::AbstractFill) = length(x) * f(getindex_value(x))
 sum(x::Zeros) = getindex_value(x)
 
-cumsum(x::AbstractFill{<:Any,1}) = range(getindex_value(x); step=getindex_value(x),
-                                                    length=length(x))
+# needed to support infinite case
+steprangelen(st...) = StepRangeLen(st...)
+cumsum(x::AbstractFill{<:Any,1}) = steprangelen(getindex_value(x), getindex_value(x), length(x))
 
 cumsum(x::ZerosVector) = x
 cumsum(x::ZerosVector{Bool}) = x
 cumsum(x::OnesVector{II}) where II<:Integer = convert(AbstractVector{II}, oneto(length(x)))
 cumsum(x::OnesVector{Bool}) = oneto(length(x))
-cumsum(x::AbstractFillVector{Bool}) = cumsum(AbstractFill{Int}(x))
 
 
 #########
@@ -717,4 +762,6 @@ Base.@propagate_inbounds function view(A::AbstractFill{<:Any,N}, I::Vararg{Real,
     fillsimilar(A)
 end
 
+include("oneelement.jl")
+
 end # module
diff --git a/src/fillalgebra.jl b/src/fillalgebra.jl
index 2dec1b61..0e35d1f9 100644
--- a/src/fillalgebra.jl
+++ b/src/fillalgebra.jl
@@ -86,7 +86,6 @@ end
 *(a::ZerosMatrix, b::AbstractMatrix) = mult_zeros(a, b)
 *(a::AbstractMatrix, b::ZerosVector) = mult_zeros(a, b)
 *(a::AbstractMatrix, b::ZerosMatrix) = mult_zeros(a, b)
-*(a::ZerosVector, b::AbstractVector) = mult_zeros(a, b)
 *(a::ZerosMatrix, b::AbstractVector) = mult_zeros(a, b)
 *(a::AbstractVector, b::ZerosMatrix) = mult_zeros(a, b)
 
@@ -160,38 +159,22 @@ function *(a::Transpose{T, <:AbstractVector{T}}, b::ZerosVector{T}) where T<:Rea
 end
 *(a::Transpose{T, <:AbstractMatrix{T}}, b::ZerosVector{T}) where T<:Real = mult_zeros(a, b)
 
-# treat zero separately to support ∞-vectors
-function _zero_dot(a, b)
-    axes(a) == axes(b) || throw(DimensionMismatch("dot product arguments have lengths $(length(a)) and $(length(b))"))
-    zero(promote_type(eltype(a),eltype(b)))
-end
-
-_fill_dot(a::Zeros, b::Zeros) = _zero_dot(a, b)
-_fill_dot(a::Zeros, b) = _zero_dot(a, b)
-_fill_dot(a, b::Zeros) = _zero_dot(a, b)
-_fill_dot(a::Zeros, b::AbstractFill) = _zero_dot(a, b)
-_fill_dot(a::AbstractFill, b::Zeros) = _zero_dot(a, b)
-
-function _fill_dot(a::AbstractFill, b::AbstractFill)
-    axes(a) == axes(b) || throw(DimensionMismatch("dot product arguments have lengths $(length(a)) and $(length(b))"))
-    getindex_value(a)getindex_value(b)*length(b)
-end
-
 # support types with fast sum
-function _fill_dot(a::AbstractFill, b)
+# infinite cases should be supported in InfiniteArrays.jl
+# type issues of Bool dot are ignored at present.
+function _fill_dot(a::AbstractFillVector{T}, b::AbstractVector{V}) where {T,V}
     axes(a) == axes(b) || throw(DimensionMismatch("dot product arguments have lengths $(length(a)) and $(length(b))"))
-    getindex_value(a)sum(b)
+    dot(getindex_value(a), sum(b))
 end
 
-function _fill_dot(a, b::AbstractFill)
+function _fill_dot_rev(a::AbstractVector{T}, b::AbstractFillVector{V}) where {T,V}
     axes(a) == axes(b) || throw(DimensionMismatch("dot product arguments have lengths $(length(a)) and $(length(b))"))
-    sum(a)getindex_value(b)
+    dot(sum(a), getindex_value(b))
 end
 
-
 dot(a::AbstractFillVector, b::AbstractFillVector) = _fill_dot(a, b)
 dot(a::AbstractFillVector, b::AbstractVector) = _fill_dot(a, b)
-dot(a::AbstractVector, b::AbstractFillVector) = _fill_dot(a, b)
+dot(a::AbstractVector, b::AbstractFillVector) = _fill_dot_rev(a, b)
 
 function dot(u::AbstractVector, E::Eye, v::AbstractVector)
     length(u) == size(E,1) && length(v) == size(E,2) ||
diff --git a/src/fillbroadcast.jl b/src/fillbroadcast.jl
index 574fb145..82ba41e2 100644
--- a/src/fillbroadcast.jl
+++ b/src/fillbroadcast.jl
@@ -102,10 +102,6 @@ _broadcasted_zeros(f, a, b) = Zeros{Base.Broadcast.combine_eltypes(f, (a, b))}(b
 _broadcasted_ones(f, a, b) = Ones{Base.Broadcast.combine_eltypes(f, (a, b))}(broadcast_shape(axes(a), axes(b)))
 _broadcasted_nan(f, a, b) = Fill(convert(Base.Broadcast.combine_eltypes(f, (a, b)), NaN), broadcast_shape(axes(a), axes(b)))
 
-# TODO: remove at next breaking version
-_broadcasted_zeros(a, b) = _broadcasted_zeros(+, a, b)
-_broadcasted_ones(a, b) = _broadcasted_ones(+, a, b)
-
 broadcasted(::DefaultArrayStyle, ::typeof(+), a::Zeros, b::Zeros) = _broadcasted_zeros(+, a, b)
 broadcasted(::DefaultArrayStyle, ::typeof(+), a::Ones, b::Zeros) = _broadcasted_ones(+, a, b)
 broadcasted(::DefaultArrayStyle, ::typeof(+), a::Zeros, b::Ones) = _broadcasted_ones(+, a, b)
@@ -247,3 +243,8 @@ broadcasted(::DefaultArrayStyle{N}, ::typeof(Base.literal_pow), ::Base.RefValue{
 broadcasted(::DefaultArrayStyle{N}, ::typeof(Base.literal_pow), ::Base.RefValue{typeof(^)}, r::Ones{T,N}, ::Base.RefValue{Val{k}}) where {T,N,k} = Ones{T}(axes(r))
 broadcasted(::DefaultArrayStyle{N}, ::typeof(Base.literal_pow), ::Base.RefValue{typeof(^)}, r::Zeros{T,N}, ::Base.RefValue{Val{0}}) where {T,N} = Ones{T}(axes(r))
 broadcasted(::DefaultArrayStyle{N}, ::typeof(Base.literal_pow), ::Base.RefValue{typeof(^)}, r::Zeros{T,N}, ::Base.RefValue{Val{k}}) where {T,N,k} = Zeros{T}(axes(r))
+
+# supports structured broadcast
+if isdefined(LinearAlgebra, :fzero)
+    LinearAlgebra.fzero(x::Zeros) = zero(eltype(x))
+end
\ No newline at end of file
diff --git a/src/oneelement.jl b/src/oneelement.jl
new file mode 100644
index 00000000..abd37f5e
--- /dev/null
+++ b/src/oneelement.jl
@@ -0,0 +1,51 @@
+"""
+    OneElement(val, ind, axesorsize) <: AbstractArray
+
+Represents an array with the specified axes (if its a tuple of `AbstractUnitRange`s)
+or size (if its a tuple of `Integer`s), with a single entry set to `val` and all others equal to zero,
+specified by `ind``.
+"""
+struct OneElement{T,N,I,A} <: AbstractArray{T,N}
+  val::T
+  ind::I
+  axes::A
+  OneElement(val::T, ind::I, axes::A) where {T<:Number, I<:NTuple{N,Int}, A<:NTuple{N,AbstractUnitRange}} where {N} = new{T,N,I,A}(val, ind, axes)
+end
+
+OneElement(val, inds::NTuple{N,Int}, sz::NTuple{N,Integer}) where N = OneElement(val, inds, oneto.(sz))
+"""
+    OneElement(val, ind::Int, n::Int)
+
+Creates a length `n` vector where the `ind` entry is equal to `val`, and all other entries are zero.
+"""
+OneElement(val, ind::Int, len::Int) = OneElement(val, (ind,), (len,))
+"""
+    OneElement(ind::Int, n::Int)
+
+Creates a length `n` vector where the `ind` entry is equal to `1`, and all other entries are zero.
+"""
+OneElement(inds::Int, sz::Int) = OneElement(1, inds, sz)
+OneElement{T}(val, inds::NTuple{N,Int}, sz::NTuple{N,Integer}) where {T,N} = OneElement(convert(T,val), inds, oneto.(sz))
+OneElement{T}(val, inds::Int, sz::Int) where T = OneElement{T}(val, (inds,), (sz,))
+
+"""
+    OneElement{T}(val, ind::Int, n::Int)
+
+Creates a length `n` vector where the `ind` entry is equal to `one(T)`, and all other entries are zero.
+"""
+OneElement{T}(inds::Int, sz::Int) where T = OneElement(one(T), inds, sz)
+
+Base.size(A::OneElement) = map(length, A.axes)
+Base.axes(A::OneElement) = A.axes
+function Base.getindex(A::OneElement{T,N}, kj::Vararg{Int,N}) where {T,N}
+    @boundscheck checkbounds(A, kj...)
+    ifelse(kj == A.ind, A.val, zero(T))
+end
+
+Base.replace_in_print_matrix(o::OneElement{<:Any,2}, k::Integer, j::Integer, s::AbstractString) =
+    o.ind == (k,j) ? s : Base.replace_with_centered_mark(s)
+
+function Base.setindex(A::Zeros{T,N}, v, kj::Vararg{Int,N}) where {T,N}
+    @boundscheck checkbounds(A, kj...)
+    OneElement(convert(T, v), kj, axes(A))
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 81ba61c5..808ec402 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -20,6 +20,7 @@ include("infinitearrays.jl")
 
             for T in (Int, Float64)
                 Z = $Typ{T}(5)
+                @test $Typ(T, 5) ≡ Z
                 @test eltype(Z) == T
                 @test Array(Z) == $funcs(T,5)
                 @test Array{T}(Z) == $funcs(T,5)
@@ -28,12 +29,14 @@ include("infinitearrays.jl")
                 @test convert(AbstractArray,Z) ≡ Z
                 @test convert(AbstractArray{T},Z) ≡ AbstractArray{T}(Z) ≡ Z
                 @test convert(AbstractVector{T},Z) ≡ AbstractVector{T}(Z) ≡ Z
+                @test convert(AbstractFill{T},Z) ≡ AbstractFill{T}(Z) ≡ Z
 
                 @test $Typ{T,1}(2ones(T,5)) == Z
                 @test $Typ{T}(2ones(T,5)) == Z
                 @test $Typ(2ones(T,5)) == Z
 
                 Z = $Typ{T}(5, 5)
+                @test $Typ(T, 5, 5) ≡ Z
                 @test eltype(Z) == T
                 @test Array(Z) == $funcs(T,5,5)
                 @test Array{T}(Z) == $funcs(T,5,5)
@@ -51,6 +54,9 @@ include("infinitearrays.jl")
                 @test convert(Fill{Int}, Ones(20)) ≡ Fill(1, 20)
                 @test convert(Fill{Int,1}, Ones(20)) ≡ Fill(1, 20)
                 @test convert(Fill{Int,1,Tuple{Base.OneTo{Int}}}, Ones(20)) ≡ Fill(1, 20)
+                @test convert(AbstractFill{Int}, Ones(20)) ≡ AbstractFill{Int}(Ones(20)) ≡ Ones{Int}(20)
+                @test convert(AbstractFill{Int,1}, Ones(20)) ≡ AbstractFill{Int,1}(Ones(20)) ≡ Ones{Int}(20)
+                @test convert(AbstractFill{Int,1,Tuple{Base.OneTo{Int}}}, Ones(20)) ≡ AbstractFill{Int,1,Tuple{Base.OneTo{Int}}}(Ones(20)) ≡ Ones{Int}(20)
 
                 @test $Typ{T,2}(2ones(T,5,5)) ≡ $Typ{T}(5,5)
                 @test $Typ{T}(2ones(T,5,5)) ≡ $Typ{T}(5,5)
@@ -194,6 +200,23 @@ include("infinitearrays.jl")
         A = FillArrays.RectDiagonal(Int[], (1:0, 1:4))
         @test !(0 in A)
     end
+
+    @testset "promotion" begin
+        Z = Zeros{Int}(5)
+        Zf = Zeros(5)
+        O = Ones{Int}(5)
+        Of = Ones{Float64}(5)
+        @test [Z,O] isa Vector{Fill{Int,1,Tuple{Base.OneTo{Int}}}}
+        @test [Z,Of] isa Vector{Fill{Float64,1,Tuple{Base.OneTo{Int}}}}
+        @test [O,O] isa Vector{Ones{Int,1,Tuple{Base.OneTo{Int}}}}
+        @test [O,Of] isa Vector{Ones{Float64,1,Tuple{Base.OneTo{Int}}}}
+        @test [Z,Zf] isa Vector{Zeros{Float64,1,Tuple{Base.OneTo{Int}}}}
+
+        @test convert(Ones{Int}, Of) ≡ convert(Ones{Int,1}, Of) ≡ convert(typeof(O), Of) ≡ O
+        @test convert(Zeros{Int}, Zf) ≡ convert(Zeros{Int,1}, Zf) ≡ convert(typeof(Z), Zf) ≡ Z
+
+        @test_throws MethodError convert(Zeros{SVector{2,Int}}, Zf)
+    end
 end
 
 @testset "indexing" begin
@@ -306,20 +329,32 @@ end
 # type, and produce numerically correct results.
 as_array(x::AbstractArray) = Array(x)
 as_array(x::UniformScaling) = x
-function test_addition_and_subtraction(As, Bs, Tout::Type)
+equal_or_undef(a::Number, b::Number) = (a == b) || isequal(a, b)
+equal_or_undef(a, b) = all(equal_or_undef.(a, b))
+function test_addition_subtraction_dot(As, Bs, Tout::Type)
     for A in As, B in Bs
-        @testset "$(typeof(A)) ± $(typeof(B))" begin
+        @testset "$(typeof(A)) and $(typeof(B))" begin
             @test A + B isa Tout{promote_type(eltype(A), eltype(B))}
-            @test as_array(A + B) == as_array(A) + as_array(B)
+            @test equal_or_undef(as_array(A + B), as_array(A) + as_array(B))
 
             @test A - B isa Tout{promote_type(eltype(A), eltype(B))}
-            @test as_array(A - B) == as_array(A) - as_array(B)
+            @test equal_or_undef(as_array(A - B), as_array(A) - as_array(B))
 
             @test B + A isa Tout{promote_type(eltype(B), eltype(A))}
-            @test as_array(B + A) == as_array(B) + as_array(A)
+            @test equal_or_undef(as_array(B + A), as_array(B) + as_array(A))
 
             @test B - A isa Tout{promote_type(eltype(B), eltype(A))}
-            @test as_array(B - A) == as_array(B) - as_array(A)
+            @test equal_or_undef(as_array(B - A), as_array(B) - as_array(A))
+            
+            # Julia 1.6 doesn't support dot(UniformScaling)
+            if VERSION < v"1.6.0" || VERSION >= v"1.8.0"
+                d1 = dot(A, B)
+                d2 = dot(as_array(A), as_array(B))
+                d3 = dot(B, A)
+                d4 = dot(as_array(B), as_array(A))
+                @test d1 ≈ d2 || d1 ≡ d2
+                @test d3 ≈ d4 || d3 ≡ d4
+            end
         end
     end
 end
@@ -349,24 +384,24 @@ end
     @test -A_fill === Fill(-A_fill.value, 5)
 
     # FillArray +/- FillArray should construct a new FillArray.
-    test_addition_and_subtraction((A_fill, B_fill), (A_fill, B_fill), Fill)
+    test_addition_subtraction_dot((A_fill, B_fill), (A_fill, B_fill), Fill)
     test_addition_and_subtraction_dim_mismatch(A_fill, Fill(randn(rng), 5, 2))
 
     # FillArray + Array (etc) should construct a new Array using `getindex`.
-    A_dense, B_dense = randn(rng, 5), [5, 4, 3, 2, 1]
-    test_addition_and_subtraction((A_fill, B_fill), (A_dense, B_dense), Array)
+    B_dense = (randn(rng, 5), [5, 4, 3, 2, 1], fill(Inf, 5), fill(NaN, 5))
+    test_addition_subtraction_dot((A_fill, B_fill), B_dense, Array)
     test_addition_and_subtraction_dim_mismatch(A_fill, randn(rng, 5, 2))
 
     # FillArray + StepLenRange / UnitRange (etc) should yield an AbstractRange.
     A_ur, B_ur = 1.0:5.0, 6:10
-    test_addition_and_subtraction((A_fill, B_fill), (A_ur, B_ur), AbstractRange)
+    test_addition_subtraction_dot((A_fill, B_fill), (A_ur, B_ur), AbstractRange)
     test_addition_and_subtraction_dim_mismatch(A_fill, 1.0:6.0)
     test_addition_and_subtraction_dim_mismatch(A_fill, 5:10)
 
     # FillArray + UniformScaling should yield a Matrix in general
     As_fill_square = (Fill(randn(rng, Float64), 3, 3), Fill(5, 4, 4))
     Bs_us = (UniformScaling(2.3), UniformScaling(3))
-    test_addition_and_subtraction(As_fill_square, Bs_us, Matrix)
+    test_addition_subtraction_dot(As_fill_square, Bs_us, Matrix)
     As_fill_nonsquare = (Fill(randn(rng, Float64), 3, 2), Fill(5, 3, 4))
     for A in As_fill_nonsquare, B in Bs_us
         test_addition_and_subtraction_dim_mismatch(A, B)
@@ -374,12 +409,12 @@ end
 
     # FillArray + StaticArray should not have ambiguities
     A_svec, B_svec = SVector{5}(rand(5)), SVector(1, 2, 3, 4, 5)
-    test_addition_and_subtraction((A_fill, B_fill, Zeros(5)), (A_svec, B_svec), SVector{5})
+    test_addition_subtraction_dot((A_fill, B_fill, Zeros(5)), (A_svec, B_svec), SVector{5})
 
     # Issue #224
     A_matmat, B_matmat = Fill(rand(3,3),5), [rand(3,3) for n=1:5]
-    test_addition_and_subtraction((A_matmat,), (A_matmat,), Fill)
-    test_addition_and_subtraction((B_matmat,), (A_matmat,), Vector)
+    test_addition_subtraction_dot((A_matmat,), (A_matmat,), Fill)
+    test_addition_subtraction_dot((B_matmat,), (A_matmat,), Vector)
 
     # Optimizations for Zeros and RectOrDiagonal{<:Any, <:AbstractFill}
     As_special_square = (
@@ -389,7 +424,7 @@ end
         RectDiagonal(Fill(randn(rng, Float64), 3), 3, 3), RectDiagonal(Fill(3, 4), 4, 4)
     )
     DiagonalAbstractFill{T} = Diagonal{T, <:AbstractFill{T, 1}}
-    test_addition_and_subtraction(As_special_square, Bs_us, DiagonalAbstractFill)
+    test_addition_subtraction_dot(As_special_square, Bs_us, DiagonalAbstractFill)
     As_special_nonsquare = (
         Zeros(3, 2), Zeros{Int}(3, 4),
         Eye(3, 2), Eye{Int}(3, 4),
@@ -508,13 +543,13 @@ end
     @test_throws MethodError [1,2,3]*Zeros(3) # Not defined for [1,2,3]*[0,0,0] either
 
     @testset "Check multiplication by Adjoint vectors works as expected." begin
-        @test randn(4, 3)' * Zeros(4) === Zeros(3)
-        @test randn(4)' * Zeros(4) === zero(Float64)
-        @test [1, 2, 3]' * Zeros{Int}(3) === zero(Int)
+        @test randn(4, 3)' * Zeros(4) ≡ Zeros(3)
+        @test randn(4)' * Zeros(4) ≡ transpose(randn(4)) * Zeros(4) ≡ zero(Float64)
+        @test [1, 2, 3]' * Zeros{Int}(3) ≡ zero(Int)
         @test [SVector(1,2)', SVector(2,3)', SVector(3,4)']' * Zeros{Int}(3) === SVector(0,0)
         @test_throws DimensionMismatch randn(4)' * Zeros(3)
         @test Zeros(5)' * randn(5,3) ≡ Zeros(5)'*Zeros(5,3) ≡ Zeros(5)'*Ones(5,3) ≡ Zeros(3)'
-        @test Zeros(5)' * randn(5) ≡ Zeros(5)' * Zeros(5) ≡ Zeros(5)' * Ones(5) ≡ 0.0
+        @test abs(Zeros(5)' * randn(5)) ≡ abs(Zeros(5)' * Zeros(5)) ≡ abs(Zeros(5)' * Ones(5)) ≡ 0.0
         @test Zeros(5) * Zeros(6)' ≡ Zeros(5,1) * Zeros(6)' ≡ Zeros(5,6)
         @test randn(5) * Zeros(6)' ≡ randn(5,1) * Zeros(6)' ≡ Zeros(5,6)
         @test Zeros(5) * randn(6)' ≡ Zeros(5,6)
@@ -529,10 +564,11 @@ end
         @test transpose([1, 2, 3]) * Zeros{Int}(3) === zero(Int)
         @test_throws DimensionMismatch transpose(randn(4)) * Zeros(3)
         @test transpose(Zeros(5)) * randn(5,3) ≡ transpose(Zeros(5))*Zeros(5,3) ≡ transpose(Zeros(5))*Ones(5,3) ≡ transpose(Zeros(3))
-        @test transpose(Zeros(5)) * randn(5) ≡ transpose(Zeros(5)) * Zeros(5) ≡ transpose(Zeros(5)) * Ones(5) ≡ 0.0
+        @test abs(transpose(Zeros(5)) * randn(5)) ≡ abs(transpose(Zeros(5)) * Zeros(5)) ≡ abs(transpose(Zeros(5)) * Ones(5)) ≡ 0.0
         @test randn(5) * transpose(Zeros(6)) ≡ randn(5,1) * transpose(Zeros(6)) ≡ Zeros(5,6)
         @test Zeros(5) * transpose(randn(6)) ≡ Zeros(5,6)
         @test transpose(randn(5)) * Zeros(5) ≡ 0.0
+        @test transpose(randn(5) .+ im) * Zeros(5) ≡ 0.0 + 0im
 
         @test transpose([[1,2]]) * Zeros{SVector{2,Int}}(1) ≡ 0
         @test_broken transpose([[1,2,3]]) * Zeros{SVector{2,Int}}(1)
@@ -547,13 +583,13 @@ end
         @test +(z1) === z1
         @test -(z1) === z1
 
-        test_addition_and_subtraction((z1, z2), (z1, z2), Zeros)
+        test_addition_subtraction_dot((z1, z2), (z1, z2), Zeros)
         test_addition_and_subtraction_dim_mismatch(z1, Zeros{Float64}(4, 2))
     end
 
     # `Zeros` +/- `Fill`s should yield `Fills`.
     fill1, fill2 = Fill(5.0, 4), Fill(5, 4)
-    test_addition_and_subtraction((z1, z2), (fill1, fill2), Fill)
+    test_addition_subtraction_dot((z1, z2), (fill1, fill2), Fill)
     test_addition_and_subtraction_dim_mismatch(z1, Fill(5, 5))
 
     X = randn(3, 5)
@@ -600,12 +636,13 @@ end
 
     @testset "Tests for ranges." begin
         X = randn(5)
-        @test !(Zeros(5) + X === X)
-        @test Zeros{Int}(5) + (1:5) === (1:5) && (1:5) + Zeros{Int}(5) === (1:5)
-        @test Zeros(5) + (1:5) === (1.0:1.0:5.0) && (1:5) + Zeros(5) === (1.0:1.0:5.0)
-        @test (1:5) - Zeros{Int}(5) === (1:5)
-        @test Zeros{Int}(5) - (1:5) === -1:-1:-5
-        @test Zeros(5) - (1:5) === -1.0:-1.0:-5.0
+        @test !(Zeros(5) + X ≡ X)
+        @test Zeros{Int}(5) + (1:5) ≡ (1:5) + Zeros{Int}(5) ≡ (1:5)
+        @test Zeros(5) + (1:5) ≡ (1:5) + Zeros(5) ≡ (1.0:1.0:5.0)
+        @test (1:5) - Zeros{Int}(5) ≡ (1:5)
+        @test Zeros{Int}(5) - (1:5) ≡ -1:-1:-5
+        @test Zeros(5) - (1:5) ≡ -1.0:-1.0:-5.0
+        @test Zeros{Int}(5) + (1.0:5) ≡ (1.0:5) + Zeros{Int}(5) ≡ 1.0:5
     end
 
     @testset "test Base.zero" begin
@@ -634,11 +671,11 @@ end
     @test sum(Fill(3,10)) ≡ 30
     @test reduce(+, Fill(3,10)) ≡ 30
     @test sum(x -> x + 1, Fill(3,10)) ≡ 40
-    @test cumsum(Fill(3,10)) ≡ 3:3:30
+    @test cumsum(Fill(3,10)) ≡ StepRangeLen(3,3,10)
 
     @test sum(Ones(10)) ≡ 10.0
     @test sum(x -> x + 1, Ones(10)) ≡ 20.0
-    @test cumsum(Ones(10)) ≡ 1.0:10.0
+    @test cumsum(Ones(10)) ≡ StepRangeLen(1.0, 1.0, 10)
 
     @test sum(Ones{Int}(10)) ≡ 10
     @test sum(x -> x + 1, Ones{Int}(10)) ≡ 20
@@ -654,7 +691,7 @@ end
 
     @test cumsum(Zeros{Bool}(10)) ≡ Zeros{Bool}(10)
     @test cumsum(Ones{Bool}(10)) ≡ Base.OneTo{Int}(10)
-    @test cumsum(Fill(true,10)) ≡ 1:1:10
+    @test cumsum(Fill(true,10)) ≡ StepRangeLen(true, true, 10)
 
     @test diff(Fill(1,10)) ≡ Zeros{Int}(9)
     @test diff(Ones{Float64}(10)) ≡ Zeros{Float64}(9)
@@ -759,6 +796,8 @@ end
         @test broadcast(*, rnge, Zeros(10, 10)) ≡ Zeros{Float64}(10, 10)
         @test broadcast(*, Ones{Int}(10), rnge) ≡ rnge
         @test broadcast(*, rnge, Ones{Int}(10)) ≡ rnge
+        @test broadcast(*, Ones(10), -5:4) ≡ broadcast(*, -5:4, Ones(10)) ≡ rnge
+        @test broadcast(*, Ones(10), -5:1:4) ≡ broadcast(*, -5:1:4, Ones(10)) ≡ rnge
         @test_throws DimensionMismatch broadcast(*, Fill(5.0, 11), rnge)
         @test broadcast(*, rnge, Fill(5.0, 10)) == broadcast(*, rnge, 5.0)
         @test_throws DimensionMismatch broadcast(*, rnge, Fill(5.0, 11))
@@ -1299,17 +1338,19 @@ end
     Random.seed!(5)
     u = rand(n)
     v = rand(n)
+    c = rand(ComplexF16, n)
 
     @test dot(u, D, v) == dot(u, v)
     @test dot(u, 2D, v) == 2dot(u, v)
     @test dot(u, Z, v) == 0
 
-    @test dot(Zeros(5), Zeros{ComplexF16}(5)) ≡ zero(ComplexF64)
-    @test dot(Zeros(5), Ones{ComplexF16}(5)) ≡ zero(ComplexF64)
-    @test dot(Ones{ComplexF16}(5), Zeros(5)) ≡ zero(ComplexF64)
-    @test dot(randn(5), Zeros{ComplexF16}(5)) ≡ dot(Zeros{ComplexF16}(5), randn(5)) ≡ zero(ComplexF64)
+    @test @inferred(dot(Zeros(5), Zeros{ComplexF16}(5))) ≡ zero(ComplexF64)
+    @test @inferred(dot(Zeros(5), Ones{ComplexF16}(5))) ≡ zero(ComplexF64)
+    @test abs(@inferred(dot(Ones{ComplexF16}(5), Zeros(5)))) ≡ abs(@inferred(dot(randn(5), Zeros{ComplexF16}(5)))) ≡ abs(@inferred(dot(Zeros{ComplexF16}(5), randn(5)))) ≡ zero(Float64) # 0.0 !≡ -0.0
+    @test @inferred(dot(c, Fill(1 + im, 15))) ≡ (@inferred(dot(Fill(1 + im, 15), c)))' ≈ @inferred(dot(c, fill(1 + im, 15)))
 
-    @test dot(Fill(1,5), Fill(2.0,5)) ≡ 10.0
+    @test @inferred(dot(Fill(1,5), Fill(2.0,5))) ≡ 10.0
+    @test_skip dot(Fill(true,5), Fill(Int8(1),5)) isa Int8 # not working at present
 
     let N = 2^big(1000) # fast dot for fast sum
         @test dot(Fill(2,N),1:N) == dot(Fill(2,N),1:N) == dot(1:N,Fill(2,N)) == 2*sum(1:N)
@@ -1338,6 +1379,8 @@ end
     @test stringmime("text/plain", Fill(7,2,3)) == "2×3 Fill{$Int}, with entries equal to 7"
     @test stringmime("text/plain", Fill(8.0,1)) == "1-element Fill{Float64}, with entry equal to 8.0"
     @test stringmime("text/plain", Eye(5)) == "5×5 Eye{Float64}"
+    # used downstream in LazyArrays.jl to deduce sparsity
+    @test Base.replace_in_print_matrix(Zeros(5,3), 1, 2, "0.0") == " ⋅ "
 
     # 2-arg show, compact printing
     @test repr(Zeros(3)) == "Zeros(3)"
@@ -1469,3 +1512,40 @@ end
     @test cor(Fill(3, 4, 5)) ≈ cor(fill(3, 4, 5)) nans=true
     @test cor(Fill(3, 4, 5), dims=2) ≈ cor(fill(3, 4, 5), dims=2) nans=true
 end
+
+@testset "Structured broadcast" begin
+    D = Diagonal(1:5)
+    @test D + Zeros(5,5) isa Diagonal
+    @test D - Zeros(5,5) isa Diagonal
+    @test D .+ Zeros(5,5) isa Diagonal
+    @test D .- Zeros(5,5) isa Diagonal
+    @test D .* Zeros(5,5) isa Diagonal
+    @test Zeros(5,5) .* D isa Diagonal
+    @test Zeros(5,5) - D isa Diagonal
+    @test Zeros(5,5) + D isa Diagonal
+    @test Zeros(5,5) .- D isa Diagonal
+    @test Zeros(5,5) .+ D isa Diagonal
+    f = (x,y) -> x+1
+    @test f.(D, Zeros(5,5)) isa Matrix
+end
+
+@testset "OneElement" begin
+    e₁ = OneElement(2, 5)
+    @test e₁ == [0,1,0,0,0]
+    @test_throws BoundsError e₁[6]
+
+    e₁ = OneElement{Float64}(2, 5)
+    @test e₁ == [0,1,0,0,0]
+
+    v = OneElement{Float64}(2, 3, 4)
+    @test v == [0,0,2,0]
+
+    V = OneElement(2, (2,3), (3,4))
+    @test V == [0 0 0 0; 0 0 2 0; 0 0 0 0]
+
+    @test stringmime("text/plain", V) == "3×4 OneElement{$Int, 2, Tuple{$Int, $Int}, Tuple{Base.OneTo{$Int}, Base.OneTo{$Int}}}:\n ⋅  ⋅  ⋅  ⋅\n ⋅  ⋅  2  ⋅\n ⋅  ⋅  ⋅  ⋅"
+
+    @test Base.setindex(Zeros(5), 2, 2) ≡ OneElement(2.0, 2, 5)
+    @test Base.setindex(Zeros(5,3), 2, 2, 3) ≡ OneElement(2.0, (2,3), (5,3))
+    @test_throws BoundsError Base.setindex(Zeros(5), 2, 6)
+end
\ No newline at end of file