diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index f9b562ef..16e12508 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -13,8 +13,7 @@ end function f(g, x, y, n) t = x for i=1:n - t = g(t,x) - t = g(t,y) + t = g(x,y) end t end @@ -22,14 +21,14 @@ function f(g, x, n) t = x for i=1:n t = g(t) - t = g(t) end t end n = 100 +dims = (5,20) -S = SUITE["+"] = BenchmarkGroup() -for i in (5,10,20,50) +S = SUITE["arithmetic"] = BenchmarkGroup() +for i in dims r = rand(i) x = Taylor1(r) x2 = Taylor1(r) @@ -40,20 +39,31 @@ for i in (5,10,20,50) S["STaylor1{$i,Float64} + Float64"] = @benchmarkable f(+, $q, $y, $n) S["Taylor1{$i,Float64} - Float64"] = @benchmarkable f(-, $x, $y, $n) S["STaylor1{$i,Float64} - Float64"] = @benchmarkable f(-, $q, $y, $n) + S["Taylor1{$i,Float64} / Float64"] = @benchmarkable f(/, $x, $y, $n) + S["STaylor1{$i,Float64} / Float64"] = @benchmarkable f(/, $q, $y, $n) + S["Taylor1{$i,Float64} * Float64"] = @benchmarkable f(*, $x, $y, $n) + S["STaylor1{$i,Float64} * Float64"] = @benchmarkable f(*, $q, $y, $n) S["Taylor1{$i,Float64} + Taylor1{$i,Float64}"] = @benchmarkable f(+, $x, $x2, $n) S["STaylor1{$i,Float64} + STaylor1{$i,Float64}"] = @benchmarkable f(+, $q, $q2, $n) S["Taylor1{$i,Float64} - Taylor1{$i,Float64}"] = @benchmarkable f(-, $x, $x2, $n) S["STaylor1{$i,Float64} - STaylor1{$i,Float64}"] = @benchmarkable f(-, $q, $q2, $n) + end S = SUITE["functions"] = BenchmarkGroup() -for i in (5,10,20,50) +for i in dims r = rand(i) x = Taylor1(r) q = STaylor1(r) y = rand() - S["Taylor1{$i,Float64}"] = @benchmarkable f(exp, $x, $n) - S["STaylor1{$i,Float64}"] = @benchmarkable f(exp, $q, $n) + S["exp(Taylor1{Float64}), len = $i"] = @benchmarkable f(exp, $x, $n) + S["exp(STaylor1{$i,Float64}),"] = @benchmarkable f(exp, $q, $n) + S["zero(Taylor1{Float64}), len = $i"] = @benchmarkable f(zero, $x, $n) + S["zero(STaylor1{$i,Float64}),"] = @benchmarkable f(zero, $q, $n) + S["one(Taylor1{Float64}), len = $i"] = @benchmarkable f(one, $x, $n) + S["one(STaylor1{$i,Float64}),"] = @benchmarkable f(one, $q, $n) + S["iszero(Taylor1{Float64}), len = $i"] = @benchmarkable f(iszero, $x, $n) + S["iszero(STaylor1{$i,Float64}),"] = @benchmarkable f(iszero, $q, $n) end #= diff --git a/src/arithmetic.jl b/src/arithmetic.jl index cec59002..b4a2da50 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -305,6 +305,14 @@ end end end + +function *(a::STaylor1{N,T}, b::T) where {N, T<:Number} + STaylor1{N,T}(scale_tuple(a.coeffs, b)) +end +function *(b::T, a::STaylor1{N,T}) where {N, T<:Number} + STaylor1{N,T}(scale_tuple(a.coeffs, b)) +end + for T in (:HomogeneousPolynomial, :TaylorN) @eval begin @@ -508,6 +516,10 @@ function /(a::Taylor1{T}, b::Taylor1{T}) where {T<:Number} return c end +function /(a::STaylor1{N,T}, b::T) where {N, T<:Number} + STaylor1{N,T}(div_tuple_by_scalar(a.coeffs, b)) +end + /(a::TaylorN{T}, b::TaylorN{S}) where {T<:NumberNotSeriesN, S<:NumberNotSeriesN} = /(promote(a,b)...) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 470726d5..946979c5 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -286,6 +286,18 @@ function Base.findlast(a::Taylor1{T}) where {T<:Number} return last-1 end +function Base.findfirst(a::STaylor1{N,T}) where {N, T<:Number} + first = findfirst(x->!iszero(x), a.coeffs) + isa(first, Nothing) && return -1 + return first-1 +end +# Finds the last non-zero entry; extended to Taylor1 +function Base.findlast(a::STaylor1{N,T}) where {N, T<:Number} + last = findlast(x->!iszero(x), a.coeffs) + isa(last, Nothing) && return -1 + return last-1 +end + ## copyto! ## # Inspired from base/abstractarray.jl, line 665 diff --git a/src/broadcasting.jl b/src/broadcasting.jl index 4c7790ea..2679010e 100644 --- a/src/broadcasting.jl +++ b/src/broadcasting.jl @@ -12,7 +12,7 @@ import .Broadcast: BroadcastStyle, Broadcasted, broadcasted # BroadcastStyle definitions and basic precedence rules struct Taylor1Style{T} <: Base.Broadcast.AbstractArrayStyle{0} end -Taylor1Style{T}(::Val{N}) where {T, N}= Base.Broadcast.DefaultArrayStyle{N}() +Taylor1Style{T}(::Val{N}) where {T, N} = Base.Broadcast.DefaultArrayStyle{N}() BroadcastStyle(::Type{<:Taylor1{T}}) where {T} = Taylor1Style{T}() BroadcastStyle(::Taylor1Style{T}, ::Base.Broadcast.DefaultArrayStyle{0}) where {T} = Taylor1Style{T}() BroadcastStyle(::Taylor1Style{T}, ::Base.Broadcast.DefaultArrayStyle{1}) where {T} = Base.Broadcast.DefaultArrayStyle{1}() diff --git a/src/constructors.jl b/src/constructors.jl index 22ba4054..9c7cef3f 100644 --- a/src/constructors.jl +++ b/src/constructors.jl @@ -72,12 +72,31 @@ Taylor1(::Type{T}, order::Int) where {T<:Number} = Taylor1( [zero(T), one(T)], o Taylor1(order::Int) = Taylor1(Float64, order) ######################### STaylor1 +""" + Taylor1{N,T<:Number} <: AbstractSeries{T} + +DataType for polynomial expansions in one independent variable. + +**Fields:** + +- `coeffs :: NTuple{N,T}` Expansion coefficients; the ``i``-th + component is the coefficient of degree ``i-1`` of the expansion. + +Note that `STaylor1` variables are callable. For more information, see +[`evaluate`](@ref). +""" struct STaylor1{N,T<:Number} <: AbstractSeries{T} coeffs::NTuple{N,T} end ## Outer constructors ## -@inline STaylor1(x::STaylor1{N,T}) where {N,T<:Number} = x + +""" + STaylor1(x::T, v::Val{N}) + +Shortcut to define the independent variable of a `STaylor1{N,T}` polynomial of +given `N` with constant term equal to `x`. +""" @generated function STaylor1(x::T, v::Val{N}) where {N,T<:Number} y = Any[zero(T) for i=1:N] tup = :((x,)) @@ -90,6 +109,7 @@ end @inline function STaylor1(coeffs::Vector{T}) where {T<:Number} STaylor1{length(coeffs),T}(tuple(coeffs...)) end +@inline STaylor1(x::STaylor1{N,T}) where {N,T<:Number} = x ######################### HomogeneousPolynomial """ diff --git a/test/onevariable.jl b/test/onevariable.jl index 148c6178..8ffe2874 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -581,6 +581,10 @@ end t1 = STaylor1([1.1, 2.1, 3.1]) t2 = Taylor1([1.1, 2.1, 3.1]) @test test_vs_Taylor1(exp(t1), exp(t2)) + + a = STaylor1([0.0, 1.2, 2.3, 4.5, 0.0]) + @test findfirst(a) == 1 + @test findlast(a) == 3 end #=