diff --git a/src/higherorderfns.jl b/src/higherorderfns.jl
index 04166135..8ff81948 100644
--- a/src/higherorderfns.jl
+++ b/src/higherorderfns.jl
@@ -15,7 +15,8 @@ using ..SparseArrays: SparseVector, SparseMatrixCSC, FixedSparseCSC, SparseMatri
                       SparseMatrixCSCColumnSubset, SparseColumnView, SparseVectorView, SparseVectorPartialView,
                       indtype, fixed, move_fixed, nnz, nzrange, spzeros, ftranspose,
                       nonzeroinds, nonzeros, rowvals, getcolptr, getrowval, getnzval, widelength,
-                      _iszero, _isnotzero, _is_fixed, _checkbuffers, _densestructure!, _copyto_fixed!,
+                      _iszero, _isnotzero, _is_fixed, _checkbuffers, _densestructure!, _checkdensestructure,
+                      _copyto_fixed!, _densezero,
                       ReadOnly, _storagebuffers, _unaliasedcopy
 using Base.Broadcast: BroadcastStyle, Broadcasted, flatten
 using LinearAlgebra
@@ -197,13 +198,54 @@ _aliasesbuffers(Cbuffers, b) = isempty(b) ? any(c -> c === b, Cbuffers) : any(c
     # Avoid calculating f(zero) unless necessary as it may fail.
     let r = if _haszeros(A) && all(_haszeros, Bs)
             fofzeros = f(_zeros_eltypes(A, Bs...)...)
-            _iszero(fofzeros) ? _map_zeropres!(f, C, A, Bs...) :
-                _map_notzeropres!(f, fofzeros, C, A, Bs...)
+            _iszeroresult(fofzeros) ? _guardpattern(_map_zeropres!, f, C, A, Bs...) :
+                _guardfill(_map_notzeropres!, f, fofzeros, C, A, Bs...)
         else
-            _map_zeropres!(f, C, A, Bs...)
+            _guardpattern(_map_zeropres!, f, C, A, Bs...)
         end
     return r
 end
+# A zero-preserving kernel rebuilds the pattern of C as it goes. When `f` or the conversion
+# of its result throws part-way, a C that had the pattern of the result is left partly
+# written, as a dense destination is, and any other C structurally empty rather than with a
+# pattern that is half its own. A fixed C keeps its pattern throughout.
+function _guardpattern(kernel!::K, f::Tf, C::SparseVecOrMat, As::Vararg{SparseVecOrMat,N}) where {K,Tf,N}
+    try
+        return kernel!(f, C, As...)
+    catch
+        _is_fixed(C) || _keepspattern(C, As...) || _emptypattern!(C)
+        rethrow()
+    end
+end
+# The result of one argument of the shape of C has the pattern of that argument, so C held
+# it all along if it still does. The pattern of another result depends on its values.
+_keepspattern(C, A) = _aresameshape(C, A) && _samepattern(C, A)
+_keepspattern(C, As...) = false
+_emptypattern!(C) = (trimstorage!(C, 0); _finishempty!(C))
+# A kernel that fills C stores every entry of it first. What can be checked is checked before
+# that, the conversion of the fill value included, and leaves C as it was. A later failure
+# leaves a C that is not fixed structurally empty, as its values are no longer its own.
+function _guardfill(kernel!::K, f::Tf, fillvalue, C::SparseVecOrMat, As::Vararg{SparseVecOrMat,N}) where {K,Tf,N}
+    _checkdensestructure(C)
+    convert(eltype(C), fillvalue)
+    try
+        return kernel!(f, fillvalue, C, As...)
+    catch
+        _is_fixed(C) || _emptypattern!(C)
+        rethrow()
+    end
+end
+# Whether `f` at the zeros is a zero. A value of a type that has none, such as a `String`,
+# a `Tuple` or `nothing`, is not one, and fills the result like any other nonzero.
+_iszeroresult(x::Union{Number,AbstractArray,Missing}) = _iszero(x)
+function _iszeroresult(x)
+    try
+        return _iszero(x)
+    catch e
+        e isa MethodError || rethrow()
+        return false
+    end
+end
 
 function _noshapecheck_map(f::Tf, A::SparseVecOrMat, Bs::Vararg{SparseVecOrMat,N}) where {Tf,N}
     # Avoid calculating f(zero) unless necessary as it may fail.
@@ -211,7 +253,7 @@ function _noshapecheck_map(f::Tf, A::SparseVecOrMat, Bs::Vararg{SparseVecOrMat,N
     indextypeC = _promote_indtype(A, Bs...)
     r = if _haszeros(A) && all(_haszeros, Bs)
             fofzeros = f(_zeros_eltypes(A, Bs...)...)
-            fpreszeros = _iszero(fofzeros)
+            fpreszeros = _iszeroresult(fofzeros)
             maxnnzC = Int(fpreszeros ? min(widelength(A), _sumnnzs(A, Bs...)) : widelength(A))
             C = _allocres(size(A), indextypeC, entrytypeC, maxnnzC)
             fpreszeros ? _map_zeropres!(f, C, A, Bs...) :
@@ -232,13 +274,20 @@ copy(bc::SpBroadcasted1) = _asksthezeros(bc) && _fillszeros(bc.f, bc.args[1]) ?
     isempty(C) && return _finishempty!(C)
     f = bc.f
     fofnoargs = f()
-    if _iszero(fofnoargs) # f() is zero, so empty C
+    if _iszeroresult(fofnoargs) # f() is zero, so empty C
         trimstorage!(C, 0)
         _finishempty!(C)
     else # f() is nonzero, so densify C and fill with independent calls to f()
-        _densestructure!(C)
-        storedvals(C)[1] = fofnoargs
-        broadcast!(f, view(storedvals(C), 2:length(storedvals(C))))
+        _checkdensestructure(C)
+        firstvalue = convert(eltype(C), fofnoargs)
+        try
+            _densestructure!(C)
+            storedvals(C)[1] = firstvalue
+            broadcast!(f, view(storedvals(C), 2:length(storedvals(C))))
+        catch
+            _is_fixed(C) || _emptypattern!(C)
+            rethrow()
+        end
     end
     return _checkbuffers(C)
 end
@@ -254,7 +303,7 @@ function _diffshape_broadcast(f::Tf, A::SparseVecOrMat, Bs::Vararg{SparseVecOrMa
     # Avoid calculating f(zero) unless necessary as it may fail.
     r = if _haszeros(A) && all(_haszeros, Bs)
             fofzeros = f(_zeros_eltypes(A, Bs...)...)
-            _iszero(fofzeros) ?
+            _iszeroresult(fofzeros) ?
                 _broadcast_zeropres!(f, _allocres(shapeC, indextypeC, entrytypeC, maxnnzC), A, Bs...) :
                 _broadcast_notzeropres!(f, fofzeros, _allocres(shapeC, indextypeC, entrytypeC, _densennz(shapeC)), A, Bs...)
         else
@@ -328,12 +377,15 @@ Densifies `C`, storing `fillvalue` in place of each unstored entry in `A` and
 function _map_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseVecOrMat) where Tf
     # Build dense matrix structure in C, expanding storage if necessary
     _densestructure!(C)
+    # Below, and in the other `notzeropres` kernels, a value is left as filled only when it is
+    # identical to `fillvalue`: one that compares equal may be of another type or sign, and
+    # the comparison itself may be `missing`.
     # Populate values
     fill!(storedvals(C), fillvalue)
     @inbounds for (j, jo) in zip(columns(C), _densecoloffsets(C))
         for Ak in colrange(A, j)
             Cx = f(storedvals(A)[Ak])
-            Cx != fillvalue && (storedvals(C)[jo + storedinds(A)[Ak]] = Cx)
+            Cx !== fillvalue && (storedvals(C)[jo + storedinds(A)[Ak]] = Cx)
         end
     end
     # NOTE: Combining the fill! above into the loop above to avoid multiple sweeps over /
@@ -419,7 +471,7 @@ function _map_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseVecOrMa
                 Cx, Ci = f(zero(eltype(A)), storedvals(B)[Bk]), Bi
                 Bk += oneunit(Bk); Bi = Bk < stopBk ? W(storedinds(B)[Bk]) : rowsentinelB
             end
-            Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx)
+            Cx !== fillvalue && (storedvals(C)[jo + Ci] = Cx)
         end
     end
     return _checkbuffers(C)
@@ -475,7 +527,7 @@ function _map_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, As::Vararg{Spars
         while activerow < rowsentinel
             vals, ks, rows = _fusedupdate_all(rowsentinel, activerow, rows, ks, stopks, As)
             Cx = f(vals...)
-            Cx != fillvalue && (storedvals(C)[jo + activerow] = Cx)
+            Cx !== fillvalue && (storedvals(C)[jo + activerow] = Cx)
             activerow = min(rows...)
         end
     end
@@ -615,7 +667,7 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV
             bccolrangejA = numcols(A) == 1 ? colrange(A, 1) : colrange(A, j)
             for Ak in bccolrangejA
                 Cx, Ci = f(storedvals(A)[Ak]), storedinds(A)[Ak]
-                Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx)
+                Cx !== fillvalue && (storedvals(C)[jo + Ci] = Cx)
             end
         end
     # Cases with vertical expansion
@@ -625,7 +677,7 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV
             Ak, stopAk = numcols(A) == 1 ? (colstartind(A, 1), colboundind(A, 1)) : (colstartind(A, j), colboundind(A, j))
             Ax = Ak < stopAk ? svA[Ak] : zero(eltype(A))
             fofAx = f(Ax)
-            if fofAx != fillvalue
+            if fofAx !== fillvalue
                 for i in (jo + 1):(jo + numrows(C))
                     svC[i] = fofAx
                 end
@@ -842,7 +894,7 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV
                     Ak += oneunit(Ak); Ai = Ak < stopAk ? W(storedinds(A)[Ak]) : rowsentinelA
                     Bk += oneunit(Bk); Bi = Bk < stopBk ? W(storedinds(B)[Bk]) : rowsentinelB
                 end
-                Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx)
+                Cx !== fillvalue && (storedvals(C)[jo + Ci] = Cx)
             end
         end
     # Cases with vertical expansion
@@ -853,7 +905,7 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV
             Ax = Ak < stopAk ? storedvals(A)[Ak] : zero(eltype(A))
             Bx = Bk < stopBk ? storedvals(B)[Bk] : zero(eltype(B))
             Cx = f(Ax, Bx)
-            if Cx != fillvalue
+            if Cx !== fillvalue
                 for Ck::Int in (jo + 1):(jo + numrows(C))
                     storedvals(C)[Ck] = Cx
                 end
@@ -865,10 +917,10 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV
             Bk, stopBk = numcols(B) == 1 ? (colstartind(B, 1), colboundind(B, 1)) : (colstartind(B, j), colboundind(B, j))
             Ax = Ak < stopAk ? storedvals(A)[Ak] : zero(eltype(A))
             fvAzB = f(Ax, zero(eltype(B)))
-            if fvAzB == fillvalue
+            if fvAzB === fillvalue
                 while Bk < stopBk
                     Cx = f(Ax, storedvals(B)[Bk])
-                    Cx != fillvalue && (storedvals(C)[jo + storedinds(B)[Bk]] = Cx)
+                    Cx !== fillvalue && (storedvals(C)[jo + storedinds(B)[Bk]] = Cx)
                     Bk += oneunit(Bk)
                 end
             else
@@ -880,7 +932,7 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV
                     else
                         Cx = fvAzB
                     end
-                    Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx)
+                    Cx !== fillvalue && (storedvals(C)[jo + Ci] = Cx)
                 end
             end
         end
@@ -890,10 +942,10 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV
             Bk, stopBk = numcols(B) == 1 ? (colstartind(B, 1), colboundind(B, 1)) : (colstartind(B, j), colboundind(B, j))
             Bx = Bk < stopBk ? storedvals(B)[Bk] : zero(eltype(B))
             fzAvB = f(zero(eltype(A)), Bx)
-            if fzAvB == fillvalue
+            if fzAvB === fillvalue
                 while Ak < stopAk
                     Cx = f(storedvals(A)[Ak], Bx)
-                    Cx != fillvalue && (storedvals(C)[jo + storedinds(A)[Ak]] = Cx)
+                    Cx !== fillvalue && (storedvals(C)[jo + storedinds(A)[Ak]] = Cx)
                     Ak += oneunit(Ak)
                 end
             else
@@ -905,7 +957,7 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV
                     else
                         Cx = fzAvB
                     end
-                    Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx)
+                    Cx !== fillvalue && (storedvals(C)[jo + Ci] = Cx)
                 end
             end
         end
