From c6eacfdf45eb9622682e0603e78a21151e66d866 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Fri, 3 Jan 2025 17:35:31 -0500 Subject: [PATCH 1/5] evalpoly for matrix polynomials --- src/generic.jl | 80 +++++++++++++++++++++++++++++++++++++++++++++++++ test/generic.jl | 24 +++++++++++++++ 2 files changed, 104 insertions(+) diff --git a/src/generic.jl b/src/generic.jl index 2b03b249..ddab4cb6 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -2091,3 +2091,83 @@ end function copytrito!(B::StridedMatrixStride1{T}, A::StridedMatrixStride1{T}, uplo::AbstractChar) where {T<:BlasFloat} LAPACK.lacpy!(B, A, uplo) end + +# non-inplace fallback for evalpoly(X, p) +function _evalpoly(X::AbstractMatrix, p) + Base.require_one_based_indexing(p) + p0 = isempty(p) ? Base.reduce_empty_iter(+, p) : p[end] + Xone = one(X) + S = Base.promote_op(*, typeof(Xone), typeof(Xone))(Xone) * p0 + for i = length(p)-1:-1:1 + S = X * S + @inbounds(p[i] isa AbstractMatrix ? p[i] : p[i] * I) + end + return S +end + +_scalarval(x::Number) = x +_scalarval(x::UniformScaling) = x.λ + +""" + evalpoly!(Y::AbstractMatrix, X::AbstractMatrix, p) + +Evaluate the matrix polynomial ``Y = \\sum_k X^{k-1} p[k]``, storing the result +in-place in `Y`, for the coefficients `p[k]` (a vector or tuple). The coefficients +can be scalars, matrices, or [`UniformScaling`](@ref). + +Similar to `evalpoly`, but may be more efficient by working more in-place. (Some +allocations may still be required, however.) +""" +function evalpoly!(Y::AbstractMatrix, X::AbstractMatrix, p::Union{AbstractVector,Tuple}) + @boundscheck axes(Y,1) == axes(Y,2) == axes(X,1) == axes(X,2) + Base.require_one_based_indexing(p) + + N = length(p) + pN = iszero(N) ? Base.reduce_empty_iter(+, p) : p[N] + if pN isa AbstractMatrix + Y .= pN + elseif N > 1 && p[N-1] isa Union{Number,UniformScaling} + # initialize Y to p[N-1] I + X p[N], in-place + Y .= X .* _scalarval(pN) + for i in axes(Y,1) + @inbounds Y[i,i] += p[N-1] * I + end + N -= 1 + else + # initialize Y to one(Y) * pN in-place + for i in axes(Y,1) + for j in axes(Y,2) + @inbounds Y[i,j] = zero(Y[i,j]) + end + @inbounds Y[i,i] += one(Y[i,i]) * pN + end + end + if N > 1 + Z = similar(Y) # workspace for mul! + for i = N-1:-1:1 + mul!(Z, X, Y) + if p[i] isa AbstractMatrix + Y .= p[i] .+ Z + else + # Y = p[i] * I + Z, in-place + Y .= Z + for j in axes(Y,1) + @inbounds Y[j,j] += p[i] * I + end + end + end + end + return Y +end + +Base.evalpoly(X::AbstractMatrix, p::Tuple) = _evalpoly(X, p) +Base.evalpoly(X::AbstractMatrix, ::Tuple{}) = zero(one(X)) # dimensionless zero, i.e. 0 * x^0 +Base.evalpoly(X::AbstractMatrix, p::AbstractVector) = _evalpoly(X, p) + +Base.evalpoly(X::StridedMatrix{<:Number}, p::Tuple{Union{Number, UniformScaling}, Vararg{Union{Number, UniformScaling}}}) = + evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), typeof(_scalarval(p[begin])))), X, p) +Base.evalpoly(X::StridedMatrix{<:Number}, p::Tuple{AbstractMatrix{<:Number}, Vararg{AbstractMatrix{<:Number}}}) = + evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), eltype(p[begin]))), X, p) +Base.evalpoly(X::StridedMatrix{<:Number}, p::AbstractVector{<:Union{Number, UniformScaling}}) = + isempty(p) ? _evalpoly(X, p) : evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), typeof(_scalarval(p[begin])))), X, p) +Base.evalpoly(X::StridedMatrix{<:Number}, p::AbstractVector{<:AbstractMatrix{<:Number}}) = + isempty(p) ? _evalpoly(X, p) : evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), eltype(p[begin]))), X, p) diff --git a/test/generic.jl b/test/generic.jl index 6d11ec82..efa86adf 100644 --- a/test/generic.jl +++ b/test/generic.jl @@ -837,4 +837,28 @@ end end end +using LinearAlgebra: _evalpoly +naive_evalpoly(X, p) = sum(X^(i-1) * p[i] for i=1:length(p)) + +@testset "evalpoly" begin + for X in ([1 2 3;4 5 6;7 8 9], UpperTriangular([1 2 3;0 5 6;0 0 9]), SymTridiagonal([1,2,3],[4,5])) + @test @inferred(evalpoly(X, ())) == zero(X) == evalpoly(X, Int[]) + @test @inferred(evalpoly(X, (17,))) == one(X) * 17 + @test @inferred(_evalpoly(X, [1,2,3,4])) == @inferred(evalpoly(X, [1,2,3,4])) == + @inferred(evalpoly(X, (1,2,3,4))) == + naive_evalpoly(X, [1,2,3,4]) == 1*one(X) + 2*X + 3X^2 + 4X^3 + @test typeof(evalpoly(X, [1,2,3])) == typeof(evalpoly(X, (1,2,3))) == typeof(_evalpoly(X, [1,2,3])) == + typeof(X * X) + + for N in (1,2,4), p in (rand(-10:10, N), UniformScaling.(rand(-10:10, N)), [rand(-5:5,3,3) for _ = 1:N]) + @test _evalpoly(X, p) == evalpoly(X, p) == evalpoly(X, Tuple(p)) == naive_evalpoly(X, p) + end + for N in (1,2,4), p in (rand(N), UniformScaling.(rand(N)), [rand(3,3) for _ = 1:N]) + @test _evalpoly(X, p) ≈ evalpoly(X, p) ≈ evalpoly(X, Tuple(p)) ≈ naive_evalpoly(X, p) + end + + @test_throws MethodError evalpoly(X, []) + end +end + end # module TestGeneric From d5be83bdc9e0b5c2d1140a586d2d60e9f6e0badf Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Fri, 3 Jan 2025 17:41:29 -0500 Subject: [PATCH 2/5] tweak --- src/generic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/generic.jl b/src/generic.jl index ddab4cb6..808024ba 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -2160,7 +2160,7 @@ function evalpoly!(Y::AbstractMatrix, X::AbstractMatrix, p::Union{AbstractVector end Base.evalpoly(X::AbstractMatrix, p::Tuple) = _evalpoly(X, p) -Base.evalpoly(X::AbstractMatrix, ::Tuple{}) = zero(one(X)) # dimensionless zero, i.e. 0 * x^0 +Base.evalpoly(X::AbstractMatrix, ::Tuple{}) = zero(one(X)) # dimensionless zero, i.e. 0 * X^0 Base.evalpoly(X::AbstractMatrix, p::AbstractVector) = _evalpoly(X, p) Base.evalpoly(X::StridedMatrix{<:Number}, p::Tuple{Union{Number, UniformScaling}, Vararg{Union{Number, UniformScaling}}}) = From dbe178777e8556e10fff510f628f104449a4faaf Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Fri, 3 Jan 2025 18:24:34 -0500 Subject: [PATCH 3/5] fix load-order issues by moving evalpoly methods to a separate file --- src/LinearAlgebra.jl | 2 ++ src/evalpoly.jl | 81 ++++++++++++++++++++++++++++++++++++++++++++ src/generic.jl | 80 ------------------------------------------- 3 files changed, 83 insertions(+), 80 deletions(-) create mode 100644 src/evalpoly.jl diff --git a/src/LinearAlgebra.jl b/src/LinearAlgebra.jl index fc1081e0..2f5f79d9 100644 --- a/src/LinearAlgebra.jl +++ b/src/LinearAlgebra.jl @@ -561,6 +561,8 @@ include("schur.jl") include("structuredbroadcast.jl") include("deprecated.jl") +include("evalpoly.jl") + const ⋅ = dot const × = cross export ⋅, × diff --git a/src/evalpoly.jl b/src/evalpoly.jl new file mode 100644 index 00000000..b9fa6388 --- /dev/null +++ b/src/evalpoly.jl @@ -0,0 +1,81 @@ +# matrix methods for evalpoly(X, p) = ∑ₖ Xᵏ⁻¹ p[k] + +# non-inplace fallback for evalpoly(X, p) +function _evalpoly(X::AbstractMatrix, p) + Base.require_one_based_indexing(p) + p0 = isempty(p) ? Base.reduce_empty_iter(+, p) : p[end] + Xone = one(X) + S = Base.promote_op(*, typeof(Xone), typeof(Xone))(Xone) * p0 + for i = length(p)-1:-1:1 + S = X * S + @inbounds(p[i] isa AbstractMatrix ? p[i] : p[i] * I) + end + return S +end + +_scalarval(x::Number) = x +_scalarval(x::UniformScaling) = x.λ + +""" + evalpoly!(Y::AbstractMatrix, X::AbstractMatrix, p) + +Evaluate the matrix polynomial ``Y = \\sum_k X^{k-1} p[k]``, storing the result +in-place in `Y`, for the coefficients `p[k]` (a vector or tuple). The coefficients +can be scalars, matrices, or [`UniformScaling`](@ref). + +Similar to `evalpoly`, but may be more efficient by working more in-place. (Some +allocations may still be required, however.) +""" +function evalpoly!(Y::AbstractMatrix, X::AbstractMatrix, p::Union{AbstractVector,Tuple}) + @boundscheck axes(Y,1) == axes(Y,2) == axes(X,1) == axes(X,2) + Base.require_one_based_indexing(p) + + N = length(p) + pN = iszero(N) ? Base.reduce_empty_iter(+, p) : p[N] + if pN isa AbstractMatrix + Y .= pN + elseif N > 1 && p[N-1] isa Union{Number,UniformScaling} + # initialize Y to p[N-1] I + X p[N], in-place + Y .= X .* _scalarval(pN) + for i in axes(Y,1) + @inbounds Y[i,i] += p[N-1] * I + end + N -= 1 + else + # initialize Y to one(Y) * pN in-place + for i in axes(Y,1) + for j in axes(Y,2) + @inbounds Y[i,j] = zero(Y[i,j]) + end + @inbounds Y[i,i] += one(Y[i,i]) * pN + end + end + if N > 1 + Z = similar(Y) # workspace for mul! + for i = N-1:-1:1 + mul!(Z, X, Y) + if p[i] isa AbstractMatrix + Y .= p[i] .+ Z + else + # Y = p[i] * I + Z, in-place + Y .= Z + for j in axes(Y,1) + @inbounds Y[j,j] += p[i] * I + end + end + end + end + return Y +end + +Base.evalpoly(X::AbstractMatrix, p::Tuple) = _evalpoly(X, p) +Base.evalpoly(X::AbstractMatrix, ::Tuple{}) = zero(one(X)) # dimensionless zero, i.e. 0 * X^0 +Base.evalpoly(X::AbstractMatrix, p::AbstractVector) = _evalpoly(X, p) + +Base.evalpoly(X::StridedMatrix{<:Number}, p::Tuple{Union{Number, UniformScaling}, Vararg{Union{Number, UniformScaling}}}) = + evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), typeof(_scalarval(p[begin])))), X, p) +Base.evalpoly(X::StridedMatrix{<:Number}, p::Tuple{AbstractMatrix{<:Number}, Vararg{AbstractMatrix{<:Number}}}) = + evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), eltype(p[begin]))), X, p) +Base.evalpoly(X::StridedMatrix{<:Number}, p::AbstractVector{<:Union{Number, UniformScaling}}) = + isempty(p) ? _evalpoly(X, p) : evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), typeof(_scalarval(p[begin])))), X, p) +Base.evalpoly(X::StridedMatrix{<:Number}, p::AbstractVector{<:AbstractMatrix{<:Number}}) = + isempty(p) ? _evalpoly(X, p) : evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), eltype(p[begin]))), X, p) diff --git a/src/generic.jl b/src/generic.jl index 808024ba..2b03b249 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -2091,83 +2091,3 @@ end function copytrito!(B::StridedMatrixStride1{T}, A::StridedMatrixStride1{T}, uplo::AbstractChar) where {T<:BlasFloat} LAPACK.lacpy!(B, A, uplo) end - -# non-inplace fallback for evalpoly(X, p) -function _evalpoly(X::AbstractMatrix, p) - Base.require_one_based_indexing(p) - p0 = isempty(p) ? Base.reduce_empty_iter(+, p) : p[end] - Xone = one(X) - S = Base.promote_op(*, typeof(Xone), typeof(Xone))(Xone) * p0 - for i = length(p)-1:-1:1 - S = X * S + @inbounds(p[i] isa AbstractMatrix ? p[i] : p[i] * I) - end - return S -end - -_scalarval(x::Number) = x -_scalarval(x::UniformScaling) = x.λ - -""" - evalpoly!(Y::AbstractMatrix, X::AbstractMatrix, p) - -Evaluate the matrix polynomial ``Y = \\sum_k X^{k-1} p[k]``, storing the result -in-place in `Y`, for the coefficients `p[k]` (a vector or tuple). The coefficients -can be scalars, matrices, or [`UniformScaling`](@ref). - -Similar to `evalpoly`, but may be more efficient by working more in-place. (Some -allocations may still be required, however.) -""" -function evalpoly!(Y::AbstractMatrix, X::AbstractMatrix, p::Union{AbstractVector,Tuple}) - @boundscheck axes(Y,1) == axes(Y,2) == axes(X,1) == axes(X,2) - Base.require_one_based_indexing(p) - - N = length(p) - pN = iszero(N) ? Base.reduce_empty_iter(+, p) : p[N] - if pN isa AbstractMatrix - Y .= pN - elseif N > 1 && p[N-1] isa Union{Number,UniformScaling} - # initialize Y to p[N-1] I + X p[N], in-place - Y .= X .* _scalarval(pN) - for i in axes(Y,1) - @inbounds Y[i,i] += p[N-1] * I - end - N -= 1 - else - # initialize Y to one(Y) * pN in-place - for i in axes(Y,1) - for j in axes(Y,2) - @inbounds Y[i,j] = zero(Y[i,j]) - end - @inbounds Y[i,i] += one(Y[i,i]) * pN - end - end - if N > 1 - Z = similar(Y) # workspace for mul! - for i = N-1:-1:1 - mul!(Z, X, Y) - if p[i] isa AbstractMatrix - Y .= p[i] .+ Z - else - # Y = p[i] * I + Z, in-place - Y .= Z - for j in axes(Y,1) - @inbounds Y[j,j] += p[i] * I - end - end - end - end - return Y -end - -Base.evalpoly(X::AbstractMatrix, p::Tuple) = _evalpoly(X, p) -Base.evalpoly(X::AbstractMatrix, ::Tuple{}) = zero(one(X)) # dimensionless zero, i.e. 0 * X^0 -Base.evalpoly(X::AbstractMatrix, p::AbstractVector) = _evalpoly(X, p) - -Base.evalpoly(X::StridedMatrix{<:Number}, p::Tuple{Union{Number, UniformScaling}, Vararg{Union{Number, UniformScaling}}}) = - evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), typeof(_scalarval(p[begin])))), X, p) -Base.evalpoly(X::StridedMatrix{<:Number}, p::Tuple{AbstractMatrix{<:Number}, Vararg{AbstractMatrix{<:Number}}}) = - evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), eltype(p[begin]))), X, p) -Base.evalpoly(X::StridedMatrix{<:Number}, p::AbstractVector{<:Union{Number, UniformScaling}}) = - isempty(p) ? _evalpoly(X, p) : evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), typeof(_scalarval(p[begin])))), X, p) -Base.evalpoly(X::StridedMatrix{<:Number}, p::AbstractVector{<:AbstractMatrix{<:Number}}) = - isempty(p) ? _evalpoly(X, p) : evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), eltype(p[begin]))), X, p) From b43e5e6cad25fc4f141cbc2a977332c4cd712d5c Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Fri, 3 Jan 2025 18:59:47 -0500 Subject: [PATCH 4/5] fix heterogeneous and abstract cases --- src/evalpoly.jl | 16 ++++++++++------ test/generic.jl | 23 ++++++++++++++++------- 2 files changed, 26 insertions(+), 13 deletions(-) diff --git a/src/evalpoly.jl b/src/evalpoly.jl index b9fa6388..9639d605 100644 --- a/src/evalpoly.jl +++ b/src/evalpoly.jl @@ -67,15 +67,19 @@ function evalpoly!(Y::AbstractMatrix, X::AbstractMatrix, p::Union{AbstractVector return Y end +# fallback cases: call out-of-place _evalpoly Base.evalpoly(X::AbstractMatrix, p::Tuple) = _evalpoly(X, p) Base.evalpoly(X::AbstractMatrix, ::Tuple{}) = zero(one(X)) # dimensionless zero, i.e. 0 * X^0 Base.evalpoly(X::AbstractMatrix, p::AbstractVector) = _evalpoly(X, p) -Base.evalpoly(X::StridedMatrix{<:Number}, p::Tuple{Union{Number, UniformScaling}, Vararg{Union{Number, UniformScaling}}}) = - evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), typeof(_scalarval(p[begin])))), X, p) -Base.evalpoly(X::StridedMatrix{<:Number}, p::Tuple{AbstractMatrix{<:Number}, Vararg{AbstractMatrix{<:Number}}}) = - evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), eltype(p[begin]))), X, p) +# optimized in-place cases, limited to types like homogeneous tuples with length > 1 +# where we can reliably deduce the output type (= type of X * p[2]), +# and restricted to StridedMatrix (for now) so that we can be more confident that this is a performance win: +Base.evalpoly(X::StridedMatrix{<:Number}, p::Tuple{T, T, Vararg{T}}) where {T<:Union{Number, UniformScaling}} = + evalpoly!(similar(X, Base.promote_op(*, eltype(X), typeof(_scalarval(p[2])))), X, p) +Base.evalpoly(X::StridedMatrix{<:Number}, p::Tuple{AbstractMatrix{T}, AbstractMatrix{T}, Vararg{AbstractMatrix{T}}}) where {T<:Number} = + evalpoly!(similar(X, Base.promote_op(*, eltype(X), T)), X, p) Base.evalpoly(X::StridedMatrix{<:Number}, p::AbstractVector{<:Union{Number, UniformScaling}}) = - isempty(p) ? _evalpoly(X, p) : evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), typeof(_scalarval(p[begin])))), X, p) + length(p) < 2 ? _evalpoly(X, p) : evalpoly!(similar(X, Base.promote_op(*, eltype(X), typeof(_scalarval(p[begin+1])))), X, p) Base.evalpoly(X::StridedMatrix{<:Number}, p::AbstractVector{<:AbstractMatrix{<:Number}}) = - isempty(p) ? _evalpoly(X, p) : evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), eltype(p[begin]))), X, p) + length(p) < 2 ? _evalpoly(X, p) : evalpoly!(similar(X, Base.promote_op(*, eltype(X), eltype(p[begin+1]))), X, p) diff --git a/test/generic.jl b/test/generic.jl index efa86adf..49528b93 100644 --- a/test/generic.jl +++ b/test/generic.jl @@ -837,23 +837,32 @@ end end end -using LinearAlgebra: _evalpoly -naive_evalpoly(X, p) = sum(X^(i-1) * p[i] for i=1:length(p)) +using LinearAlgebra: _evalpoly # fallback routine, which we'll test explicitly + +# naive sum, a little complicated since X^0 fails if eltype(X) is abstract: +naive_evalpoly(X, p) = length(p) == 1 ? one(X) * p[1] : one(X) * p[1] + sum(X^(i-1) * p[i] for i=2:length(p)) @testset "evalpoly" begin - for X in ([1 2 3;4 5 6;7 8 9], UpperTriangular([1 2 3;0 5 6;0 0 9]), SymTridiagonal([1,2,3],[4,5])) + for X in ([1 2 3;4 5 6;7 8 9], UpperTriangular([1 2 3;0 5 6;0 0 9]), + SymTridiagonal([1,2,3],[4,5]), Real[1 2 3;4 5 6;7 8 9]) @test @inferred(evalpoly(X, ())) == zero(X) == evalpoly(X, Int[]) @test @inferred(evalpoly(X, (17,))) == one(X) * 17 - @test @inferred(_evalpoly(X, [1,2,3,4])) == @inferred(evalpoly(X, [1,2,3,4])) == - @inferred(evalpoly(X, (1,2,3,4))) == + @test _evalpoly(X, [1,2,3,4]) == evalpoly(X, [1,2,3,4]) == @inferred(evalpoly(X, (1,2,3,4))) == naive_evalpoly(X, [1,2,3,4]) == 1*one(X) + 2*X + 3X^2 + 4X^3 @test typeof(evalpoly(X, [1,2,3])) == typeof(evalpoly(X, (1,2,3))) == typeof(_evalpoly(X, [1,2,3])) == typeof(X * X) - for N in (1,2,4), p in (rand(-10:10, N), UniformScaling.(rand(-10:10, N)), [rand(-5:5,3,3) for _ = 1:N]) + # _evalpoly is not type-stable if eltype(X) is abstract + # because one(Real[...]) returns a Matrix{Int} + if isconcretetype(eltype(X)) + @inferred evalpoly(X, [1,2,3,4]) + @inferred _evalpoly(X, [1,2,3,4]) + end + + for N in (1,2,4), p in (Real[1,2], rand(-10:10, N), UniformScaling.(rand(-10:10, N)), [rand(-5:5,3,3) for _ = 1:N]) @test _evalpoly(X, p) == evalpoly(X, p) == evalpoly(X, Tuple(p)) == naive_evalpoly(X, p) end - for N in (1,2,4), p in (rand(N), UniformScaling.(rand(N)), [rand(3,3) for _ = 1:N]) + for N in (1,2,4), p in ((5, 6.7), rand(N), UniformScaling.(rand(N)), [rand(3,3) for _ = 1:N]) @test _evalpoly(X, p) ≈ evalpoly(X, p) ≈ evalpoly(X, Tuple(p)) ≈ naive_evalpoly(X, p) end From 1b43a2e26974227cf7926457db23c104a1c66cd0 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Fri, 10 Jan 2025 09:51:31 -0500 Subject: [PATCH 5/5] Update src/evalpoly.jl Co-authored-by: Martin Holters --- src/evalpoly.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/evalpoly.jl b/src/evalpoly.jl index 9639d605..40d87492 100644 --- a/src/evalpoly.jl +++ b/src/evalpoly.jl @@ -5,7 +5,7 @@ function _evalpoly(X::AbstractMatrix, p) Base.require_one_based_indexing(p) p0 = isempty(p) ? Base.reduce_empty_iter(+, p) : p[end] Xone = one(X) - S = Base.promote_op(*, typeof(Xone), typeof(Xone))(Xone) * p0 + S = convert(Base.promote_op(*, typeof(Xone), typeof(Xone)), Xone) * p0 for i = length(p)-1:-1:1 S = X * S + @inbounds(p[i] isa AbstractMatrix ? p[i] : p[i] * I) end