@@ -955,7 +1007,14 @@ _finishempty!(C::AbstractCompressedVector) = C
 _finishempty!(C::AbstractSparseMatrixCSC) = (fill!(getcolptr(C), 1); C)
 
 # special case - vector outer product
-_copy(f::typeof(*), x::SparseVectorOrView, y::AdjOrTransSparseVectorOrView) = _outer(x, y)
+_copy(f::typeof(*), x::SparseVectorOrView, y::AdjOrTransSparseVectorOrView) =
+    _outerkeepszeros(x, y) ? _outer(x, y) : _copy(f, _sparsifystructured(x), _sparsifystructured(y))
+# `_outer` multiplies stored entries with stored entries alone, which gives the product when
+# every stored entry times a zero is a zero. An `Inf` or a `NaN` times a zero is not.
+function _outerkeepszeros(x, y)
+    zx, zy, op = _densezero(x), _densezero(y), wrapperop(y)
+    return all(v -> _iszero(v * zy), nonzeros(x)) && all(v -> _iszero(zx * op(v)), nonzeros(parent(y)))
+end
 @inline _outer(x::SparseVectorOrView, y::Adjoint) = return _outer(conj, x, parent(y))
 @inline _outer(x::SparseVectorOrView, y::Transpose) = return _outer(identity, x, parent(y))
 function _outer(trans::Tf, x, y) where Tf
@@ -1079,11 +1138,11 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, As::Vararg
         rows = _initrowforcol_all(j, rowsentinel, isemptys, expandsverts, ks, As)
         defaultCx = f(defargs...)
         activerow = min(rows...)
-        if defaultCx == fillvalue # fillvalue-preserving column scan
+        if defaultCx === fillvalue # fillvalue-preserving column scan
             while activerow < rowsentinel
                 args, ks, rows = _fusedupdatebc_all(rowsentinel, activerow, rows, defargs, ks, stopks, As)
                 Cx = f(args...)
-                Cx != fillvalue && (storedvals(C)[jo + activerow] = Cx)
+                Cx !== fillvalue && (storedvals(C)[jo + activerow] = Cx)
                 activerow = min(rows...)
             end
         else # fillvalue-non-preserving column scan
@@ -1095,7 +1154,7 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, As::Vararg
                 else
                     Cx = defaultCx
                 end
-                Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx)
+                Cx !== fillvalue && (storedvals(C)[jo + Ci] = Cx)
             end
         end
     end
@@ -1192,15 +1251,14 @@ end
     Base.Broadcast.check_broadcast_axes(axes(dest), As...)
     _is_fixed(dest) && return _copyto_fixed_broadcast!(f, dest, As...)
     # Avoid calculating f(zero) unless necessary as it may fail.
-    all(_haszeros, As) || return _broadcast_zeropres!(f, dest, As...)
+    all(_haszeros, As) || return _guardpattern(_broadcast_zeropres!, f, dest, As...)
     fofzeros = f(_zeros_eltypes(As...)...)
-    if _iszero(fofzeros)
-        return _broadcast_zeropres!(f, dest, As...)
+    if _iszeroresult(fofzeros)
+        return _guardpattern(_broadcast_zeropres!, f, dest, As...)
     else
-        return _broadcast_notzeropres!(f, fofzeros, dest, As...)
+        return _guardfill(_broadcast_notzeropres!, f, fofzeros, dest, As...)
     end
 end
-
 # A fixed `dest` whose pattern the kernels cannot write as they go gets the result from a
 # writable array of its shape through `_copyto_fixed!`, which checks the pattern and the
 # values before it writes. Every nonzero of the result must lie in the pattern of `dest`,
@@ -1335,11 +1393,13 @@ _fillszeros(f, args...) = false
 # A function that throws at the zeros, or an element type without a zero, is left to the
 # sparse kernels, which evaluate the function there only if some entry is not stored.
 function _preserveszeros(f::F, args::Vararg{Any,N}) where {F,N}
+    local fofzeros
     try
-        return _iszero(f(_zeros_eltypes(args...)...))
+        fofzeros = f(_zeros_eltypes(args...)...)
     catch
         return true
     end
+    return _iszeroresult(fofzeros)
 end
 
 # The sparse kernels evaluate the function once at the zeros and then at the stored entries
@@ -1349,8 +1409,12 @@ end
 function _densenonzero(bc::Broadcasted, f::F, args::Vararg{Any,N}) where {F,N}
     T = Base.promote_typejoin_union(Base.promote_op(f, map(eltype, args)...))
     (isconcretetype(T) && T !== Bool) || return _densebroadcast(bc)
-    return _storedasdense(_copy(f, map(_sparsifystructured, args)...))
+    As = map(_sparsifystructured, args)
+    # the kernels store every entry, which the index type of the arguments may not address
+    _addresses(_promote_indtype(As...), map(length, axes(bc))) || return _densebroadcast(bc)
+    return _storedasdense(_copy(f, As...))
 end
+_addresses(::Type{Ti}, shape) where {Ti} = !(Ti <: Base.BitInteger) || prod(Int64, shape) < typemax(Ti)
 _storedasdense(C::Union{SparseVector,SparseMatrixCSC}) =
     length(getnzval(C)) == nnz(C) == length(C) ? reshape(getnzval(C), size(C)) : Array(C)
 _storedasdense(C) = Array(C)
@@ -1387,6 +1451,12 @@ _densifysparse(x) = x
 _sparsifystructured(M::AbstractMatrix) = _isdenselike(M) ? _fullystored(M) : SparseMatrixCSC(M)
 _sparsifystructured(V::AbstractVector) = _isdenselike(V) ? _fullystored(V) : SparseVector(V)
 _sparsifystructured(S::SparseVecOrMat) = S
+# stores the whole diagonal, as for the other banded types: a kernel that read a `-0.0` there
+# as a structural zero would lose its sign
+function _sparsifystructured(D::Diagonal)
+    n = size(D, 1)
+    return SparseMatrixCSC(n, n, collect(1:(n + 1)), collect(1:n), copyto!(Vector{eltype(D)}(undef, n), D.diag))
+end
 _sparsifystructured(S::SparseViewOfColumns) = copy(S)
 _sparsifystructured(x) = x
 
diff --git a/src/reductions.jl b/src/reductions.jl
index 3930d243..0cdb613d 100644
--- a/src/reductions.jl
+++ b/src/reductions.jl
@@ -61,17 +61,86 @@ end
 function Base._mapreduce(f::F, op::G, ::Base.IndexCartesian, A::SparseMatrixCSCOrColumnSubset{T}) where {F,G,T}
     z = nnz(A)
     n = widelength(A)
-    if z == 0
-        if n == 0
-            Base.mapreduce_empty(f, op, T)
-        else
-            _mapreducezeros(f, op, T, n-z-1, f(zero(T)))
-        end
+    if n == 0
+        Base.mapreduce_empty(f, op, T)
+    elseif !_commutes(f, op, T)
+        # every entry stored: they are folded in their order, without `f(0)`, which might throw
+        z == n ? Base._mapreduce(f, op, nzvalview(A)) : _mapreduce_ordered(f, op, A)
+    elseif z == 0
+        _mapreducezeros(f, op, T, n-z-1, f(zero(T)))
     else
         _mapreducezeros(f, op, T, n-z, Base._mapreduce(f, op, nzvalview(A)))
     end
 end
 
+# Whether `op` may fold the mapped values out of their order. It is taken to be associative,
+# and to be commutative only where it is known to be: a product commutes when its factors are
+# real or complex numbers, which is asked of inference.
+# `f` takes an element of each of the types given, one for each array that is reduced.
+_commutes(f, ::Union{typeof(+),typeof(Base.add_sum),typeof(max),typeof(min),typeof(|),typeof(&),
+                     typeof(Base.or_any),typeof(Base.and_all),typeof(Base._extrema_rf)}, ::Vararg{Type,N}) where {N} = true
+_commutes(f::F, ::Union{typeof(*),typeof(Base.mul_prod)}, Ts::Vararg{Type,N}) where {F,N} =
+    Base.promote_op(f, Ts...) <: Union{Real,Complex}
+_commutes(f, op, ::Vararg{Type,N}) where {N} = false
+
+# An `op` not known to commute is folded over the positions in their order. `A` has an
+# unstored entry, so `f(0)` is taken once. A run of unstored entries is folded into the stored
+# entry after it, which an associative `op` allows and which keeps the accumulator one call
+# of `op` away from its next value.
+function _mapreduce_ordered(f::F, op::G, A::SparseMatrixCSCOrColumnSubset{T}) where {F,G,T}
+    rf = Base.BottomRF(op)
+    rowval, nzval = getrowval(A), getnzval(A)
+    m = size(A, 1)
+    z, idem = _unstoredvalue(f, op, T)
+    v = Base._InitialValue()
+    unstored = Int64(0)   # unstored positions since the last stored one
+    @inbounds for col in axes(A, 2)
+        prev = 0
+        for k in getnzrange(A, col)
+            row = Int(rowval[k])
+            v = rf(v, _afterunstored(op, z, idem, unstored + (row - prev - 1), f(nzval[k])))
+            unstored, prev = Int64(0), row
+        end
+        unstored += m - prev
+    end
+    return unstored == 0 ? v : rf(v, _foldunstored(op, z, idem, unstored))
+end
+function _mapreduce_ordered(f::F, op::G, A::SparseVectorOrView{T}) where {F,G,T}
+    rf = Base.BottomRF(op)
+    nzind, nzval = nonzeroinds(A), nonzeros(A)
+    z, idem = _unstoredvalue(f, op, T)
+    v = Base._InitialValue()
+    prev = 0
+    @inbounds for k in eachindex(nzind, nzval)
+        i = Int(nzind[k])
+        v = rf(v, _afterunstored(op, z, idem, i - prev - 1, f(nzval[k])))
+        prev = i
+    end
+    n = length(A) - prev
+    return n == 0 ? v : rf(v, _foldunstored(op, z, idem, n))
+end
+# `f(0)` as `op` first sees it, and whether folding it with itself leaves it as it is, in
+# value and in type, which makes a run of any length that one value
+function _unstoredvalue(f::F, op::G, ::Type{T}) where {F,G,T}
+    z = Base.reduce_first(op, f(zero(T)))
+    zz = op(z, z)
+    return z, typeof(zz) === typeof(z) && isequal(zz, z)
+end
+# `z` folded over `n > 0` consecutive positions, in O(log(n)) calls of `op`
+@inline function _foldunstored(op::G, z, idem::Bool, n) where {G}
+    (idem || n == 1) && return z
+    p = r = z
+    n -= 1
+    while n > 0
+        isodd(n) && (r = op(r, p))
+        n >>= 1
+        n > 0 && (p = op(p, p))
+    end
+    return r
+end
+# the stored value `x` with the `n` unstored positions before it folded in
+@inline _afterunstored(op::G, z, idem::Bool, n, x) where {G} = n == 0 ? x : op(_foldunstored(op, z, idem, n), x)
+
 # Specialized mapreduce for +/*/min/max/_extrema_rf. The entries a sparse array does not
 # store are folded in at once, as `reduce_first` of `f(0)` (an `Int` for `sum` and `prod` of
 # small integers, as for dense input) times or to the power of their count, in the type the
@@ -93,10 +162,15 @@ _mapreducezeros(f::Base.ExtremaMap, op::typeof(Base._extrema_rf), ::Type{T}, nze
     nzeros == 0 ? v0 : op(v0, f(zero(T)))
 
 # Specialized mapreduce for any and all
+# Each value of the predicate is checked to be Boolean, as the dense reductions check it:
+# `|` and `&` would absorb another value next to a `missing`, or combine two of them.
+_boolean(v::Union{Bool,Missing}) = v
+_boolean(v) = v::Bool
+_booleanpredicate(f::F) where {F} = _boolean ∘ f
 Base._any(f, A::SparseMatrixCSCOrColumnSubset, ::Colon) =
-    iszero(widelength(A)) ? false : Base._mapreduce(f, |, IndexCartesian(), A)
+    iszero(widelength(A)) ? false : Base._mapreduce(_booleanpredicate(f), |, IndexCartesian(), A)
 Base._all(f, A::SparseMatrixCSCOrColumnSubset, ::Colon) =
-    iszero(widelength(A)) ? true  : Base._mapreduce(f, &, IndexCartesian(), A)
+    iszero(widelength(A)) ? true  : Base._mapreduce(_booleanpredicate(f), &, IndexCartesian(), A)
 
 function Base._mapreduce(f::F, op::Union{typeof(Base.mul_prod),typeof(*)}, ::Base.IndexCartesian, A::SparseMatrixCSCOrColumnSubset{T}) where {F,T}
     nnzA = nnz(A)
@@ -104,10 +178,16 @@ function Base._mapreduce(f::F, op::Union{typeof(Base.mul_prod),typeof(*)}, ::Bas
     if nzeros == 0
         # No zeros, so don't compute f(0) since it might throw
         Base._mapreduce(f, op, nzvalview(A))
+    elseif !_commutes(f, op, T)
+        _mapreduce_ordered(f, op, A)
     else
         v = Base.reduce_first(op, f(zero(T)))^nzeros
-        # Bail out early if initial reduction value is zero or if there are no stored elements
-        (_iszero(v) || nnzA == 0) ? v : v*Base._mapreduce(f, op, nzvalview(A))
+        (nnzA == 0 || (_iszero(v) && Base.promote_op(f, T) <: Integer)) && return v
+        r = v*Base._mapreduce(f, op, nzvalview(A))
+        # A zero ends the product unless a stored `Inf` or `NaN` makes it `NaN`. The product
+        # of the stored entries may also have overflowed to one; folding them into the zero
+        # instead cannot.
+        (_iszero(r) || !_iszero(v)) ? r : mapfoldl(f, op, nzvalview(A); init=v)
     end
 end
 
@@ -153,12 +233,6 @@ for T in (:SparseVectorOrView, :SparseMatrixCSCOrColumnSubset)
                                                          _mapreduce_ordered(f, op, init, A, B)
     end
 end
-# A product commutes when its factors are real or complex numbers, which is asked of inference:
-# the answer picks between two kernels that agree for such factors.
-_commutes(f, ::Union{typeof(+),typeof(Base.add_sum),typeof(max),typeof(min)}, ::Type, ::Type) = true
-_commutes(f::F, ::Union{typeof(*),typeof(Base.mul_prod)}, ::Type{Ta}, ::Type{Tb}) where {F,Ta,Tb} =
-    Base.promote_op(f, Ta, Tb) <: Union{Real,Complex}
-_commutes(f, op, ::Type, ::Type) = false
 
 # An `op` not known to commute is folded over the positions in their order. A zero is taken
 # only where an entry is not stored, as an array that stores every entry need not have one.
@@ -339,6 +413,8 @@ function _mapreduce_dim_sparse(f, op, init, A::SparseMatrixCSCOrColumnSubset{T,T
             getcolptr(R)[2] = 2
         end
         return R
+    elseif (rm == 1 || rn == 1) && !_inorder(f, op, A)
+        return _mapreduce_dim_sparse_ordered(f, op, init, A, dims, Tr, rm == 1)
     elseif rm == 1
         return _mapreducerows_sparse!(f, op, init, spzeros(Tr, Ti, 1, n), A)
     elseif rn == 1
@@ -349,6 +425,39 @@ function _mapreduce_dim_sparse(f, op, init, A::SparseMatrixCSCOrColumnSubset{T,T
     end
 end
 
+# The two kernels below fold each slice's stored entries first. For an `op` not known to
+# commute, the dense reduction, which keeps the order of the elements and takes time
+# proportional to nnz(A) + length(result) as well, is given the pattern they produce.
+function _mapreduce_dim_sparse_ordered(f, op, init, A::SparseMatrixCSCOrColumnSubset{T,Ti}, dims, ::Type{Tr},
+                                       alongrows::Bool) where {T,Ti,Tr}
+    m, n = size(A)
+    D = Base._mapreduce_dim(f, op, init, A, dims)
+    stores = fill(false, alongrows ? n : m)
+    if alongrows
+        @inbounds for col in 1:n
+            stores[col] = !isempty(getnzrange(A, col))
+        end
+    else
+        @inbounds for i in view(getrowval(A), _storedinds(A))
+            stores[i] = true
+        end
+    end
+    # a slice that stores nothing is kept when its reduction is not a zero
+    z = zero(Tr)
+    @inbounds for i in eachindex(stores)
+        stores[i] || (stores[i] = !isequal(D[i], z))
+    end
+    inds = convert(Vector{Ti}, findall(stores))
+    vals = Tr[D[i] for i in inds]
+    alongrows || return SparseMatrixCSC(m, 1, Ti[1, length(inds) + 1], inds, vals)
+    colptr = Vector{Ti}(undef, n + 1)
+    colptr[1] = 1
+    @inbounds for col in 1:n
+        colptr[col + 1] = colptr[col] + stores[col]
+    end
+    return SparseMatrixCSC(1, n, colptr, fill!(inds, 1), vals)
+end
+
 # `R` is a structurally empty `1 x n` sparse matrix: its columns are built in order
 function _mapreducerows_sparse!(f, op, init, R::SparseMatrixCSC, A::SparseMatrixCSCOrColumnSubset{T}) where T
     nzval = getnzval(A)
@@ -475,6 +584,44 @@ function _mapreducecols!(f, op, R::AbstractArray, A::SparseMatrixCSCOrColumnSubs
     R
 end
 
+# The kernels above fold each slice's stored entries first and its unstored ones last, which
+# is the order of the elements when `op` commutes or every entry is stored. Otherwise the
+# unstored entries are folded where they lie, with `f(0)` taken once.
+_inorder(f, op, A::AbstractArray{T}) where {T} = _commutes(f, op, T) || nnz(A) == length(A)
+function _mapreducerows_ordered!(f::F, op::G, R::AbstractArray, A::SparseMatrixCSCOrColumnSubset{T}) where {F,G,T}
+    require_one_based_indexing(A, R)
+    rowval, nzval = getrowval(A), getnzval(A)
+    m = size(A, 1)
+    z, idem = _unstoredvalue(f, op, T)
+    @inbounds for col in axes(A, 2)
+        r = R[1, col]
+        prev = 0
+        for j in getnzrange(A, col)
+            row = Int(rowval[j])
+            r = op(r, _afterunstored(op, z, idem, row - prev - 1, f(nzval[j])))
+            prev = row
+        end
+        R[1, col] = prev == m ? r : op(r, _foldunstored(op, z, idem, m - prev))
+    end
+    R
+end
+function _mapreducecols_ordered!(f::F, op::G, R::AbstractArray, A::SparseMatrixCSCOrColumnSubset{T}) where {F,G,T}
+    require_one_based_indexing(A, R)
+    rowval, nzval = getrowval(A), getnzval(A)
+    m, n = size(A)
+    z, idem = _unstoredvalue(f, op, T)
+    lastcol = zeros(Int, m)   # the last column folded into each row
+    @inbounds for col in axes(A, 2), j in getnzrange(A, col)
+        row = rowval[j]
+        R[row, 1] = op(R[row, 1], _afterunstored(op, z, idem, col - lastcol[row] - 1, f(nzval[j])))
+        lastcol[row] = col
+    end
+    @inbounds for i in 1:m
+        lastcol[i] == n || (R[i, 1] = op(R[i, 1], _foldunstored(op, z, idem, n - lastcol[i])))
+    end
+    R
+end
+
 function Base._mapreducedim!(f::F, op::G, R::AbstractArray, A::SparseMatrixCSCOrColumnSubset{T}) where {F,G,T}
     require_one_based_indexing(A, R)
     lsiz = Base.check_reducedims(R,A)
@@ -485,10 +632,10 @@ function Base._mapreducedim!(f::F, op::G, R::AbstractArray, A::SparseMatrixCSCOr
         R[1, 1] = op(R[1, 1], mapreduce(f, op, A))
     elseif size(R, 1) == 1
         # Reduction along rows
-        _mapreducerows!(f, op, R, A)
+        _inorder(f, op, A) ? _mapreducerows!(f, op, R, A) : _mapreducerows_ordered!(f, op, R, A)
     elseif size(R, 2) == 1
         # Reduction along columns
-        _mapreducecols!(f, op, R, A)
+        _inorder(f, op, A) ? _mapreducecols!(f, op, R, A) : _mapreducecols_ordered!(f, op, R, A)
     else
         # Reduction along a dimension > 2
         # Compute op(R, f(A))
@@ -582,7 +729,8 @@ function _mapreducecols!(f, op::typeof(+), R::AbstractArray, A::SparseMatrixCSCO
                 end
             end
             for i = 1:m
-                R[i, 1] += _addzeros(zeroval, rownz[i])
+                # a row stored in full has no `f(0)` to add, which need not be finite
+                rownz[i] == 0 || (R[i, 1] += _addzeros(zeroval, rownz[i]))
             end
         end
     end
@@ -590,37 +738,35 @@ function _mapreducecols!(f, op::typeof(+), R::AbstractArray, A::SparseMatrixCSCO
 end
 
 # any(pred, A, dims = 1) => mapreduce(pred, |, A, dims = 1)
-# final argument `post` is to allow post-mapping each columnar mapreduce
+# `R` holds the initial value of each column, which the column's reduction is folded into.
+# With `negate`, `pred` is the negated predicate of an `&` reduction, see below.
 function _mapreducerows!(pred::P, ::typeof(|), R::AbstractMatrix{Bool}, A::SparseMatrixCSCOrColumnSubset{Tv},
-                         post::F = identity) where {P, F, Tv}
+                         negate::Bool = false) where {P, Tv}
     nzval = getnzval(A)
     m, n = size(A)
     @inbounds for ii in axes(A,2)
+        # an initial value that already decides the column
+        R[1, ii] == negate || continue
         rng = getnzrange(A, ii)
         len = length(rng)
-        # An empty column is trivial
-        if len == 0
-            R[1, ii] = post(pred(zero(Tv)))
-            continue
-        end
-        # If predicate on zero is true, then sparse column can be short-circuited
-        if pred(zero(Tv)) && len < m
-            R[1, ii] = post(true)
-            continue
-        end
-        # Otherwise reduce over the stored values
-        r = false
-        for jj in rng
-            r = pred(nzval[jj])
-            r && break
+        if len < m && pred(zero(Tv))
+            # an unstored entry satisfies the predicate
+            r = true
+        else
+            # reduce over the stored values
+            r = false
+            for jj in rng
+                r = pred(nzval[jj])
+                r && break
+            end
         end
-        R[1, ii] = post(r)
+        R[1, ii] = r != negate
     end
     return R
 end
 # all(pred, A, dims = 1) => mapreduce(pred, &, A, dims = 1) == .!mapreduce(!pred, |, A, dims = 1)
 _mapreducerows!(pred::P, ::typeof(&), R::AbstractMatrix{Bool},
-                A::SparseMatrixCSCOrColumnSubset) where {P} = _mapreducerows!(!pred, |, R, A, !)
+                A::SparseMatrixCSCOrColumnSubset) where {P} = _mapreducerows!(!pred, |, R, A, true)
 
 # findmax/min and argmax/min methods
 # find the first entry of the sparse matrix that is not stored - return its cartesian index
@@ -724,8 +870,10 @@ function _findr(op, A::AbstractSparseMatrixCSC{Tv}, region) where {Tv}
     end
 end
 
-_isless_fm(a, b)    =  b == b && ( a != a || isless(a, b) )
-_isgreater_fm(a, b) =  b == b && ( a != a || isless(b, a) )
+# whether `a` replaces the best value so far `b`, by the test of Base's dense `findmin` and
+# `findmax`, so that the result agrees with dense on `NaN`, `missing` and signed zeros
+_isless_fm(a, b)    = Base.isgreater(b, a)
+_isgreater_fm(a, b) = isless(b, a)
 
 findmin(A::AbstractSparseMatrixCSC{Tv}, region::Union{Integer,Tuple{Integer},NTuple{2,Integer}}) where {Tv} =
     _findr(_isless_fm, A, region)
diff --git a/src/sparsevector.jl b/src/sparsevector.jl
index b56a3d0f..36d5e523 100644
--- a/src/sparsevector.jl
+++ b/src/sparsevector.jl
@@ -1516,6 +1516,7 @@ function Base._mapreduce(f::F, op::G, ::IndexCartesian, A::SparseVectorOrView) w
     T = eltype(A)
     isempty(A) && return Base.mapreduce_empty(f, op, T)
     z = nnz(A)
+    (z == length(A) || _commutes(f, op, T)) || return _mapreduce_ordered(f, op, A)
     rest, ini = if z == 0
         length(A)-z-1, f(zero(T))
     else
@@ -1525,9 +1526,9 @@ function Base._mapreduce(f::F, op::G, ::IndexCartesian, A::SparseVectorOrView) w
 end
 
 Base._any(f, A::SparseVectorOrView, ::Colon) =
-    iszero(length(A)) ? false : Base._mapreduce(f, |, IndexCartesian(), A)
+    iszero(length(A)) ? false : Base._mapreduce(_booleanpredicate(f), |, IndexCartesian(), A)
 Base._all(f, A::SparseVectorOrView, ::Colon) =
-    iszero(length(A)) ? true  : Base._mapreduce(f, &, IndexCartesian(), A)
+    iszero(length(A)) ? true  : Base._mapreduce(_booleanpredicate(f), &, IndexCartesian(), A)
 
 function Base.mapreducedim!(f::F, op::G, R::AbstractVector, A::SparseVectorOrView) where {F,G}
     isempty(A) && return R
@@ -1954,9 +1955,17 @@ _densifiable(A::SparseVecOrMat) = !_is_fixed(A) || nnz(A) == widelength(A)
 _checkdensifiable(A::SparseVecOrMat) = _densifiable(A) ||
     throw(ArgumentError("cannot store a nonzero f(0) into a $(nameof(typeof(A))), its sparsity pattern is read-only"))
 
+# What `_densestructure!` checks before it writes anything.
+_checkdensestructure(A::AbstractCompressedVector) = _checkdensifiable(A)
+function _checkdensestructure(A::AbstractSparseMatrixCSC)
+    _checkdensifiable(A)
+    _is_fixed(A) || convert(indtype(A), size(A, 1) * size(A, 2) + 1)   # the last column pointer
+    return true
+end
+
 # Store every position of `A`; the values beyond the former `nnz(A)` are left for the caller to fill.
 function _densestructure!(A::AbstractCompressedVector)
-    _checkdensifiable(A)
+    _checkdensestructure(A)
     _is_fixed(A) && return A   # passed the check above, so already full
     n = length(A)
     resize!(nonzeroinds(A), n)
@@ -1965,7 +1974,7 @@ function _densestructure!(A::AbstractCompressedVector)
     return A
 end
 function _densestructure!(A::AbstractSparseMatrixCSC)
-    _checkdensifiable(A)
+    _checkdensestructure(A)
     _is_fixed(A) && return A
     m, n = size(A)
     resize!(getrowval(A), m * n)
diff --git a/test/higherorderfns.jl b/test/higherorderfns.jl
index f64a71a7..68851196 100644
--- a/test/higherorderfns.jl
+++ b/test/higherorderfns.jl
@@ -1188,4 +1188,94 @@ end
     end
 end
 
+@static if COMPREHENSIVE
+@testset "map and broadcast of a function whose values are `missing`, have no zero or equal f(0) in another type" begin
+    A, Ad = fixturepair(Float64, 4, 3)
+    f = x -> x == 0 ? 1.0 : missing
+    @test isequal(map(f, A), map(f, Ad))
+    @test isequal(A .== missing, Ad .== missing)
+    g = (x, y) -> x == 0 ? 1.0 : missing
+    @test isequal(map(g, A, A), map(g, Ad, Ad))
+    @test isequal(broadcast(g, A, A[:, 1]), broadcast(g, Ad, Ad[:, 1]))
+    g3 = (x, y, z) -> x == 0 ? 1.0 : missing
+    @test isequal(broadcast(g3, A, A, A[:, 1]), broadcast(g3, Ad, Ad, Ad[:, 1]))
+    h = x -> x == 0 ? 1 : 1.0
+    @test typeof.(map(h, A)) == typeof.(map(h, Ad))
+    for k in (string, x -> (x,), x -> nothing, x -> x => 1)
+        @test map(k, A) == map(k, Ad)
+    end
+    @test string.(A) == string.(Ad) && string.(A) isa Matrix{String}
+    @test (A .=> 1) == (Ad .=> 1)
+end
+
+@testset "broadcast! of another shape into a fixed destination checks the result first" begin
+    F = SparseArrays.fixed(sparse([1.0 0; 1.0 2.0]))
+    Fd, nz = Array(F), copy(nonzeros(F))
+    @test_throws ArgumentError F .= F .+ sparsevec([3.0, 0.0])   # a nonzero at (1, 2)
+    @test nonzeros(F) == nz
+    @test_throws ArgumentError F .= F .+ [1.0 2.0]
+    @test nonzeros(F) == nz
+    Fi = SparseArrays.fixed(sparse([1 0; 1 2]))
+    @test_throws InexactError Fi .= Fi .* [0.5, 1.0]
+    @test nonzeros(Fi) == [1, 1, 2]
+    # a result that fits the pattern is stored, whatever zeros the kernel kept
+    @test (F .= F .* sparsevec([3.0, 0.0])) === F
+    @test F == Fd .* [3.0, 0.0] && rowvals(F) == [1, 2, 2]
+    copyto!(nonzeros(F), nz)
+    F .= F .* [3.0, 2.0]
+    @test F == Fd .* [3.0, 2.0]
+end
+
+@testset "broadcast! and map! leave a valid destination when `f` or a conversion throws" begin
+    A, _ = fixturepair(Float64, 4, 3)
+    top = maximum(nonzeros(A))   # not a zero: `f` is evaluated at the zeros before anything is written
+    @test top > 0
+    boom = (x, y) -> x == top ? error("boom") : x + y
+    for B in (copy(A), A[:, 1])   # the `map!` kernel, and the kernel for another shape
+        C = copy(A)
+        @test_throws ErrorException C .= boom.(A, B)
+        @test nnz(C) == 0 && iszero(Array(C))
+        C = copy(A)
+        @test_throws ErrorException C .= boom.(C, B)
+        @test nnz(C) == 0 && iszero(Array(C))
+    end
+    Ci = sparse([1 0; 0 2])
+    @test_throws InexactError Ci .= Ci .* sparse([0.5 1.0; 1.0 1.0])
+    @test nnz(Ci) == 0 && iszero(Array(Ci))
+    # one argument with the pattern of the destination: partly written, as a dense one is
+    C = sparse([1.0 0 0; 0 2.0 0; 0 0 3.0])
+    @test_throws ErrorException C .= (x -> x == 2 ? error("boom") : 2x).(C)
+    @test nnz(C) == 3 && mismatch(C, [2.0 0 0; 0 2.0 0; 0 0 3.0]) === nothing
+    # a function that fills the zeros: untouched when the fill value does not convert, and
+    # empty when `f` throws after every entry was stored
+    Ci = sparse([1 0; 0 2])
+    @test_throws InexactError Ci .= Ci .+ 0.5
+    @test nnz(Ci) == 2 && mismatch(Ci, [1 0; 0 2]) === nothing
+    for C in (sparse([1 0; 0 2]), sparse([1 0; 0 2]))
+        @test_throws ErrorException map!(x -> x == 2 ? error("boom") : x + 1, C, C)
+        @test nnz(C) == 0 && iszero(Array(C))
+    end
+    calls = Ref(0)
+    @test_throws ErrorException Ci .= (() -> (calls[] += 1) > 2 ? error("boom") : 1).()
+    @test nnz(Ci) == 0 && iszero(Array(Ci))
+    # every entry does not fit the index type: found before anything is written
+    S = sparse(Int8[1], Int8[1], [1.0], 12, 12)
+    S0 = copy(S)
+    @test_throws InexactError S .= 1.0
+    @test_throws InexactError map!(x -> x + 1, S, S0)
+    @test mismatch(S, Array(S0); Ti=Int8) === nothing
+    # a dense result does not go through that index type
+    @test cos.(S) == cos.(Array(S)) && cos.(S) isa Matrix{Float64}
+end
+
+@testset "broadcast keeps a non-finite value times an unstored zero, and the sign of a zero on a `Diagonal`" begin
+    for (x, y) in ((sparsevec([Inf, 0.0, 2.0]), sparsevec([0.0, 1.0, 3.0])),
+                   (sparsevec([1.0, 0.0, 2.0]), sparsevec([0.0, NaN, 3.0])),
+                   (sparsevec([Inf + 0im, 0, 2im]), sparsevec([0, 1im, 3.0 + 0im]))), w in (adjoint, transpose)
+        @test isequal(Array(x .* w(y)), Array(x) .* w(Array(y)))
+    end
+    @test isequal(Array(sparse(ones(2, 2)) ./ Diagonal([-0.0, 1.0])), ones(2, 2) ./ Diagonal([-0.0, 1.0]))
+end
+end
+
 end # module
diff --git a/test/reductions.jl b/test/reductions.jl
index 24f310da..af7b8fed 100644
--- a/test/reductions.jl
+++ b/test/reductions.jl
@@ -737,4 +737,89 @@ end
     end
 end
 
+@static if COMPREHENSIVE
+@testset "reductions agree with dense on non-finite values, `init`, the order of the elements and `missing`" begin
+    # a stored `Inf` or `NaN` is not absorbed by an unstored zero
+    for D in ([0.0 Inf; 1.0 2.0], [0.0 NaN; 1.0 2.0]), dims in (:, 1, 2, (1, 2))
+        @test isequal(prod(sparse(D); dims), prod(D; dims))
+    end
+    @test prod(sparse([0 3; 1 2])) === 0
+    # nor is an unstored zero lost when the product of the stored entries overflows
+    @test prod(sparse([0.0 1e200; 1e200 1.0])) === 0.0
+    @test prod(x -> x == 0 ? 0.0 : 1e200, sparse([0 1; 1 1])) === 0.0
+    # `f(0)` is not added to a row that stores every entry
+    f = x -> x == 0 ? Inf : x
+    D = [1.0 2.0; 0.0 3.0]
+    @test isequal(mapreduce(f, +, sparse(D); dims=2), mapreduce(f, +, D; dims=2))
+    # `|` and `&` along the rows fold into `init`
+    for (pred, op) in ((>(5), |), (!iszero, &)), init in (true, false),
+        D in ([1.0 2.0; 0.0 3.0], [1.0 2.0; 3.0 4.0], [0.0 6.0; 0.0 0.0])
+        @test mapreduce(pred, op, sparse(D); dims=1, init) == mapreduce(pred, op, D; dims=1, init)
+    end
+    # an `op` that is not known to commute folds the elements in their order
+    takelast = (a, b) -> b
+    A, Ad = fixturepair(Float64, 4, 3)
+    x, xd = fixturepair(Float64, 5)
+    B = sparse([1.0 2.0; 3.0 4.0])   # stores every entry
+    for (S, D) in ((A, Ad), (x, xd), (view(A, :, [3, 1]), Ad[:, [3, 1]]), (view(A, :, 2), Ad[:, 2]), (B, Array(B)))
+        @test mapreduce(v -> [v], vcat, S) == mapreduce(v -> [v], vcat, D)
+        @test mapreduce(string, *, S) == mapreduce(string, *, D)
+        @test reduce(takelast, S) === reduce(takelast, D)
+        S isa AbstractMatrix || continue
+        for dims in (1, 2, (1, 2))
+            rd = mapreduce(identity, takelast, D; dims, init=7.0)
+            @test mapreduce(identity, takelast, S; dims, init=7.0) == rd
+            @test mismatch(mapreduce(identity, takelast, S; dims, init=7.0, sparse=true), rd) === nothing
+            @test mapreduce(string, *, S; dims, init="") == mapreduce(string, *, D; dims, init="")
+        end
+    end
+    # the ordered fold takes `f(0)` once and each stored entry once
+    calls = Ref(0)
+    mapreduce(v -> (calls[] += 1; v), takelast, A)
+    @test calls[] == nnz(A) + 1
+    calls[] = 0
+    mapreduce(v -> (calls[] += 1; v), takelast, A; dims=1, init=7.0, sparse=true)
+    @test calls[] == nnz(A) + 1
+    # a run of unstored entries has the type `op` gives it, not that of `f(0)`
+    plus = (a, b) -> a + b
+    for S in (spzeros(2, 2), spzeros(4), sparse([0.0 0.0 0.0 1.0]))
+        @test mapreduce(!iszero, plus, S) === mapreduce(!iszero, plus, Array(S))
+    end
+    # `|` and `&` commute for one array and for two
+    @test SparseArrays._commutes(!iszero, |, Float64) && SparseArrays._commutes(>, &, Float64, Float64)
+    @test mapreduce(>, |, A, 2A) === mapreduce(>, |, Ad, 2Ad)
+    # a stored `-0.0` ranks below the zero of an unstored entry
+    for S in (sparse([1, 2], [1, 1], [-0.0, 0.0], 2, 2), sparse([2], [1], [-0.0], 2, 2), sparse([1], [2], [-0.0], 2, 2)),
+        fm in (findmin, findmax)
+        D = Array(S)
+        @test isequal(fm(S), fm(D))
+        @test isequal(fm(S; dims=1), fm(D; dims=1)) && isequal(fm(S; dims=2), fm(D; dims=2))
+    end
+    # `missing` in `findmin` and `findmax`, which dense ranks above `NaN`
+    for D in (Union{Missing,Float64}[1.0 missing; 0.0 2.0], Union{Missing,Float64}[NaN missing; 0.0 2.0],
+              Union{Missing,Float64}[0.0 NaN; -1.0 0.0]), fm in (findmin, findmax)
+        S = sparse(D)
+        @test isequal(fm(S), fm(D))
+        @test isequal(fm(S; dims=1), fm(D; dims=1)) && isequal(fm(S; dims=2), fm(D; dims=2))
+    end
+    # a predicate that is not Boolean is rejected, as for dense
+    for S in (sparse([1 0; 0 0]), sparsevec([1, 0]))
+        @test_throws TypeError any(S)
+        @test_throws TypeError all(S)
+    end
+    # also next to a `missing`, which `|` and `&` would let absorb the other value
+    for S in (sparse(reshape(Union{Missing,Int}[missing, 1], 2, 1)), sparse(Union{Missing,Int}[missing, 1]),
+              sparse(Union{Missing,Int}[missing 0; 1 2]))
+        @test_throws TypeError any(S)
+        @test_throws TypeError all(S)
+        @test_throws TypeError any(Array(S))
+        @test_throws TypeError all(Array(S))
+    end
+    M = sparse(Union{Missing,Bool}[missing false; false false])
+    @test any(M) === missing && all(M) === false
+    @test any(iszero, sparse([1.0 0; 2 3])) === true && all(iszero, sparse([1.0 0; 2 3])) === false
+    @test any(ismissing, M) === true && all(x -> x === false, sparsevec([false, false])) === true
+end
+end
+
 end # module