From b05db52443f5a51812f73c0c70e2a704920938e6 Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Fri, 3 Apr 2020 18:27:35 -0400 Subject: [PATCH 01/22] add STaylor1 definition, a few arithmetic operators, and benchmarking --- benchmark/REQUIRE | 2 ++ benchmark/benchmarks.jl | 28 +++++++++++++++++++++++++ benchmark/run_benchmark.jl | 12 +++++++++++ benchmark/tune.json | 1 + src/TaylorSeries.jl | 2 +- src/arithmetic.jl | 22 +++++++++++++++++++ src/auxiliary.jl | 43 ++++++++++++++++++++++++++++++++++++++ src/constructors.jl | 19 +++++++++++++++++ src/printing.jl | 20 ++++++++++++++++++ test/onevariable.jl | 9 ++++++++ 10 files changed, 157 insertions(+), 1 deletion(-) create mode 100644 benchmark/REQUIRE create mode 100644 benchmark/benchmarks.jl create mode 100644 benchmark/run_benchmark.jl create mode 100644 benchmark/tune.json diff --git a/benchmark/REQUIRE b/benchmark/REQUIRE new file mode 100644 index 00000000..6ce08119 --- /dev/null +++ b/benchmark/REQUIRE @@ -0,0 +1,2 @@ +BenchmarkTools 0.3.1 +PkgBenchmark 0.1.1 diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl new file mode 100644 index 00000000..a3523d3c --- /dev/null +++ b/benchmark/benchmarks.jl @@ -0,0 +1,28 @@ +using BenchmarkTools, TaylorSeries + +const SUITE = BenchmarkGroup() + +S = SUITE["Constructors"] = BenchmarkGroup() +for i in (05,10,50) + S["STaylor1(x::Array[len = $i])"] = @benchmarkable STaylor1(x) setup = (x = rand($i)) + S["STaylor1(x::Float64, Val($i))"] = @benchmarkable STaylor1(x, Val($i)) setup=(x = rand()) +end + +function f(g, x, y, n) + t = x + for i=1:n + t = g(t,x) + t = g(t,y) + end + t +end + +#= +S = SUITE["STaylor1 Arithmetic"] = BenchmarkGroup() +for i in (5,10,5) + S["STaylor1{$i,Float64} Scalar Addition"] + S["STaylor1{$i,Float64} Scalar Multiplication"] + S["STaylor1{$i,Float64} Addition"] + S["STaylor1{$i,Float64} Multiplication"] +end +=# diff --git a/benchmark/run_benchmark.jl b/benchmark/run_benchmark.jl new file mode 100644 index 00000000..98dd433e --- /dev/null +++ b/benchmark/run_benchmark.jl @@ -0,0 +1,12 @@ +using PkgBenchmark +results = benchmarkpkg("TaylorSeries") +show(results) + +#= +# specify tag and uncommit to benchmark versus prior tagged version +tag = +results = judge("TaylorSeries", tag) +show(results) +=# + +export_markdown("..\\results.md", results) diff --git a/benchmark/tune.json b/benchmark/tune.json new file mode 100644 index 00000000..229325de --- /dev/null +++ b/benchmark/tune.json @@ -0,0 +1 @@ +[{"Julia":"1.3.1","BenchmarkTools":"0.4.3"},[["BenchmarkGroup",{"data":{"Constructors":["BenchmarkGroup",{"data":{"STaylor1(x::Float64, Val(5))":["BenchmarkTools.Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":12,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"STaylor1(x::Array[length = 10])":["BenchmarkTools.Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":15,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"STaylor1(x::Float64, Val(10))":["BenchmarkTools.Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":24,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"STaylor1(x::Float64, Val(50))":["BenchmarkTools.Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":22,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"STaylor1(x::Array[length = 5])":["BenchmarkTools.Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":10,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"STaylor1(x::Array[length = 50])":["BenchmarkTools.Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":10,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}]},"tags":[]}]},"tags":[]}]]] \ No newline at end of file diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl index afcdb0ae..9f8f53b5 100644 --- a/src/TaylorSeries.jl +++ b/src/TaylorSeries.jl @@ -42,7 +42,7 @@ import Base: zero, one, zeros, ones, isinf, isnan, iszero, power_by_squaring, rtoldefault, isfinite, isapprox, rad2deg, deg2rad -export Taylor1, TaylorN, HomogeneousPolynomial, AbstractSeries +export Taylor1, TaylorN, STaylor1, HomogeneousPolynomial, AbstractSeries export getcoeff, derivative, integrate, differentiate, evaluate, evaluate!, inverse, set_taylor1_varname, diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 7a3ec714..5b3a980b 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -205,6 +205,10 @@ end -(a::Taylor1{T}, b::TaylorN{S}) where {T<:NumberNotSeries,S<:NumberNotSeries} = -(promote(a,b)...) +@inline +(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(add_tuples(a.coeffs, b.coeffs)) +@inline -(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(sub_tuples(a.coeffs, b.coeffs)) +@inline -(partials::STaylor1) = STaylor1(minus_tuple(partials.values)) + ## Multiplication ## for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @@ -225,6 +229,24 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) end end +@generated function Base.:*(p1::STaylor1{N, T}, p2::STaylor1{N, T}) where {N, T<:Number} + exprs = Any[nothing for i in 1:N] + for i in 0 : N-1 # order is N-1 + for j in 0:N-1 + k = i + j + 1 # setindex does not have offset + if exprs[k] === nothing + exprs[k] = :(p1[$i] * p2[$j]) + else + exprs[k] = :(muladd(p1[$i], p2[$j], $(exprs[k]))) + end + end + end + return quote + Base.@_inline_meta + STaylor1{N,T}(tuple($(exprs...))) + end +end + for T in (:HomogeneousPolynomial, :TaylorN) @eval begin diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 98741c3b..ebb7ce3e 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -317,3 +317,46 @@ linear_polynomial(a::TaylorN) = TaylorN([a[1]]) linear_polynomial(a::Vector{T}) where {T<:Number} = linear_polynomial.(a) linear_polynomial(a::Number) = a + + +#= +Tuple operations taken from ForwardDiff.jl + +Copyright (c) 2015: Jarrett Revels, Theodore Papamarkou, Miles Lubin, and other contributors + +Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: +The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + +=# +function tupexpr(f, N) + ex = Expr(:tuple, [f(i) for i=1:N]...) + return quote + $(Expr(:meta, :inline)) + @inbounds return $ex + end +end + +@generated function scale_tuple(tup::NTuple{N}, x) where N + return tupexpr(i -> :(tup[$i] * x), N) +end + +@generated function div_tuple_by_scalar(tup::NTuple{N}, x) where N + return tupexpr(i -> :(tup[$i] / x), N) +end + +@generated function add_tuples(a::NTuple{N}, b::NTuple{N}) where N + return tupexpr(i -> :(a[$i] + b[$i]), N) +end + +@generated function sub_tuples(a::NTuple{N}, b::NTuple{N}) where N + return tupexpr(i -> :(a[$i] - b[$i]), N) +end + +@generated function minus_tuple(tup::NTuple{N}) where N + return tupexpr(i -> :(-tup[$i]), N) +end + +@generated function mul_tuples(a::NTuple{N}, b::NTuple{N}, afactor, bfactor) where N + return tupexpr(i -> :((afactor * a[$i]) + (bfactor * b[$i])), N) +end diff --git a/src/constructors.jl b/src/constructors.jl index cbac7c17..22ba4054 100644 --- a/src/constructors.jl +++ b/src/constructors.jl @@ -71,6 +71,25 @@ julia> Taylor1(Rational{Int}, 4) Taylor1(::Type{T}, order::Int) where {T<:Number} = Taylor1( [zero(T), one(T)], order) Taylor1(order::Int) = Taylor1(Float64, order) +######################### STaylor1 +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 +@generated function STaylor1(x::T, v::Val{N}) where {N,T<:Number} + y = Any[zero(T) for i=1:N] + tup = :((x,)) + push!(tup.args, y...) + return quote + Base.@_inline_meta + STaylor1{(N+1),T}($tup) + end +end +@inline function STaylor1(coeffs::Vector{T}) where {T<:Number} + STaylor1{length(coeffs),T}(tuple(coeffs...)) +end ######################### HomogeneousPolynomial """ diff --git a/src/printing.jl b/src/printing.jl index 99319b76..23b9104c 100644 --- a/src/printing.jl +++ b/src/printing.jl @@ -229,3 +229,23 @@ function show(io::IO, a::Union{Taylor1, HomogeneousPolynomial, TaylorN}) return print(io, pretty_print(a)) end end + +# printing for static taylor +function coeffstring(t::STaylor1, i, variable=:t) + if i == 1 # order 0 + return string(t.coeffs[i]) + end + if i == 2 # order 1 + return string(t.coeffs[i], variable) + end + return string(t.coeffs[i], variable, "^", i-1) +end + +function print_taylor(io::IO, t::STaylor1, variable=:t) + print(io, "(" * join([coeffstring(t, i, variable) for i in 1:length(t.coeffs)], " + ") * ")") +end + +Base.show(io::IO, t::STaylor1) = print_taylor(io, t) +function Base.show(io::IO, t::STaylor1{N,T}) where {N, T<:STaylor1} + print_taylor(io, t, :t) +end diff --git a/test/onevariable.jl b/test/onevariable.jl index 0a1687f3..2cffa002 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -523,6 +523,15 @@ eeuler = Base.MathConstants.e @test Taylor1{Int}(false) == Taylor1([0]) end +@testset "Tests for STaylor1 expansions" begin + + @test STaylor1 <: AbstractSeries + @test STaylor1{1,Float64} <: AbstractSeries{Float64} + @test STaylor1([1.0, 2.0]) == STaylor1((1.0, 2.0)) + @test STaylor1(STaylor1((1.0, 2.0))) == STaylor1((1.0, 2.0)) + @test STaylor1(1.0, Val(2)) == STaylor1((1.0, 0.0, 0.0)) +end + @testset "Test `inv` for `Matrix{Taylor1{Float64}}``" begin t = Taylor1(5) a = Diagonal(rand(0:10,3)) + rand(3, 3) From 20d0c47ce862c13b481cf00c8204a2d5da9db753 Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Sat, 4 Apr 2020 11:07:33 -0400 Subject: [PATCH 02/22] Update *, add getindex --- src/arithmetic.jl | 38 ++++++++++++++++++++++---------------- src/auxiliary.jl | 8 ++++++++ 2 files changed, 30 insertions(+), 16 deletions(-) diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 5b3a980b..93e0d134 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -229,22 +229,28 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) end end -@generated function Base.:*(p1::STaylor1{N, T}, p2::STaylor1{N, T}) where {N, T<:Number} - exprs = Any[nothing for i in 1:N] - for i in 0 : N-1 # order is N-1 - for j in 0:N-1 - k = i + j + 1 # setindex does not have offset - if exprs[k] === nothing - exprs[k] = :(p1[$i] * p2[$j]) - else - exprs[k] = :(muladd(p1[$i], p2[$j], $(exprs[k]))) - end - end - end - return quote - Base.@_inline_meta - STaylor1{N,T}(tuple($(exprs...))) - end +@generated function *(x::STaylor1{N,T}, y::STaylor1{N,T}) where {T<:Number,N} + ex_calc = quote end + append!(ex_calc.args, Any[nothing for i in 1:N]) + syms = Symbol[Symbol("c$i") for i in 1:N] + for j = 1:N + ex_line = :(x.coeffs[1]*y.coeffs[$j]) + for k = 2:j + ex_line = :($ex_line + x.coeffs[$k]*y.coeffs[$(j-k+1)]) + end + sym = syms[j] + ex_line = :($sym = $ex_line) + ex_calc.args[j] = ex_line + end + exout = :(($(syms[1]),)) + for i = 2:N + push!(exout.args, syms[i]) + end + return quote + Base.@_inline_meta + $ex_calc + return STaylor1{N,T}($exout) + end end for T in (:HomogeneousPolynomial, :TaylorN) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index ebb7ce3e..9468d8ff 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -104,6 +104,12 @@ function setindex!(a::Taylor1{T}, x::Array{T,1}, u::StepRange{Int,Int}) where {T end +# getindex STaylor1 +getindex(a::STaylor1, n::Int) = a.coeffs[n+1] +getindex(a::STaylor1, u::UnitRange{Int}) = a.coeffs[u .+ 1] +getindex(a::STaylor1, c::Colon) = a.coeffs[c] +getindex(a::STaylor1, u::StepRange{Int,Int}) = a.coeffs[u .+ 1] + """ getcoeff(a, v) @@ -297,6 +303,8 @@ and `TaylorN`. The fallback behavior is to return `a` itself if """ constant_term(a::Taylor1) = a[0] +constant_term(a::STaylor1) = a[0] + constant_term(a::TaylorN) = a[0][1] constant_term(a::Vector{T}) where {T<:Number} = constant_term.(a) From fcbe1897f799230ea0ff1f8c5ad55748ed8507eb Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Sat, 4 Apr 2020 11:20:58 -0400 Subject: [PATCH 03/22] Add size, length, iteration functions --- src/auxiliary.jl | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 9468d8ff..470726d5 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -242,6 +242,15 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @inline axes(a::$T) = () end end +@inline iterate(a::STaylor1{N,T}, state=0) where {N, T<:Number} = state > N-1 ? nothing : (a.coeffs[state+1], state+1) +@inline firstindex(a::STaylor1) = 0 +@inline lastindex(a::STaylor1{N,T}) where {N, T<:Number} = N-1 +@inline eachindex(s::STaylor1{N,T}) where {N, T<:Number} = UnitRange(0, N-1) +@inline size(s::STaylor1{N,T}) where {N, T<:Number} = N +@inline length(s::STaylor1{N,T}) where {N, T<:Number} = N +@inline get_order(s::STaylor1{N,T}) where {N, T<:Number} = N - 1 +@inline eltype(s::STaylor1{N,T}) where {N, T<:Number} = T +@inline axes(a::STaylor1) = () ## fixorder ## From 63195c93e42920b12ae8327fd874fb1a17b11159 Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Sat, 4 Apr 2020 12:10:06 -0400 Subject: [PATCH 04/22] STaylor1: add promotions, conversions, real, img, conj, isinf, isnan --- src/arithmetic.jl | 2 +- src/conversion.jl | 23 +++++++++++++++++++++++ src/other_functions.jl | 6 ++++++ 3 files changed, 30 insertions(+), 1 deletion(-) diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 93e0d134..92f48632 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -207,7 +207,7 @@ end @inline +(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(add_tuples(a.coeffs, b.coeffs)) @inline -(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(sub_tuples(a.coeffs, b.coeffs)) -@inline -(partials::STaylor1) = STaylor1(minus_tuple(partials.values)) +@inline -(a::STaylor1) = STaylor1(minus_tuple(a.coeffs)) ## Multiplication ## diff --git a/src/conversion.jl b/src/conversion.jl index 9b82e5d6..4282dd44 100644 --- a/src/conversion.jl +++ b/src/conversion.jl @@ -28,6 +28,21 @@ convert(::Type{Taylor1{T}}, b::T) where {T<:Number} = Taylor1([b], 0) convert(::Type{Taylor1}, a::T) where {T<:Number} = Taylor1(a, 0) +# Conversion for STaylor1 +function convert(::Type{STaylor1{N,Rational{T}}}, a::STaylor1{N,S}) where {N,T<:Integer, S<:AbstractFloat} + STaylor1{N,T}(rationalize.(a[:])) +end +function convert(::Type{STaylor1{N,T}}, b::Array{T,1}) where {N,T<:Number} + @assert N == length(b) + STaylor1{N,T}(b) +end +function convert(::Type{STaylor1{N,T}}, b::Array{S,1}) where {N,T<:Number, S<:Number} + @assert N == length(b) + STaylor1{N,T}(convert(Array{T,1},b)) +end +convert(::Type{STaylor1{N,T}}, a::STaylor1{N,T}) where {N,T<:Number} = a +convert(::Type{STaylor1{N,T}}, b::S) where {N,T<:Number, S<:Number} = STaylor1(b, Val(N)) +convert(::Type{STaylor1{N,T}}, b::T) where {N, T<:Number} = STaylor1(b, Val(N)) convert(::Type{HomogeneousPolynomial{T}}, a::HomogeneousPolynomial) where {T<:Number} = HomogeneousPolynomial(convert(Array{T,1}, a.coeffs), a.order) @@ -171,6 +186,14 @@ promote_rule(::Type{Taylor1{T}}, ::Type{S}) where {T<:Number, S<:Number} = promote_rule(::Type{Taylor1{Taylor1{T}}}, ::Type{Taylor1{T}}) where {T<:Number} = Taylor1{Taylor1{T}} +promote_rule(::Type{STaylor1{N,T}}, ::Type{STaylor1{N,T}}) where {N, T<:Number} = STaylor1{N,T} +promote_rule(::Type{STaylor1{N,T}}, ::Type{Taylor1{N,S}}) where {N, T<:Number, S<:Number} = STaylor1{N, promote_type(T,S)} +promote_rule(::Type{STaylor1{N,T}}, ::Type{Array{T,1}}) where {N, T<:Number} = STaylor1{N,T} +promote_rule(::Type{STaylor1{N,T}}, ::Type{Array{S,1}}) where {N, T<:Number, S<:Number} = STaylor1{N,promote_type(T,S)} +promote_rule(::Type{STaylor1{N,T}}, ::Type{T}) where {N, T<:Number} = STaylor1{N,T} +promote_rule(::Type{STaylor1{N,T}}, ::Type{S}) where {N, T<:Number, S<:Number} = STaylor1{N,promote_type(T,S)} +promote_rule(::Type{STaylor1{N,T}{STaylor1{Q,R}}}, ::Type{STaylor1{Q,R}}) where {N, Q, R<:Number, T<:Number} = STaylor1{N,T}{STaylor1{Q,R}} + promote_rule(::Type{HomogeneousPolynomial{T}}, ::Type{HomogeneousPolynomial{S}}) where {T<:Number, S<:Number} = diff --git a/src/other_functions.jl b/src/other_functions.jl index e34acbb9..bacdfdc3 100644 --- a/src/other_functions.jl +++ b/src/other_functions.jl @@ -20,7 +20,13 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @eval isnan(a::$T) = any( isnan.(a.coeffs) ) end +for f in (:real, :imag, :conj) + @eval ($f)(a::STaylor1{N,T}) where {N,T<:Number} = STaylor1(($f).(a.coeffs), Val(N)) +end +adjoint(a::STaylor1) = conj(a) +isinf(a::STaylor1) = any(isinf.(a.coeffs)) +isnan(a::STaylor1) = any(isnan.(a.coeffs)) ## Division functions: rem and mod ## for op in (:mod, :rem) From 51e80c1446b8765ab14f7454e6ccda9ea3baa7cc Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Sat, 4 Apr 2020 11:20:58 -0400 Subject: [PATCH 05/22] STaylor1: Add size, length, iteration functions --- src/auxiliary.jl | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 9468d8ff..470726d5 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -242,6 +242,15 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @inline axes(a::$T) = () end end +@inline iterate(a::STaylor1{N,T}, state=0) where {N, T<:Number} = state > N-1 ? nothing : (a.coeffs[state+1], state+1) +@inline firstindex(a::STaylor1) = 0 +@inline lastindex(a::STaylor1{N,T}) where {N, T<:Number} = N-1 +@inline eachindex(s::STaylor1{N,T}) where {N, T<:Number} = UnitRange(0, N-1) +@inline size(s::STaylor1{N,T}) where {N, T<:Number} = N +@inline length(s::STaylor1{N,T}) where {N, T<:Number} = N +@inline get_order(s::STaylor1{N,T}) where {N, T<:Number} = N - 1 +@inline eltype(s::STaylor1{N,T}) where {N, T<:Number} = T +@inline axes(a::STaylor1) = () ## fixorder ## From b4eafcab38798bce0bfd14843878e2e62f84b25c Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Sat, 4 Apr 2020 12:10:06 -0400 Subject: [PATCH 06/22] STaylor1: add promotions, conversions, real, img, conj, isinf, isnan --- src/arithmetic.jl | 2 +- src/conversion.jl | 23 +++++++++++++++++++++++ src/other_functions.jl | 6 ++++++ 3 files changed, 30 insertions(+), 1 deletion(-) diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 93e0d134..92f48632 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -207,7 +207,7 @@ end @inline +(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(add_tuples(a.coeffs, b.coeffs)) @inline -(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(sub_tuples(a.coeffs, b.coeffs)) -@inline -(partials::STaylor1) = STaylor1(minus_tuple(partials.values)) +@inline -(a::STaylor1) = STaylor1(minus_tuple(a.coeffs)) ## Multiplication ## diff --git a/src/conversion.jl b/src/conversion.jl index 9b82e5d6..4282dd44 100644 --- a/src/conversion.jl +++ b/src/conversion.jl @@ -28,6 +28,21 @@ convert(::Type{Taylor1{T}}, b::T) where {T<:Number} = Taylor1([b], 0) convert(::Type{Taylor1}, a::T) where {T<:Number} = Taylor1(a, 0) +# Conversion for STaylor1 +function convert(::Type{STaylor1{N,Rational{T}}}, a::STaylor1{N,S}) where {N,T<:Integer, S<:AbstractFloat} + STaylor1{N,T}(rationalize.(a[:])) +end +function convert(::Type{STaylor1{N,T}}, b::Array{T,1}) where {N,T<:Number} + @assert N == length(b) + STaylor1{N,T}(b) +end +function convert(::Type{STaylor1{N,T}}, b::Array{S,1}) where {N,T<:Number, S<:Number} + @assert N == length(b) + STaylor1{N,T}(convert(Array{T,1},b)) +end +convert(::Type{STaylor1{N,T}}, a::STaylor1{N,T}) where {N,T<:Number} = a +convert(::Type{STaylor1{N,T}}, b::S) where {N,T<:Number, S<:Number} = STaylor1(b, Val(N)) +convert(::Type{STaylor1{N,T}}, b::T) where {N, T<:Number} = STaylor1(b, Val(N)) convert(::Type{HomogeneousPolynomial{T}}, a::HomogeneousPolynomial) where {T<:Number} = HomogeneousPolynomial(convert(Array{T,1}, a.coeffs), a.order) @@ -171,6 +186,14 @@ promote_rule(::Type{Taylor1{T}}, ::Type{S}) where {T<:Number, S<:Number} = promote_rule(::Type{Taylor1{Taylor1{T}}}, ::Type{Taylor1{T}}) where {T<:Number} = Taylor1{Taylor1{T}} +promote_rule(::Type{STaylor1{N,T}}, ::Type{STaylor1{N,T}}) where {N, T<:Number} = STaylor1{N,T} +promote_rule(::Type{STaylor1{N,T}}, ::Type{Taylor1{N,S}}) where {N, T<:Number, S<:Number} = STaylor1{N, promote_type(T,S)} +promote_rule(::Type{STaylor1{N,T}}, ::Type{Array{T,1}}) where {N, T<:Number} = STaylor1{N,T} +promote_rule(::Type{STaylor1{N,T}}, ::Type{Array{S,1}}) where {N, T<:Number, S<:Number} = STaylor1{N,promote_type(T,S)} +promote_rule(::Type{STaylor1{N,T}}, ::Type{T}) where {N, T<:Number} = STaylor1{N,T} +promote_rule(::Type{STaylor1{N,T}}, ::Type{S}) where {N, T<:Number, S<:Number} = STaylor1{N,promote_type(T,S)} +promote_rule(::Type{STaylor1{N,T}{STaylor1{Q,R}}}, ::Type{STaylor1{Q,R}}) where {N, Q, R<:Number, T<:Number} = STaylor1{N,T}{STaylor1{Q,R}} + promote_rule(::Type{HomogeneousPolynomial{T}}, ::Type{HomogeneousPolynomial{S}}) where {T<:Number, S<:Number} = diff --git a/src/other_functions.jl b/src/other_functions.jl index e34acbb9..bacdfdc3 100644 --- a/src/other_functions.jl +++ b/src/other_functions.jl @@ -20,7 +20,13 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @eval isnan(a::$T) = any( isnan.(a.coeffs) ) end +for f in (:real, :imag, :conj) + @eval ($f)(a::STaylor1{N,T}) where {N,T<:Number} = STaylor1(($f).(a.coeffs), Val(N)) +end +adjoint(a::STaylor1) = conj(a) +isinf(a::STaylor1) = any(isinf.(a.coeffs)) +isnan(a::STaylor1) = any(isnan.(a.coeffs)) ## Division functions: rem and mod ## for op in (:mod, :rem) From 6ca0206c7953abd5eba6b0321f2ea194bafc1427 Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Sat, 4 Apr 2020 14:28:12 -0400 Subject: [PATCH 07/22] STaylor1: add iszero, zero, one, == --- src/arithmetic.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 92f48632..1748e0e7 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -23,6 +23,9 @@ for T in (:Taylor1, :TaylorN) end end +==(a::STaylor1{N,T}, b::STaylor1{N,S}) where {N, T<:Number, S<:Number} = ==(promote(a,b)...) +==(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = (a.coeffs == b.coeffs) + function ==(a::HomogeneousPolynomial, b::HomogeneousPolynomial) a.order == b.order && return a.coeffs == b.coeffs return iszero(a.coeffs) && iszero(b.coeffs) @@ -31,11 +34,14 @@ end for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @eval iszero(a::$T) = iszero(a.coeffs) end +iszero(a::STaylor1) = iszero(a.coeffs) ## zero and one ## for T in (:Taylor1, :TaylorN), f in (:zero, :one) @eval ($f)(a::$T) = $T(($f)(a[0]), a.order) end +zero(a::STaylor1{N,T}) where {N, T<:Number} = STaylor1(zero(T),Val(N)) +one(a::STaylor1{N,T}) where {N, T<:Number} = STaylor1(one(T),Val(N)) function zero(a::HomogeneousPolynomial{T}) where {T<:Number} v = zero.(a.coeffs) @@ -207,6 +213,7 @@ end @inline +(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(add_tuples(a.coeffs, b.coeffs)) @inline -(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(sub_tuples(a.coeffs, b.coeffs)) +@inline +(a::STaylor1) = a @inline -(a::STaylor1) = STaylor1(minus_tuple(a.coeffs)) From d4a5117b99e5568daf8ace8d2d466d5d46e1903f Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Sat, 4 Apr 2020 15:24:21 -0400 Subject: [PATCH 08/22] STaylor1: +(STaylor,Number), -(STaylor,Number), ... and additional tests --- src/arithmetic.jl | 37 +++++++++++++++++++++++++++++++++++++ src/conversion.jl | 4 +--- test/onevariable.jl | 24 ++++++++++++++++++++++++ 3 files changed, 62 insertions(+), 3 deletions(-) diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 1748e0e7..afbea5f0 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -216,6 +216,43 @@ end @inline +(a::STaylor1) = a @inline -(a::STaylor1) = STaylor1(minus_tuple(a.coeffs)) +@generated function add_first_tuple(tup::NTuple{N,T}, x::T) where {N, T<:Number} + ex = :((tup[1] + x),) + append!(ex.args, Any[:(tup[$i]) for i in 2:N]) + return quote + Base.@_inline_meta + @inbounds return $ex + end +end + +@generated function sub_first_tuple(tup::NTuple{N,T}, x::T) where {N, T<:Number} + ex = :((tup[1] - x),) + append!(ex.args, Any[:(tup[$i]) for i in 2:N]) + return quote + Base.@_inline_meta + @inbounds return $ex + end +end + +function +(a::STaylor1{N,T}, b::T) where {N, T<:Number} + STaylor1{N,T}(add_first_tuple(a.coeffs, b)) +end +function +(b::T, a::STaylor1{N,T}) where {N, T<:Number} + STaylor1{N,T}(add_first_tuple(a.coeffs, b)) +end +function -(a::STaylor1{N,T}, b::T) where {N, T<:Number} + STaylor1{N,T}(sub_first_tuple(a.coeffs, b)) +end +-(b::T, a::STaylor1{N,T}) where {N, T<:Number} = b + (-a) + ++(a::STaylor1{N,T}, b::STaylor1{N,S}) where {N, T<:Number, S<:Number} = +(promote(a,b)...) ++(a::STaylor1{N,T}, b::S) where {N, T<:Number, S<:Number} = +(promote(a,b)...) ++(b::S, a::STaylor1{N,T}) where {N, T<:Number, S<:Number} = +(promote(b,a)...) + +-(a::STaylor1{N,T}, b::STaylor1{N,S}) where {N, T<:Number, S<:Number} = -(promote(a,b)...) +-(a::STaylor1{N,T}, b::S) where {N, T<:Number, S<:Number} = -(promote(a,b)...) +-(b::S, a::STaylor1{N,T}) where {N, T<:Number, S<:Number} = -(promote(b,a)...) + ## Multiplication ## for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) diff --git a/src/conversion.jl b/src/conversion.jl index 4282dd44..d16f653e 100644 --- a/src/conversion.jl +++ b/src/conversion.jl @@ -187,13 +187,11 @@ promote_rule(::Type{Taylor1{Taylor1{T}}}, ::Type{Taylor1{T}}) where {T<:Number} Taylor1{Taylor1{T}} promote_rule(::Type{STaylor1{N,T}}, ::Type{STaylor1{N,T}}) where {N, T<:Number} = STaylor1{N,T} -promote_rule(::Type{STaylor1{N,T}}, ::Type{Taylor1{N,S}}) where {N, T<:Number, S<:Number} = STaylor1{N, promote_type(T,S)} +promote_rule(::Type{STaylor1{N,T}}, ::Type{STaylor1{N,S}}) where {N, T<:Number, S<:Number} = STaylor1{N, promote_type(T,S)} promote_rule(::Type{STaylor1{N,T}}, ::Type{Array{T,1}}) where {N, T<:Number} = STaylor1{N,T} promote_rule(::Type{STaylor1{N,T}}, ::Type{Array{S,1}}) where {N, T<:Number, S<:Number} = STaylor1{N,promote_type(T,S)} promote_rule(::Type{STaylor1{N,T}}, ::Type{T}) where {N, T<:Number} = STaylor1{N,T} promote_rule(::Type{STaylor1{N,T}}, ::Type{S}) where {N, T<:Number, S<:Number} = STaylor1{N,promote_type(T,S)} -promote_rule(::Type{STaylor1{N,T}{STaylor1{Q,R}}}, ::Type{STaylor1{Q,R}}) where {N, Q, R<:Number, T<:Number} = STaylor1{N,T}{STaylor1{Q,R}} - promote_rule(::Type{HomogeneousPolynomial{T}}, ::Type{HomogeneousPolynomial{S}}) where {T<:Number, S<:Number} = diff --git a/test/onevariable.jl b/test/onevariable.jl index 2cffa002..e31c6e15 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -530,6 +530,30 @@ end @test STaylor1([1.0, 2.0]) == STaylor1((1.0, 2.0)) @test STaylor1(STaylor1((1.0, 2.0))) == STaylor1((1.0, 2.0)) @test STaylor1(1.0, Val(2)) == STaylor1((1.0, 0.0, 0.0)) + + @test +STaylor1([1.0, 2.0, 3.0]) == STaylor1([1.0, 2.0, 3.0]) + @test -STaylor1([1.0, 2.0, 3.0]) == -STaylor1([1.0, 2.0, 3.0]) + @test STaylor1([1.0, 2.0, 3.0]) + STaylor1([3.0, 2.0, 3.0]) == STaylor1([4.0, 4.0, 6.0]) + @test STaylor1([1.0, 2.0, 3.0]) - STaylor1([3.0, 2.0, 4.0]) == STaylor1([-2.0, 0.0, -1.0]) + @test STaylor1([1.0, 2.0, 3.0]) + 2.0 == STaylor1([3.0, 2.0, 3.0]) + @test STaylor1([1.0, 2.0, 3.0]) - 2.0 == STaylor1([-1.0, 2.0, 3.0]) + @test 2.0 + STaylor1([1.0, 2.0, 3.0]) == STaylor1([3.0, 2.0, 3.0]) + @test 2.0 - STaylor1([1.0, 2.0, 3.0]) == STaylor1([1.0, -2.0, -3.0]) + @test 2 - STaylor1([1.0, 2.0, 3.0]) == STaylor1([1.0, -2.0, -3.0]) + @test STaylor1([1.0, 2.0, 3.0]) - 2 == STaylor1([-1.0, 2.0, 3.0]) + @test STaylor1([1.0, 2.0, 3.0]) + 2 == STaylor1([3.0, 2.0, 3.0]) + @test 2 + STaylor1([1.0, 2.0, 3.0]) == STaylor1([3.0, 2.0, 3.0]) + + @test zero(STaylor1([1.0, 2.0, 3.0])) == STaylor1([0.0, 0.0, 0.0]) + @test one(STaylor1([1.0, 2.0, 3.0])) == STaylor1([1.0, 0.0, 0.0]) + @test isinf(STaylor1([Inf, 2.0, 3.0])) && ~isinf(STaylor1([0.0, 0.0, 0.0])) + @test isnan(STaylor1([NaN, 2.0, 3.0])) && ~isnan(STaylor1([1.0, 0.0, 0.0])) + @test iszero(STaylor1([0.0, 0.0, 0.0])) && ~iszero(STaylor1([0.0, 1.0, 0.0])) + + @test length(STaylor1([0.0, 0.0, 0.0])) == 3 + @test size(STaylor1([0.0, 0.0, 0.0])) == 3 + @test firstindex(STaylor1([0.0, 0.0, 0.0])) == 0 + @test lastindex(STaylor1([0.0, 0.0, 0.0])) == 2 end @testset "Test `inv` for `Matrix{Taylor1{Float64}}``" begin From 79104621b8511bbae9060fe2761695038f578aaf Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Sun, 5 Apr 2020 23:01:37 -0400 Subject: [PATCH 09/22] STaylor1: add unit tests, exp, temp fix promotions --- benchmark/benchmarks.jl | 45 +++++++++++++++++++++++++++++++++++--- benchmark/run_benchmark.jl | 2 +- src/arithmetic.jl | 12 +++++----- src/functions.jl | 33 +++++++++++++++++++++++++++- test/onevariable.jl | 20 +++++++++++++---- 5 files changed, 97 insertions(+), 15 deletions(-) diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index a3523d3c..f9b562ef 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -3,9 +3,11 @@ using BenchmarkTools, TaylorSeries const SUITE = BenchmarkGroup() S = SUITE["Constructors"] = BenchmarkGroup() -for i in (05,10,50) - S["STaylor1(x::Array[len = $i])"] = @benchmarkable STaylor1(x) setup = (x = rand($i)) - S["STaylor1(x::Float64, Val($i))"] = @benchmarkable STaylor1(x, Val($i)) setup=(x = rand()) +for i in (05,10,20,50) + x = rand(i) + t = tuple(x...) + S["STaylor1(x::Array[len = $i])"] = @benchmarkable STaylor1($x) + S["STaylor1(x::Float64, Val($i))"] = @benchmarkable STaylor1($t) end function f(g, x, y, n) @@ -16,6 +18,43 @@ function f(g, x, y, n) end t end +function f(g, x, n) + t = x + for i=1:n + t = g(t) + t = g(t) + end + t +end +n = 100 + +S = SUITE["+"] = BenchmarkGroup() +for i in (5,10,20,50) + r = rand(i) + x = Taylor1(r) + x2 = Taylor1(r) + q = STaylor1(r) + q2 = STaylor1(r) + y = rand() + 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) + 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) +end #= S = SUITE["STaylor1 Arithmetic"] = BenchmarkGroup() diff --git a/benchmark/run_benchmark.jl b/benchmark/run_benchmark.jl index 98dd433e..c074a539 100644 --- a/benchmark/run_benchmark.jl +++ b/benchmark/run_benchmark.jl @@ -9,4 +9,4 @@ results = judge("TaylorSeries", tag) show(results) =# -export_markdown("..\\results.md", results) +export_markdown("results.md", results) diff --git a/src/arithmetic.jl b/src/arithmetic.jl index afbea5f0..8ee2e3f0 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -245,13 +245,13 @@ function -(a::STaylor1{N,T}, b::T) where {N, T<:Number} end -(b::T, a::STaylor1{N,T}) where {N, T<:Number} = b + (-a) -+(a::STaylor1{N,T}, b::STaylor1{N,S}) where {N, T<:Number, S<:Number} = +(promote(a,b)...) -+(a::STaylor1{N,T}, b::S) where {N, T<:Number, S<:Number} = +(promote(a,b)...) -+(b::S, a::STaylor1{N,T}) where {N, T<:Number, S<:Number} = +(promote(b,a)...) +#+(a::STaylor1{N,T}, b::STaylor1{N,S}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(a,b)...) +#+(a::STaylor1{N,T}, b::S) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(a,b)...) +#+(b::S, a::STaylor1{N,T}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(b,a)...) --(a::STaylor1{N,T}, b::STaylor1{N,S}) where {N, T<:Number, S<:Number} = -(promote(a,b)...) --(a::STaylor1{N,T}, b::S) where {N, T<:Number, S<:Number} = -(promote(a,b)...) --(b::S, a::STaylor1{N,T}) where {N, T<:Number, S<:Number} = -(promote(b,a)...) +#-(a::STaylor1{N,T}, b::STaylor1{N,S}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = -(promote(a,b)...) +#-(a::STaylor1{N,T}, b::S) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = -(promote(a,b)...) +#-(b::S, a::STaylor1{N,T}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = -(promote(b,a)...) ## Multiplication ## diff --git a/src/functions.jl b/src/functions.jl index f642e374..28f623b0 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -147,6 +147,38 @@ for T in (:Taylor1, :TaylorN) end end +# Functions for StaticTaylor +@generated function exp(a::STaylor1{N,T}) where {N, T <: Number} + ex_calc = quote end + append!(ex_calc.args, Any[nothing for i in 1:N]) + syms = Symbol[Symbol("c$i") for i in 1:N] + + sym = syms[1] + ex_line = :($(syms[1]) = exp(a[0])) + ex_calc.args[1] = ex_line + + for k in 1:(N-1) + kT = convert(T,k) + sym = syms[k+1] + ex_line = :($kT * a[$k] * $(syms[1])) + @inbounds for i = 1:k-1 + ex_line = :($ex_line + $(kT-i) * a[$(k-i)] * $(syms[i+1])) + end + ex_line = :(($ex_line)/$kT) + ex_line = :($sym = $ex_line) + ex_calc.args[k+1] = ex_line + end + + exout = :(($(syms[1]),)) + for i = 2:N + push!(exout.args, syms[i]) + end + return quote + Base.@_inline_meta + $ex_calc + return STaylor1{N,T}($exout) + end +end # Recursive functions (homogeneous coefficients) for T in (:Taylor1, :TaylorN) @@ -422,7 +454,6 @@ for T in (:Taylor1, :TaylorN) end end - @doc doc""" inverse(f) diff --git a/test/onevariable.jl b/test/onevariable.jl index e31c6e15..80ffeab6 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -523,6 +523,17 @@ eeuler = Base.MathConstants.e @test Taylor1{Int}(false) == Taylor1([0]) end +function test_vs_Taylor1(x,y) + flag = true + for i in 0:2 + if x[i] !== y[i] + flag = false + break + end + end + flag +end + @testset "Tests for STaylor1 expansions" begin @test STaylor1 <: AbstractSeries @@ -539,10 +550,6 @@ end @test STaylor1([1.0, 2.0, 3.0]) - 2.0 == STaylor1([-1.0, 2.0, 3.0]) @test 2.0 + STaylor1([1.0, 2.0, 3.0]) == STaylor1([3.0, 2.0, 3.0]) @test 2.0 - STaylor1([1.0, 2.0, 3.0]) == STaylor1([1.0, -2.0, -3.0]) - @test 2 - STaylor1([1.0, 2.0, 3.0]) == STaylor1([1.0, -2.0, -3.0]) - @test STaylor1([1.0, 2.0, 3.0]) - 2 == STaylor1([-1.0, 2.0, 3.0]) - @test STaylor1([1.0, 2.0, 3.0]) + 2 == STaylor1([3.0, 2.0, 3.0]) - @test 2 + STaylor1([1.0, 2.0, 3.0]) == STaylor1([3.0, 2.0, 3.0]) @test zero(STaylor1([1.0, 2.0, 3.0])) == STaylor1([0.0, 0.0, 0.0]) @test one(STaylor1([1.0, 2.0, 3.0])) == STaylor1([1.0, 0.0, 0.0]) @@ -554,6 +561,11 @@ end @test size(STaylor1([0.0, 0.0, 0.0])) == 3 @test firstindex(STaylor1([0.0, 0.0, 0.0])) == 0 @test lastindex(STaylor1([0.0, 0.0, 0.0])) == 2 + + # check that STaylor1 and Taylor yeild same result + t1 = STaylor1([1.1, 2.1, 3.1]) + t2 = Taylor1([1.1, 2.1, 3.1]) + @test test_vs_Taylor1(exp(t1), exp(t2)) end @testset "Test `inv` for `Matrix{Taylor1{Float64}}``" begin From aa20279adaa31aec83f181b045af851ead486368 Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Mon, 6 Apr 2020 10:54:41 -0400 Subject: [PATCH 10/22] STaylor1: fix zero/one, onevariable tests on STaylor1 only (temporary) --- src/arithmetic.jl | 18 +++++++++++++----- test/onevariable.jl | 6 ++++++ 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 8ee2e3f0..cec59002 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -23,8 +23,7 @@ for T in (:Taylor1, :TaylorN) end end -==(a::STaylor1{N,T}, b::STaylor1{N,S}) where {N, T<:Number, S<:Number} = ==(promote(a,b)...) -==(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = (a.coeffs == b.coeffs) +==(a::STaylor1, b::STaylor1) = (a.coeffs == b.coeffs) function ==(a::HomogeneousPolynomial, b::HomogeneousPolynomial) a.order == b.order && return a.coeffs == b.coeffs @@ -34,14 +33,23 @@ end for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @eval iszero(a::$T) = iszero(a.coeffs) end -iszero(a::STaylor1) = iszero(a.coeffs) +function iszero(a::STaylor1) + flag = true + for i in eachindex(a) + flag &= iszero(a[i]) + if ~flag + break + end + end + flag +end ## zero and one ## for T in (:Taylor1, :TaylorN), f in (:zero, :one) @eval ($f)(a::$T) = $T(($f)(a[0]), a.order) end -zero(a::STaylor1{N,T}) where {N, T<:Number} = STaylor1(zero(T),Val(N)) -one(a::STaylor1{N,T}) where {N, T<:Number} = STaylor1(one(T),Val(N)) +zero(::STaylor1{N,T}) where {N, T<:Number} = STaylor1(zero(T),Val(N-1)) +one(::STaylor1{N,T}) where {N, T<:Number} = STaylor1(one(T),Val(N-1)) function zero(a::HomogeneousPolynomial{T}) where {T<:Number} v = zero.(a.coeffs) diff --git a/test/onevariable.jl b/test/onevariable.jl index 80ffeab6..73d170e1 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -7,6 +7,7 @@ using Test using LinearAlgebra, SparseArrays eeuler = Base.MathConstants.e +#= @testset "Tests for Taylor1 expansions" begin ta(a) = Taylor1([a,one(a)],15) t = Taylor1(Int,15) @@ -522,6 +523,7 @@ eeuler = Base.MathConstants.e @test Taylor1{Int}(true) == Taylor1([1]) @test Taylor1{Int}(false) == Taylor1([0]) end +=# function test_vs_Taylor1(x,y) flag = true @@ -542,6 +544,7 @@ end @test STaylor1(STaylor1((1.0, 2.0))) == STaylor1((1.0, 2.0)) @test STaylor1(1.0, Val(2)) == STaylor1((1.0, 0.0, 0.0)) + @test +STaylor1([1.0, 2.0, 3.0]) == STaylor1([1.0, 2.0, 3.0]) @test -STaylor1([1.0, 2.0, 3.0]) == -STaylor1([1.0, 2.0, 3.0]) @test STaylor1([1.0, 2.0, 3.0]) + STaylor1([3.0, 2.0, 3.0]) == STaylor1([4.0, 4.0, 6.0]) @@ -553,6 +556,7 @@ end @test zero(STaylor1([1.0, 2.0, 3.0])) == STaylor1([0.0, 0.0, 0.0]) @test one(STaylor1([1.0, 2.0, 3.0])) == STaylor1([1.0, 0.0, 0.0]) + @test isinf(STaylor1([Inf, 2.0, 3.0])) && ~isinf(STaylor1([0.0, 0.0, 0.0])) @test isnan(STaylor1([NaN, 2.0, 3.0])) && ~isnan(STaylor1([1.0, 0.0, 0.0])) @test iszero(STaylor1([0.0, 0.0, 0.0])) && ~iszero(STaylor1([0.0, 1.0, 0.0])) @@ -568,6 +572,7 @@ end @test test_vs_Taylor1(exp(t1), exp(t2)) end +#= @testset "Test `inv` for `Matrix{Taylor1{Float64}}``" begin t = Taylor1(5) a = Diagonal(rand(0:10,3)) + rand(3, 3) @@ -637,3 +642,4 @@ end @test_throws DimensionMismatch mul!(Y[1:end-1],A,B) end end +=# From ec4940e20832b193d82cdbca5b183b958f62c99e Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Mon, 6 Apr 2020 14:42:09 -0400 Subject: [PATCH 11/22] STaylor1: Add evaluate and tests --- src/evaluate.jl | 47 +++++++++++++++++++++++++++++++++++++-------- test/onevariable.jl | 11 +++++++++++ 2 files changed, 50 insertions(+), 8 deletions(-) diff --git a/src/evaluate.jl b/src/evaluate.jl index e521be6f..199c126f 100644 --- a/src/evaluate.jl +++ b/src/evaluate.jl @@ -10,9 +10,10 @@ """ evaluate(a, [dx]) -Evaluate a `Taylor1` polynomial using Horner's rule (hand coded). If `dx` is -ommitted, its value is considered as zero. Note that the syntax `a(dx)` is -equivalent to `evaluate(a,dx)`, and `a()` is equivalent to `evaluate(a)`. +Evaluate a `Taylor1` (or `STaylor1{N,T}`) polynomial using Horner's rule (hand +coded). If `dx` is ommitted, its value is considered as zero. Note that the +syntax `a(dx)` is equivalent to `evaluate(a,dx)`, and `a()` is equivalent +to `evaluate(a)`. """ function evaluate(a::Taylor1{T}, dx::T) where {T<:Number} @inbounds suma = a[end] @@ -30,18 +31,41 @@ function evaluate(a::Taylor1{T}, dx::S) where {T<:Number, S<:Number} end evaluate(a::Taylor1{T}) where {T<:Number} = a[0] +function evaluate(a::STaylor1{N,T}, dx::T) where {N, T<:Number} + @inbounds suma = a[N-1] + @inbounds for k in (N-1):-1:0 + suma = suma*dx + a[k] + end + suma +end +#= +function evaluate(a::STaylor1{N,T}, dx::S) where {N, T<:Number, S<:Number} + suma = a[N-1]*one(dx) + @inbounds for k in (N-1):-1:0 + suma = suma*dx + a[k] + end + suma +end +=# +evaluate(a::STaylor1{N,T}) where {N, T<:Number} = a[0] + """ evaluate(x, δt) -Evaluates each element of `x::Union{ Vector{Taylor1{T}}, Matrix{Taylor1{T}} }`, -representing the dependent variables of an ODE, at *time* δt. Note that the -syntax `x(δt)` is equivalent to `evaluate(x, δt)`, and `x()` -is equivalent to `evaluate(x)`. +Evaluates each element of `x::Union{Vector{Taylor1{T}}, Matrix{Taylor1{T}}, +Vector{STaylor1{N,T}}, Matrix{STaylor1{N,T}}}`, representing the dependent +variables of an ODE, at *time* δt. Note that the syntax `x(δt)` is equivalent +to `evaluate(x, δt)`, and `x()` is equivalent to `evaluate(x)`. """ evaluate(x::Union{Array{Taylor1{T}}, SubArray{Taylor1{T}}}, δt::S) where {T<:Number, S<:Number} = evaluate.(x, δt) evaluate(a::Union{Array{Taylor1{T}}, SubArray{Taylor1{T}}}) where {T<:Number} = evaluate.(a, zero(T)) +evaluate(x::Union{Array{STaylor1{N,T}}, SubArray{STaylor1{N,T}}}, δt::S) where + {N, T<:Number, S<:Number} = evaluate.(x, δt) +evaluate(a::Union{Array{STaylor1{N,T}}, SubArray{STaylor1{N,T}}}) where + {N, T<:Number} = evaluate.(a, zero(T)) + """ evaluate!(x, δt, x0) @@ -110,15 +134,22 @@ evaluate(p::Taylor1{T}, x::Array{S}) where {T<:Number, S<:Number} = #function-like behavior for Taylor1 (p::Taylor1)(x) = evaluate(p, x) - (p::Taylor1)() = evaluate(p) +(p::STaylor1)(x) = evaluate(p, x) +(p::STaylor1)() = evaluate(p) + #function-like behavior for Vector{Taylor1} (p::Array{Taylor1{T}})(x) where {T<:Number} = evaluate.(p, x) (p::SubArray{Taylor1{T}})(x) where {T<:Number} = evaluate.(p, x) (p::Array{Taylor1{T}})() where {T<:Number} = evaluate.(p) (p::SubArray{Taylor1{T}})() where {T<:Number} = evaluate.(p) +(p::Array{STaylor1{N,T}})(x) where {N,T<:Number} = evaluate.(p, x) +(p::SubArray{STaylor1{N,T}})(x) where {N,T<:Number} = evaluate.(p, x) +(p::Array{STaylor1{N,T}})() where {N,T<:Number} = evaluate.(p) +(p::SubArray{STaylor1{N,T}})() where {N,T<:Number} = evaluate.(p) + ## Evaluation of multivariable function evaluate!(x::Array{TaylorN{T},1}, δx::Array{T,1}, x0::Array{T,1}) where {T<:Number} diff --git a/test/onevariable.jl b/test/onevariable.jl index 73d170e1..148c6178 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -566,6 +566,17 @@ end @test firstindex(STaylor1([0.0, 0.0, 0.0])) == 0 @test lastindex(STaylor1([0.0, 0.0, 0.0])) == 2 + st1 = STaylor1([1.0, 2.0, 3.0]) + @test st1(2.0) == 41.0 + @test st1() == 1.00 + st2 = typeof(st1)[st1; st1] + @test st2(2.0)[1] == st2(2.0)[2] == 41.0 + @test st2()[1] == st2()[2] == 1.0 + @test evaluate(st1,2.0) == 41.0 + @test evaluate(st1) == 1.00 + @test evaluate(st2,2.0)[1] == evaluate(st2,2.0)[2] == 41.0 + @test evaluate(st2)[1] == evaluate(st2)[2] == 1.0 + # check that STaylor1 and Taylor yeild same result t1 = STaylor1([1.1, 2.1, 3.1]) t2 = Taylor1([1.1, 2.1, 3.1]) From c495b896bf2eb41e0d93d4ec677f2bcbdd546d23 Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Wed, 8 Apr 2020 08:40:11 -0400 Subject: [PATCH 12/22] STaylor1: add findfirst, findlast, *Num, /Num, construtor docs --- benchmark/benchmarks.jl | 26 ++++++++++++++++++-------- src/arithmetic.jl | 12 ++++++++++++ src/auxiliary.jl | 12 ++++++++++++ src/broadcasting.jl | 2 +- src/constructors.jl | 22 +++++++++++++++++++++- test/onevariable.jl | 4 ++++ 6 files changed, 68 insertions(+), 10 deletions(-) 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 #= From 5e051c60c9dfcbcb032686d6cb172e163dbe7a7d Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Wed, 8 Apr 2020 10:45:20 -0400 Subject: [PATCH 13/22] STaylor1: add division, abs, rad2deg, deg2rad; fix real: imag, conj --- benchmark/benchmarks.jl | 28 +++++++++------------------- src/arithmetic.jl | 39 +++++++++++++++++++++++++++++++++++++++ src/other_functions.jl | 20 +++++++++++++++++++- test/onevariable.jl | 25 ++++++++++++++++++++++++- 4 files changed, 91 insertions(+), 21 deletions(-) diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index 16e12508..c765e4be 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -47,7 +47,10 @@ for i in dims 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) - + 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() @@ -56,22 +59,9 @@ for i in dims x = Taylor1(r) q = STaylor1(r) y = rand() - 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 - -#= -S = SUITE["STaylor1 Arithmetic"] = BenchmarkGroup() -for i in (5,10,5) - S["STaylor1{$i,Float64} Scalar Addition"] - S["STaylor1{$i,Float64} Scalar Multiplication"] - S["STaylor1{$i,Float64} Addition"] - S["STaylor1{$i,Float64} Multiplication"] + for g in (exp, abs, zero, one, real, imag, conj, adjoint, iszero, isnan, + isinf, deg2rad, rad2deg) + S["$(Symbol(g))(Taylor1{Float64}), len = $i"] = @benchmarkable f($g, $x, $n) + S["$(Symbol(g))(STaylor1{$i,Float64}),"] = @benchmarkable f($g, $q, $n) + end end -=# diff --git a/src/arithmetic.jl b/src/arithmetic.jl index b4a2da50..e4301be1 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -520,6 +520,45 @@ function /(a::STaylor1{N,T}, b::T) where {N, T<:Number} STaylor1{N,T}(div_tuple_by_scalar(a.coeffs, b)) end +@generated function /(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N,T<:Number} + + ex_calc = quote end + append!(ex_calc.args, Any[nothing for i in 1:N]) + syms = Symbol[Symbol("c$i") for i in 1:N] + + # add error + ex_line = quote + if iszero(b[0]) + throw(ArgumentError("""The 0th order STaylor1 coefficient + must be non-zero for b, (a/b)(x) is + not differentiable at x=0).""")) + end + end + ex_calc.args[1] = ex_line + + # add recursion relation + for j = 0:(N-1) + ex_line = :(a[$(j)]) + for k = 1:j + sym = syms[j-k+1] + ex_line = :($ex_line - $sym*b[$k]) + end + sym = syms[j+1] + ex_line = :($sym = ($ex_line)/b[0]) + ex_calc.args[j+2] = ex_line + end + + exout = :(($(syms[1]),)) + for i = 2:N + push!(exout.args, syms[i]) + end + return quote + Base.@_inline_meta + $ex_calc + return STaylor1{N,T}($exout) + end +end + /(a::TaylorN{T}, b::TaylorN{S}) where {T<:NumberNotSeriesN, S<:NumberNotSeriesN} = /(promote(a,b)...) diff --git a/src/other_functions.jl b/src/other_functions.jl index bacdfdc3..5dc66b1d 100644 --- a/src/other_functions.jl +++ b/src/other_functions.jl @@ -21,7 +21,7 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @eval isnan(a::$T) = any( isnan.(a.coeffs) ) end for f in (:real, :imag, :conj) - @eval ($f)(a::STaylor1{N,T}) where {N,T<:Number} = STaylor1(($f).(a.coeffs), Val(N)) + @eval ($f)(a::STaylor1{N,T}) where {N,T<:Number} = STaylor1{N,T}(($f).(a.coeffs)) end adjoint(a::STaylor1) = conj(a) @@ -101,6 +101,18 @@ for T in (:Taylor1, :TaylorN) end end +function abs(a::STaylor1{N,T}) where {N,T<:Real} + if a[0] > zero(T) + return a + elseif a[0] < zero(T) + return -a + else + throw(ArgumentError( + """The 0th order Taylor1 coefficient must be non-zero + (abs(x) is not differentiable at x=0).""")) + end +end + function mod2pi(a::TaylorN{Taylor1{T}}) where {T<:Real} coeffs = copy(a.coeffs) @inbounds coeffs[1] = mod2pi( constant_term(a) ) @@ -285,6 +297,12 @@ for T in (:Taylor1, :TaylorN) @eval rad2deg(z::$T{T}) where {T<:Real} = z * (180 / convert(float(T), pi)) end +deg2rad(z::STaylor1{N, T}) where {N, T<:AbstractFloat} = z * (convert(T, pi) / 180) +deg2rad(z::STaylor1{N, T}) where {N, T<:Real} = z * (convert(float(T), pi) / 180) + +rad2deg(z::STaylor1{N, T}) where {N, T<:AbstractFloat} = z * (180 / convert(T, pi)) +rad2deg(z::STaylor1{N, T}) where {N, T<:Real} = z * (180 / convert(float(T), pi)) + # Internal mutating deg2rad!, rad2deg! functions for T in (:Taylor1, :TaylorN) @eval @inline function deg2rad!(v::$T{T}, a::$T{T}, k::Int) where {T<:AbstractFloat} diff --git a/test/onevariable.jl b/test/onevariable.jl index 8ffe2874..7e655205 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -580,11 +580,34 @@ end # check that STaylor1 and Taylor yeild same result t1 = STaylor1([1.1, 2.1, 3.1]) t2 = Taylor1([1.1, 2.1, 3.1]) - @test test_vs_Taylor1(exp(t1), exp(t2)) + for f in (exp, abs) + @test test_vs_Taylor1(f(t1), f(t2)) + end + t1a = STaylor1([2.1, 2.1, 3.1]) + t2a = Taylor1([2.1, 2.1, 3.1]) + @test isapprox((t1/t1a)[0], (t2/t2a)[0], atol=1E-10) + @test isapprox((t1/t1a)[1], (t2/t2a)[1], atol=1E-10) + @test isapprox((t1/t1a)[2], (t2/t2a)[2], atol=1E-10) + + @test isapprox((t1*t1a)[0], (t2*t2a)[0], atol=1E-10) + @test isapprox((t1*t1a)[1], (t2*t2a)[1], atol=1E-10) + @test isapprox((t1*t1a)[2], (t2*t2a)[2], atol=1E-10) a = STaylor1([0.0, 1.2, 2.3, 4.5, 0.0]) @test findfirst(a) == 1 @test findlast(a) == 3 + + a = STaylor1([5.0, 1.2, 2.3, 4.5, 0.0]) + @test isapprox(deg2rad(a)[0], 0.087266, atol=1E-5) + @test isapprox(deg2rad(a)[2], 0.040142, atol=1E-5) + @test isapprox(rad2deg(a)[0], 286.4788975, atol=1E-5) + @test isapprox(rad2deg(a)[2], 131.7802928, atol=1E-5) + @test real(a) == STaylor1([5.0, 1.2, 2.3, 4.5, 0.0]) + @test imag(a) == STaylor1([0.0, 0.0, 0.0, 0.0, 0.0]) + @test adjoint(a) == STaylor1([5.0, 1.2, 2.3, 4.5, 0.0]) + @test conj(a) == STaylor1([5.0, 1.2, 2.3, 4.5, 0.0]) + @test a == abs(a) + @test a == abs(-a) end #= From 66fa8ce669199cd046a8d3d27db17b4c8b16ea42 Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Fri, 10 Apr 2020 00:58:06 -0400 Subject: [PATCH 14/22] Enable non-STaylor1 unit tests --- test/onevariable.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/test/onevariable.jl b/test/onevariable.jl index 7e655205..5fe6d1ca 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -7,7 +7,6 @@ using Test using LinearAlgebra, SparseArrays eeuler = Base.MathConstants.e -#= @testset "Tests for Taylor1 expansions" begin ta(a) = Taylor1([a,one(a)],15) t = Taylor1(Int,15) @@ -523,7 +522,6 @@ eeuler = Base.MathConstants.e @test Taylor1{Int}(true) == Taylor1([1]) @test Taylor1{Int}(false) == Taylor1([0]) end -=# function test_vs_Taylor1(x,y) flag = true @@ -610,7 +608,6 @@ end @test a == abs(-a) end -#= @testset "Test `inv` for `Matrix{Taylor1{Float64}}``" begin t = Taylor1(5) a = Diagonal(rand(0:10,3)) + rand(3, 3) @@ -680,4 +677,3 @@ end @test_throws DimensionMismatch mul!(Y[1:end-1],A,B) end end -=# From 71f8fc89752e165c6fa2a57f31b1356ffed1973e Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Tue, 14 Apr 2020 18:34:30 -0400 Subject: [PATCH 15/22] Simplify computations --- src/arithmetic.jl | 46 ++++++++++------------------------------------ 1 file changed, 10 insertions(+), 36 deletions(-) diff --git a/src/arithmetic.jl b/src/arithmetic.jl index e4301be1..befa15e0 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -33,16 +33,8 @@ end for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @eval iszero(a::$T) = iszero(a.coeffs) end -function iszero(a::STaylor1) - flag = true - for i in eachindex(a) - flag &= iszero(a[i]) - if ~flag - break - end - end - flag -end + +iszero(a::STaylor1) = all(iszero, a.coeffs) ## zero and one ## for T in (:Taylor1, :TaylorN), f in (:zero, :one) @@ -219,37 +211,19 @@ end -(a::Taylor1{T}, b::TaylorN{S}) where {T<:NumberNotSeries,S<:NumberNotSeries} = -(promote(a,b)...) -@inline +(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(add_tuples(a.coeffs, b.coeffs)) -@inline -(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(sub_tuples(a.coeffs, b.coeffs)) +@inline +(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(a.coeffs .+ b.coeffs) +@inline -(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(a.coeffs .+ b.coeffs) @inline +(a::STaylor1) = a @inline -(a::STaylor1) = STaylor1(minus_tuple(a.coeffs)) -@generated function add_first_tuple(tup::NTuple{N,T}, x::T) where {N, T<:Number} - ex = :((tup[1] + x),) - append!(ex.args, Any[:(tup[$i]) for i in 2:N]) - return quote - Base.@_inline_meta - @inbounds return $ex - end -end - -@generated function sub_first_tuple(tup::NTuple{N,T}, x::T) where {N, T<:Number} - ex = :((tup[1] - x),) - append!(ex.args, Any[:(tup[$i]) for i in 2:N]) - return quote - Base.@_inline_meta - @inbounds return $ex - end -end - function +(a::STaylor1{N,T}, b::T) where {N, T<:Number} - STaylor1{N,T}(add_first_tuple(a.coeffs, b)) + STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] + b : a.coeffs[i], N)) end function +(b::T, a::STaylor1{N,T}) where {N, T<:Number} - STaylor1{N,T}(add_first_tuple(a.coeffs, b)) + STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] + b : a.coeffs[i], N)) end function -(a::STaylor1{N,T}, b::T) where {N, T<:Number} - STaylor1{N,T}(sub_first_tuple(a.coeffs, b)) + STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] - b : a.coeffs[i], N)) end -(b::T, a::STaylor1{N,T}) where {N, T<:Number} = b + (-a) @@ -307,10 +281,10 @@ end function *(a::STaylor1{N,T}, b::T) where {N, T<:Number} - STaylor1{N,T}(scale_tuple(a.coeffs, b)) + STaylor1{N,T}(b .* a.coeffs) end function *(b::T, a::STaylor1{N,T}) where {N, T<:Number} - STaylor1{N,T}(scale_tuple(a.coeffs, b)) + STaylor1{N,T}(b .* a.coeffs) end for T in (:HomogeneousPolynomial, :TaylorN) @@ -517,7 +491,7 @@ function /(a::Taylor1{T}, b::Taylor1{T}) where {T<:Number} end function /(a::STaylor1{N,T}, b::T) where {N, T<:Number} - STaylor1{N,T}(div_tuple_by_scalar(a.coeffs, b)) + STaylor1{N,T}(a.coeffs ./ b) end @generated function /(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N,T<:Number} From 409c2aed883384801c033ad3d936005516ec04b9 Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Tue, 14 Apr 2020 18:47:42 -0400 Subject: [PATCH 16/22] Fix failing test --- src/arithmetic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/arithmetic.jl b/src/arithmetic.jl index befa15e0..b141bc75 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -212,7 +212,7 @@ end -(promote(a,b)...) @inline +(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(a.coeffs .+ b.coeffs) -@inline -(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(a.coeffs .+ b.coeffs) +@inline -(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(a.coeffs .- b.coeffs) @inline +(a::STaylor1) = a @inline -(a::STaylor1) = STaylor1(minus_tuple(a.coeffs)) From e8c6121d96c01b0e935c930dea27f3e34cc8f5ed Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Wed, 15 Apr 2020 09:46:42 -0400 Subject: [PATCH 17/22] Simplify Benchmarking Suite --- benchmark/benchmarks.jl | 22 ++++++---------------- 1 file changed, 6 insertions(+), 16 deletions(-) diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index c765e4be..6cdc76ce 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -35,22 +35,12 @@ for i in dims q = STaylor1(r) q2 = STaylor1(r) y = rand() - 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} * 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) - 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) + for g in (+, -, /, *) + S["Taylor1{$i,Float64} $(Symbol(g)) Float64"] = @benchmarkable f($g, $x, $y, $n) + S["STaylor1{$i,Float64} $(Symbol(g)) Float64"] = @benchmarkable f($g, $q, $y, $n) + S["Taylor1{$i,Float64} $(Symbol(g)) Taylor1{$i,Float64}"] = @benchmarkable f($g, $x, $x2, $n) + S["STaylor1{$i,Float64} $(Symbol(g)) STaylor1{$i,Float64}"] = @benchmarkable f($g, $q, $q2, $n) + end end S = SUITE["functions"] = BenchmarkGroup() From 840299711d4ca098805ca45045c0258595b1581a Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Wed, 15 Apr 2020 21:21:30 -0400 Subject: [PATCH 18/22] STaylor1: Add square function --- src/power.jl | 36 ++++++++++++++++++++++++++++++++++++ test/onevariable.jl | 2 +- 2 files changed, 37 insertions(+), 1 deletion(-) diff --git a/src/power.jl b/src/power.jl index 52a6636b..9914cbe9 100644 --- a/src/power.jl +++ b/src/power.jl @@ -283,6 +283,42 @@ function square(a::HomogeneousPolynomial) return res end +@generated function square(a::STaylor1{N,T}) where {N,T} + ex_calc = quote end + append!(ex_calc.args, Any[nothing for i in 1:N]) + syms = Symbol[Symbol("c$i") for i in 1:N] + + sym = syms[1] + ex_line = :($(syms[1]) = a[0]^2) + ex_calc.args[1] = ex_line + + for k in 1:(N-1) + two = convert(T,2) + kodd = k%2 + kend = div(k - 2 + kodd, 2) + ex_line = :(a[0] * a[$k]) + @inbounds for i = 1:kend + ex_line = :($ex_line + a[$i] * a[$(k-i)]) + end + ex_line = :($two*($ex_line)) + if kodd !== 1 + ex_line = :($ex_line + a[$(div(k,2))]^2) + end + ex_line = :($(syms[k+1]) = $ex_line) + ex_calc.args[k+1] = ex_line + end + + exout = :(($(syms[1]),)) + for i = 2:N + push!(exout.args, syms[i]) + end + return quote + Base.@_inline_meta + $ex_calc + return STaylor1{N,T}($exout) + end +end + # Homogeneous coefficients for square @doc doc""" diff --git a/test/onevariable.jl b/test/onevariable.jl index 5fe6d1ca..0044534a 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -578,7 +578,7 @@ end # check that STaylor1 and Taylor yeild same result t1 = STaylor1([1.1, 2.1, 3.1]) t2 = Taylor1([1.1, 2.1, 3.1]) - for f in (exp, abs) + for f in (exp, abs, TaylorSeries.square) @test test_vs_Taylor1(f(t1), f(t2)) end t1a = STaylor1([2.1, 2.1, 3.1]) From 38e90246dbfa7fd6ca4e326877133dfdf689fc7d Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Wed, 22 Apr 2020 20:33:00 -0400 Subject: [PATCH 19/22] Update version conversion and constructor --- Project.toml | 2 +- src/constructors.jl | 5 ++++- src/conversion.jl | 2 +- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index beddaa9c..82dfb77c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TaylorSeries" uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git" -version = "0.10.3" +version = "0.11.0" [deps] InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" diff --git a/src/constructors.jl b/src/constructors.jl index 9c7cef3f..d2d92680 100644 --- a/src/constructors.jl +++ b/src/constructors.jl @@ -98,7 +98,7 @@ 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] + y = Any[:(zero($T)) for i=1:N] tup = :((x,)) push!(tup.args, y...) return quote @@ -106,6 +106,9 @@ given `N` with constant term equal to `x`. STaylor1{(N+1),T}($tup) end end +function STaylor1(coeffs::Vector{T}, l::Val{L}, v::Val{N}) where {L,N,T<:Number} + STaylor1{(N+1),T}(ntuple(i -> (i > L+1) ? coeffs[i] : zero(T), N+1)) +end @inline function STaylor1(coeffs::Vector{T}) where {T<:Number} STaylor1{length(coeffs),T}(tuple(coeffs...)) end diff --git a/src/conversion.jl b/src/conversion.jl index d16f653e..fb191e5e 100644 --- a/src/conversion.jl +++ b/src/conversion.jl @@ -41,7 +41,7 @@ function convert(::Type{STaylor1{N,T}}, b::Array{S,1}) where {N,T<:Number, S<:Nu STaylor1{N,T}(convert(Array{T,1},b)) end convert(::Type{STaylor1{N,T}}, a::STaylor1{N,T}) where {N,T<:Number} = a -convert(::Type{STaylor1{N,T}}, b::S) where {N,T<:Number, S<:Number} = STaylor1(b, Val(N)) +convert(::Type{STaylor1{N,T}}, b::S) where {N, T<:Number, S<:Number} = STaylor1(convert(T,b), Val(N)) convert(::Type{STaylor1{N,T}}, b::T) where {N, T<:Number} = STaylor1(b, Val(N)) convert(::Type{HomogeneousPolynomial{T}}, a::HomogeneousPolynomial) where {T<:Number} = From 1778a92006d11df90f7c03803515cf266b9747e5 Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Thu, 23 Apr 2020 12:10:35 -0400 Subject: [PATCH 20/22] Fix allocations in ntuple --- src/arithmetic.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/arithmetic.jl b/src/arithmetic.jl index b141bc75..5f7d92b5 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -217,16 +217,17 @@ end @inline -(a::STaylor1) = STaylor1(minus_tuple(a.coeffs)) function +(a::STaylor1{N,T}, b::T) where {N, T<:Number} - STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] + b : a.coeffs[i], N)) + STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] + b : a.coeffs[i], Val(N))) end function +(b::T, a::STaylor1{N,T}) where {N, T<:Number} - STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] + b : a.coeffs[i], N)) + STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] + b : a.coeffs[i], Val(N))) end function -(a::STaylor1{N,T}, b::T) where {N, T<:Number} - STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] - b : a.coeffs[i], N)) + STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] - b : a.coeffs[i], Val(N))) end -(b::T, a::STaylor1{N,T}) where {N, T<:Number} = b + (-a) +#+(a::STaylor1{N,T}, b::S) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(a,b)...) #+(a::STaylor1{N,T}, b::STaylor1{N,S}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(a,b)...) #+(a::STaylor1{N,T}, b::S) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(a,b)...) #+(b::S, a::STaylor1{N,T}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(b,a)...) From ed558cd67bef73224f8649e52c627e3dedcfc02a Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Thu, 23 Apr 2020 15:13:28 -0400 Subject: [PATCH 21/22] STaylor1: add log --- src/functions.jl | 36 +++++++++++++++++++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/src/functions.jl b/src/functions.jl index 28f623b0..66034145 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -147,7 +147,7 @@ for T in (:Taylor1, :TaylorN) end end -# Functions for StaticTaylor +# Functions for STaylor1 @generated function exp(a::STaylor1{N,T}) where {N, T <: Number} ex_calc = quote end append!(ex_calc.args, Any[nothing for i in 1:N]) @@ -180,6 +180,40 @@ end end end +@generated function log(a::STaylor1{N,T}) where {N, T <: Number} + ex_calc = quote end + append!(ex_calc.args, Any[nothing for i in 1:N]) + syms = Symbol[Symbol("c$i") for i in 1:N] + + (N >= 1) && (ex_calc.args[1] = :($(syms[1]) = log(constant_term(a)))) + (N >= 2) && (ex_calc.args[2] = :($(syms[2]) = a[1]/constant_term(a))) + + for k in 2:(N-1) + ex_line = :($(k-1)*a[1]*$(syms[k])) + @inbounds for i = 2:k-1 + ex_line = :($ex_line + $(k-i)*a[$i] * $(syms[k+1-i])) + end + ex_line = :((a[$k] - ($ex_line)/$(convert(T,k)))/constant_term(a)) + ex_line = :($(syms[k+1]) = $ex_line) + ex_calc.args[k+1] = ex_line + end + + exout = :(($(syms[1]),)) + for i = 2:N + push!(exout.args, syms[i]) + end + + return quote + Base.@_inline_meta + iszero(constant_term(a)) && throw(ArgumentError(""" + The 0-th order `STaylor1` coefficient must be non-zero + in order to expand `log` around 0. + """)) + $ex_calc + return STaylor1{N,T}($exout) + end +end + # Recursive functions (homogeneous coefficients) for T in (:Taylor1, :TaylorN) @eval begin From e169d1f20217d345b8ff24ec9124b01b95b32953 Mon Sep 17 00:00:00 2001 From: Matthew Wilhelm Date: Wed, 16 Sep 2020 11:01:46 -0400 Subject: [PATCH 22/22] Remove STaylor1 --- src/TaylorSeries.jl | 2 +- src/arithmetic.jl | 107 ----------------------------------------- src/auxiliary.jl | 32 ------------ src/constructors.jl | 43 ----------------- src/conversion.jl | 23 --------- src/evaluate.jl | 31 ------------ src/functions.jl | 67 -------------------------- src/intervals.jl | 2 + src/other_functions.jl | 25 ---------- src/power.jl | 39 +-------------- src/printing.jl | 20 -------- src/static.jl | 63 ++++++++++++++++++++++++ test/onevariable.jl | 74 ---------------------------- 13 files changed, 68 insertions(+), 460 deletions(-) create mode 100644 src/static.jl diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl index 9f8f53b5..afcdb0ae 100644 --- a/src/TaylorSeries.jl +++ b/src/TaylorSeries.jl @@ -42,7 +42,7 @@ import Base: zero, one, zeros, ones, isinf, isnan, iszero, power_by_squaring, rtoldefault, isfinite, isapprox, rad2deg, deg2rad -export Taylor1, TaylorN, STaylor1, HomogeneousPolynomial, AbstractSeries +export Taylor1, TaylorN, HomogeneousPolynomial, AbstractSeries export getcoeff, derivative, integrate, differentiate, evaluate, evaluate!, inverse, set_taylor1_varname, diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 5f7d92b5..9680a9de 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -23,8 +23,6 @@ for T in (:Taylor1, :TaylorN) end end -==(a::STaylor1, b::STaylor1) = (a.coeffs == b.coeffs) - function ==(a::HomogeneousPolynomial, b::HomogeneousPolynomial) a.order == b.order && return a.coeffs == b.coeffs return iszero(a.coeffs) && iszero(b.coeffs) @@ -34,14 +32,10 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @eval iszero(a::$T) = iszero(a.coeffs) end -iszero(a::STaylor1) = all(iszero, a.coeffs) - ## zero and one ## for T in (:Taylor1, :TaylorN), f in (:zero, :one) @eval ($f)(a::$T) = $T(($f)(a[0]), a.order) end -zero(::STaylor1{N,T}) where {N, T<:Number} = STaylor1(zero(T),Val(N-1)) -one(::STaylor1{N,T}) where {N, T<:Number} = STaylor1(one(T),Val(N-1)) function zero(a::HomogeneousPolynomial{T}) where {T<:Number} v = zero.(a.coeffs) @@ -211,32 +205,6 @@ end -(a::Taylor1{T}, b::TaylorN{S}) where {T<:NumberNotSeries,S<:NumberNotSeries} = -(promote(a,b)...) -@inline +(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(a.coeffs .+ b.coeffs) -@inline -(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(a.coeffs .- b.coeffs) -@inline +(a::STaylor1) = a -@inline -(a::STaylor1) = STaylor1(minus_tuple(a.coeffs)) - -function +(a::STaylor1{N,T}, b::T) where {N, T<:Number} - STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] + b : a.coeffs[i], Val(N))) -end -function +(b::T, a::STaylor1{N,T}) where {N, T<:Number} - STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] + b : a.coeffs[i], Val(N))) -end -function -(a::STaylor1{N,T}, b::T) where {N, T<:Number} - STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] - b : a.coeffs[i], Val(N))) -end --(b::T, a::STaylor1{N,T}) where {N, T<:Number} = b + (-a) - -#+(a::STaylor1{N,T}, b::S) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(a,b)...) -#+(a::STaylor1{N,T}, b::STaylor1{N,S}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(a,b)...) -#+(a::STaylor1{N,T}, b::S) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(a,b)...) -#+(b::S, a::STaylor1{N,T}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(b,a)...) - -#-(a::STaylor1{N,T}, b::STaylor1{N,S}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = -(promote(a,b)...) -#-(a::STaylor1{N,T}, b::S) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = -(promote(a,b)...) -#-(b::S, a::STaylor1{N,T}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = -(promote(b,a)...) - - ## Multiplication ## for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @@ -256,38 +224,6 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) end end -@generated function *(x::STaylor1{N,T}, y::STaylor1{N,T}) where {T<:Number,N} - ex_calc = quote end - append!(ex_calc.args, Any[nothing for i in 1:N]) - syms = Symbol[Symbol("c$i") for i in 1:N] - for j = 1:N - ex_line = :(x.coeffs[1]*y.coeffs[$j]) - for k = 2:j - ex_line = :($ex_line + x.coeffs[$k]*y.coeffs[$(j-k+1)]) - end - sym = syms[j] - ex_line = :($sym = $ex_line) - ex_calc.args[j] = ex_line - end - exout = :(($(syms[1]),)) - for i = 2:N - push!(exout.args, syms[i]) - end - return quote - Base.@_inline_meta - $ex_calc - return STaylor1{N,T}($exout) - end -end - - -function *(a::STaylor1{N,T}, b::T) where {N, T<:Number} - STaylor1{N,T}(b .* a.coeffs) -end -function *(b::T, a::STaylor1{N,T}) where {N, T<:Number} - STaylor1{N,T}(b .* a.coeffs) -end - for T in (:HomogeneousPolynomial, :TaylorN) @eval begin @@ -491,49 +427,6 @@ 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}(a.coeffs ./ b) -end - -@generated function /(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N,T<:Number} - - ex_calc = quote end - append!(ex_calc.args, Any[nothing for i in 1:N]) - syms = Symbol[Symbol("c$i") for i in 1:N] - - # add error - ex_line = quote - if iszero(b[0]) - throw(ArgumentError("""The 0th order STaylor1 coefficient - must be non-zero for b, (a/b)(x) is - not differentiable at x=0).""")) - end - end - ex_calc.args[1] = ex_line - - # add recursion relation - for j = 0:(N-1) - ex_line = :(a[$(j)]) - for k = 1:j - sym = syms[j-k+1] - ex_line = :($ex_line - $sym*b[$k]) - end - sym = syms[j+1] - ex_line = :($sym = ($ex_line)/b[0]) - ex_calc.args[j+2] = ex_line - end - - exout = :(($(syms[1]),)) - for i = 2:N - push!(exout.args, syms[i]) - end - return quote - Base.@_inline_meta - $ex_calc - return STaylor1{N,T}($exout) - end -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 946979c5..541598e5 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -103,13 +103,6 @@ function setindex!(a::Taylor1{T}, x::Array{T,1}, u::StepRange{Int,Int}) where {T end end - -# getindex STaylor1 -getindex(a::STaylor1, n::Int) = a.coeffs[n+1] -getindex(a::STaylor1, u::UnitRange{Int}) = a.coeffs[u .+ 1] -getindex(a::STaylor1, c::Colon) = a.coeffs[c] -getindex(a::STaylor1, u::StepRange{Int,Int}) = a.coeffs[u .+ 1] - """ getcoeff(a, v) @@ -242,16 +235,6 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @inline axes(a::$T) = () end end -@inline iterate(a::STaylor1{N,T}, state=0) where {N, T<:Number} = state > N-1 ? nothing : (a.coeffs[state+1], state+1) -@inline firstindex(a::STaylor1) = 0 -@inline lastindex(a::STaylor1{N,T}) where {N, T<:Number} = N-1 -@inline eachindex(s::STaylor1{N,T}) where {N, T<:Number} = UnitRange(0, N-1) -@inline size(s::STaylor1{N,T}) where {N, T<:Number} = N -@inline length(s::STaylor1{N,T}) where {N, T<:Number} = N -@inline get_order(s::STaylor1{N,T}) where {N, T<:Number} = N - 1 -@inline eltype(s::STaylor1{N,T}) where {N, T<:Number} = T -@inline axes(a::STaylor1) = () - ## fixorder ## for T in (:Taylor1, :TaylorN) @@ -286,19 +269,6 @@ 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 for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @@ -324,8 +294,6 @@ and `TaylorN`. The fallback behavior is to return `a` itself if """ constant_term(a::Taylor1) = a[0] -constant_term(a::STaylor1) = a[0] - constant_term(a::TaylorN) = a[0][1] constant_term(a::Vector{T}) where {T<:Number} = constant_term.(a) diff --git a/src/constructors.jl b/src/constructors.jl index d2d92680..e4166f99 100644 --- a/src/constructors.jl +++ b/src/constructors.jl @@ -71,49 +71,6 @@ julia> Taylor1(Rational{Int}, 4) Taylor1(::Type{T}, order::Int) where {T<:Number} = Taylor1( [zero(T), one(T)], order) 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 ## - -""" - 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,)) - push!(tup.args, y...) - return quote - Base.@_inline_meta - STaylor1{(N+1),T}($tup) - end -end -function STaylor1(coeffs::Vector{T}, l::Val{L}, v::Val{N}) where {L,N,T<:Number} - STaylor1{(N+1),T}(ntuple(i -> (i > L+1) ? coeffs[i] : zero(T), N+1)) -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 """ HomogeneousPolynomial{T<:Number} <: AbstractSeries{T} diff --git a/src/conversion.jl b/src/conversion.jl index fb191e5e..01ffc4b9 100644 --- a/src/conversion.jl +++ b/src/conversion.jl @@ -28,22 +28,6 @@ convert(::Type{Taylor1{T}}, b::T) where {T<:Number} = Taylor1([b], 0) convert(::Type{Taylor1}, a::T) where {T<:Number} = Taylor1(a, 0) -# Conversion for STaylor1 -function convert(::Type{STaylor1{N,Rational{T}}}, a::STaylor1{N,S}) where {N,T<:Integer, S<:AbstractFloat} - STaylor1{N,T}(rationalize.(a[:])) -end -function convert(::Type{STaylor1{N,T}}, b::Array{T,1}) where {N,T<:Number} - @assert N == length(b) - STaylor1{N,T}(b) -end -function convert(::Type{STaylor1{N,T}}, b::Array{S,1}) where {N,T<:Number, S<:Number} - @assert N == length(b) - STaylor1{N,T}(convert(Array{T,1},b)) -end -convert(::Type{STaylor1{N,T}}, a::STaylor1{N,T}) where {N,T<:Number} = a -convert(::Type{STaylor1{N,T}}, b::S) where {N, T<:Number, S<:Number} = STaylor1(convert(T,b), Val(N)) -convert(::Type{STaylor1{N,T}}, b::T) where {N, T<:Number} = STaylor1(b, Val(N)) - convert(::Type{HomogeneousPolynomial{T}}, a::HomogeneousPolynomial) where {T<:Number} = HomogeneousPolynomial(convert(Array{T,1}, a.coeffs), a.order) @@ -186,13 +170,6 @@ promote_rule(::Type{Taylor1{T}}, ::Type{S}) where {T<:Number, S<:Number} = promote_rule(::Type{Taylor1{Taylor1{T}}}, ::Type{Taylor1{T}}) where {T<:Number} = Taylor1{Taylor1{T}} -promote_rule(::Type{STaylor1{N,T}}, ::Type{STaylor1{N,T}}) where {N, T<:Number} = STaylor1{N,T} -promote_rule(::Type{STaylor1{N,T}}, ::Type{STaylor1{N,S}}) where {N, T<:Number, S<:Number} = STaylor1{N, promote_type(T,S)} -promote_rule(::Type{STaylor1{N,T}}, ::Type{Array{T,1}}) where {N, T<:Number} = STaylor1{N,T} -promote_rule(::Type{STaylor1{N,T}}, ::Type{Array{S,1}}) where {N, T<:Number, S<:Number} = STaylor1{N,promote_type(T,S)} -promote_rule(::Type{STaylor1{N,T}}, ::Type{T}) where {N, T<:Number} = STaylor1{N,T} -promote_rule(::Type{STaylor1{N,T}}, ::Type{S}) where {N, T<:Number, S<:Number} = STaylor1{N,promote_type(T,S)} - promote_rule(::Type{HomogeneousPolynomial{T}}, ::Type{HomogeneousPolynomial{S}}) where {T<:Number, S<:Number} = HomogeneousPolynomial{promote_type(T,S)} diff --git a/src/evaluate.jl b/src/evaluate.jl index 199c126f..8b12d17b 100644 --- a/src/evaluate.jl +++ b/src/evaluate.jl @@ -31,24 +31,6 @@ function evaluate(a::Taylor1{T}, dx::S) where {T<:Number, S<:Number} end evaluate(a::Taylor1{T}) where {T<:Number} = a[0] -function evaluate(a::STaylor1{N,T}, dx::T) where {N, T<:Number} - @inbounds suma = a[N-1] - @inbounds for k in (N-1):-1:0 - suma = suma*dx + a[k] - end - suma -end -#= -function evaluate(a::STaylor1{N,T}, dx::S) where {N, T<:Number, S<:Number} - suma = a[N-1]*one(dx) - @inbounds for k in (N-1):-1:0 - suma = suma*dx + a[k] - end - suma -end -=# -evaluate(a::STaylor1{N,T}) where {N, T<:Number} = a[0] - """ evaluate(x, δt) @@ -61,11 +43,6 @@ evaluate(x::Union{Array{Taylor1{T}}, SubArray{Taylor1{T}}}, δt::S) where {T<:Number, S<:Number} = evaluate.(x, δt) evaluate(a::Union{Array{Taylor1{T}}, SubArray{Taylor1{T}}}) where {T<:Number} = evaluate.(a, zero(T)) -evaluate(x::Union{Array{STaylor1{N,T}}, SubArray{STaylor1{N,T}}}, δt::S) where - {N, T<:Number, S<:Number} = evaluate.(x, δt) -evaluate(a::Union{Array{STaylor1{N,T}}, SubArray{STaylor1{N,T}}}) where - {N, T<:Number} = evaluate.(a, zero(T)) - """ evaluate!(x, δt, x0) @@ -136,20 +113,12 @@ evaluate(p::Taylor1{T}, x::Array{S}) where {T<:Number, S<:Number} = (p::Taylor1)(x) = evaluate(p, x) (p::Taylor1)() = evaluate(p) -(p::STaylor1)(x) = evaluate(p, x) -(p::STaylor1)() = evaluate(p) - #function-like behavior for Vector{Taylor1} (p::Array{Taylor1{T}})(x) where {T<:Number} = evaluate.(p, x) (p::SubArray{Taylor1{T}})(x) where {T<:Number} = evaluate.(p, x) (p::Array{Taylor1{T}})() where {T<:Number} = evaluate.(p) (p::SubArray{Taylor1{T}})() where {T<:Number} = evaluate.(p) -(p::Array{STaylor1{N,T}})(x) where {N,T<:Number} = evaluate.(p, x) -(p::SubArray{STaylor1{N,T}})(x) where {N,T<:Number} = evaluate.(p, x) -(p::Array{STaylor1{N,T}})() where {N,T<:Number} = evaluate.(p) -(p::SubArray{STaylor1{N,T}})() where {N,T<:Number} = evaluate.(p) - ## Evaluation of multivariable function evaluate!(x::Array{TaylorN{T},1}, δx::Array{T,1}, x0::Array{T,1}) where {T<:Number} diff --git a/src/functions.jl b/src/functions.jl index 66034145..2bffb9de 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -147,73 +147,6 @@ for T in (:Taylor1, :TaylorN) end end -# Functions for STaylor1 -@generated function exp(a::STaylor1{N,T}) where {N, T <: Number} - ex_calc = quote end - append!(ex_calc.args, Any[nothing for i in 1:N]) - syms = Symbol[Symbol("c$i") for i in 1:N] - - sym = syms[1] - ex_line = :($(syms[1]) = exp(a[0])) - ex_calc.args[1] = ex_line - - for k in 1:(N-1) - kT = convert(T,k) - sym = syms[k+1] - ex_line = :($kT * a[$k] * $(syms[1])) - @inbounds for i = 1:k-1 - ex_line = :($ex_line + $(kT-i) * a[$(k-i)] * $(syms[i+1])) - end - ex_line = :(($ex_line)/$kT) - ex_line = :($sym = $ex_line) - ex_calc.args[k+1] = ex_line - end - - exout = :(($(syms[1]),)) - for i = 2:N - push!(exout.args, syms[i]) - end - return quote - Base.@_inline_meta - $ex_calc - return STaylor1{N,T}($exout) - end -end - -@generated function log(a::STaylor1{N,T}) where {N, T <: Number} - ex_calc = quote end - append!(ex_calc.args, Any[nothing for i in 1:N]) - syms = Symbol[Symbol("c$i") for i in 1:N] - - (N >= 1) && (ex_calc.args[1] = :($(syms[1]) = log(constant_term(a)))) - (N >= 2) && (ex_calc.args[2] = :($(syms[2]) = a[1]/constant_term(a))) - - for k in 2:(N-1) - ex_line = :($(k-1)*a[1]*$(syms[k])) - @inbounds for i = 2:k-1 - ex_line = :($ex_line + $(k-i)*a[$i] * $(syms[k+1-i])) - end - ex_line = :((a[$k] - ($ex_line)/$(convert(T,k)))/constant_term(a)) - ex_line = :($(syms[k+1]) = $ex_line) - ex_calc.args[k+1] = ex_line - end - - exout = :(($(syms[1]),)) - for i = 2:N - push!(exout.args, syms[i]) - end - - return quote - Base.@_inline_meta - iszero(constant_term(a)) && throw(ArgumentError(""" - The 0-th order `STaylor1` coefficient must be non-zero - in order to expand `log` around 0. - """)) - $ex_calc - return STaylor1{N,T}($exout) - end -end - # Recursive functions (homogeneous coefficients) for T in (:Taylor1, :TaylorN) @eval begin diff --git a/src/intervals.jl b/src/intervals.jl index ef280624..85d65a4b 100644 --- a/src/intervals.jl +++ b/src/intervals.jl @@ -146,3 +146,5 @@ function _normalize(a::TaylorN, I::IntervalBox{N,T}, ::Val{false}) where {N,T} end return a(x) end + +square(x::Interval{T}) where {T <: Real} = pow(x,2) diff --git a/src/other_functions.jl b/src/other_functions.jl index 5dc66b1d..d1f8ef3b 100644 --- a/src/other_functions.jl +++ b/src/other_functions.jl @@ -20,13 +20,6 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @eval isnan(a::$T) = any( isnan.(a.coeffs) ) end -for f in (:real, :imag, :conj) - @eval ($f)(a::STaylor1{N,T}) where {N,T<:Number} = STaylor1{N,T}(($f).(a.coeffs)) -end - -adjoint(a::STaylor1) = conj(a) -isinf(a::STaylor1) = any(isinf.(a.coeffs)) -isnan(a::STaylor1) = any(isnan.(a.coeffs)) ## Division functions: rem and mod ## for op in (:mod, :rem) @@ -101,18 +94,6 @@ for T in (:Taylor1, :TaylorN) end end -function abs(a::STaylor1{N,T}) where {N,T<:Real} - if a[0] > zero(T) - return a - elseif a[0] < zero(T) - return -a - else - throw(ArgumentError( - """The 0th order Taylor1 coefficient must be non-zero - (abs(x) is not differentiable at x=0).""")) - end -end - function mod2pi(a::TaylorN{Taylor1{T}}) where {T<:Real} coeffs = copy(a.coeffs) @inbounds coeffs[1] = mod2pi( constant_term(a) ) @@ -297,12 +278,6 @@ for T in (:Taylor1, :TaylorN) @eval rad2deg(z::$T{T}) where {T<:Real} = z * (180 / convert(float(T), pi)) end -deg2rad(z::STaylor1{N, T}) where {N, T<:AbstractFloat} = z * (convert(T, pi) / 180) -deg2rad(z::STaylor1{N, T}) where {N, T<:Real} = z * (convert(float(T), pi) / 180) - -rad2deg(z::STaylor1{N, T}) where {N, T<:AbstractFloat} = z * (180 / convert(T, pi)) -rad2deg(z::STaylor1{N, T}) where {N, T<:Real} = z * (180 / convert(float(T), pi)) - # Internal mutating deg2rad!, rad2deg! functions for T in (:Taylor1, :TaylorN) @eval @inline function deg2rad!(v::$T{T}, a::$T{T}, k::Int) where {T<:AbstractFloat} diff --git a/src/power.jl b/src/power.jl index 9914cbe9..a5cc05c9 100644 --- a/src/power.jl +++ b/src/power.jl @@ -150,7 +150,6 @@ function ^(a::TaylorN, r::S) where {S<:Real} return c end - # Homogeneous coefficients for real power @doc doc""" pow!(c, a, r::Real, k::Int) @@ -283,42 +282,8 @@ function square(a::HomogeneousPolynomial) return res end -@generated function square(a::STaylor1{N,T}) where {N,T} - ex_calc = quote end - append!(ex_calc.args, Any[nothing for i in 1:N]) - syms = Symbol[Symbol("c$i") for i in 1:N] - - sym = syms[1] - ex_line = :($(syms[1]) = a[0]^2) - ex_calc.args[1] = ex_line - - for k in 1:(N-1) - two = convert(T,2) - kodd = k%2 - kend = div(k - 2 + kodd, 2) - ex_line = :(a[0] * a[$k]) - @inbounds for i = 1:kend - ex_line = :($ex_line + a[$i] * a[$(k-i)]) - end - ex_line = :($two*($ex_line)) - if kodd !== 1 - ex_line = :($ex_line + a[$(div(k,2))]^2) - end - ex_line = :($(syms[k+1]) = $ex_line) - ex_calc.args[k+1] = ex_line - end - - exout = :(($(syms[1]),)) - for i = 2:N - push!(exout.args, syms[i]) - end - return quote - Base.@_inline_meta - $ex_calc - return STaylor1{N,T}($exout) - end -end - +#two(x::T) = convert(T, x) +square(x::T) where {T <: Real} = x^2 # Homogeneous coefficients for square @doc doc""" diff --git a/src/printing.jl b/src/printing.jl index 23b9104c..99319b76 100644 --- a/src/printing.jl +++ b/src/printing.jl @@ -229,23 +229,3 @@ function show(io::IO, a::Union{Taylor1, HomogeneousPolynomial, TaylorN}) return print(io, pretty_print(a)) end end - -# printing for static taylor -function coeffstring(t::STaylor1, i, variable=:t) - if i == 1 # order 0 - return string(t.coeffs[i]) - end - if i == 2 # order 1 - return string(t.coeffs[i], variable) - end - return string(t.coeffs[i], variable, "^", i-1) -end - -function print_taylor(io::IO, t::STaylor1, variable=:t) - print(io, "(" * join([coeffstring(t, i, variable) for i in 1:length(t.coeffs)], " + ") * ")") -end - -Base.show(io::IO, t::STaylor1) = print_taylor(io, t) -function Base.show(io::IO, t::STaylor1{N,T}) where {N, T<:STaylor1} - print_taylor(io, t, :t) -end diff --git a/src/static.jl b/src/static.jl new file mode 100644 index 00000000..76800b71 --- /dev/null +++ b/src/static.jl @@ -0,0 +1,63 @@ +==(a::STaylor1, b::STaylor1) = (a.coeffs == b.coeffs) + +iszero(a::STaylor1) = all(iszero, a.coeffs) + +zero(::STaylor1{N,T}) where {N, T<:Number} = STaylor1(zero(T), Val(N-1)) +one(::STaylor1{N,T}) where {N, T<:Number} = STaylor1(one(T), Val(N-1)) + +@inline +(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(a.coeffs .+ b.coeffs) +@inline -(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(a.coeffs .- b.coeffs) +@inline +(a::STaylor1) = a +@inline -(a::STaylor1) = STaylor1(minus_tuple(a.coeffs)) + +function +(a::STaylor1{N,T}, b::T) where {N, T<:Number} + STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] + b : a.coeffs[i], Val(N))) +end +function +(b::T, a::STaylor1{N,T}) where {N, T<:Number} + STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] + b : a.coeffs[i], Val(N))) +end +function -(a::STaylor1{N,T}, b::T) where {N, T<:Number} + STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] - b : a.coeffs[i], Val(N))) +end +-(b::T, a::STaylor1{N,T}) where {N, T<:Number} = b + (-a) + +#+(a::STaylor1{N,T}, b::S) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(a,b)...) +#+(a::STaylor1{N,T}, b::STaylor1{N,S}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(a,b)...) +#+(a::STaylor1{N,T}, b::S) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(a,b)...) +#+(b::S, a::STaylor1{N,T}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(b,a)...) + +#-(a::STaylor1{N,T}, b::STaylor1{N,S}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = -(promote(a,b)...) +#-(a::STaylor1{N,T}, b::S) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = -(promote(a,b)...) +#-(b::S, a::STaylor1{N,T}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = -(promote(b,a)...) + +@generated function *(x::STaylor1{N,T}, y::STaylor1{N,T}) where {T<:Number,N} + ex_calc = quote end + append!(ex_calc.args, Any[nothing for i in 1:N]) + syms = Symbol[Symbol("c$i") for i in 1:N] + for j = 1:N + ex_line = :(x.coeffs[1]*y.coeffs[$j]) + for k = 2:j + ex_line = :($ex_line + x.coeffs[$k]*y.coeffs[$(j-k+1)]) + end + sym = syms[j] + ex_line = :($sym = $ex_line) + ex_calc.args[j] = ex_line + end + exout = :(($(syms[1]),)) + for i = 2:N + push!(exout.args, syms[i]) + end + return quote + Base.@_inline_meta + $ex_calc + return STaylor1{N,T}($exout) + end +end + + +function *(a::STaylor1{N,T}, b::T) where {N, T<:Number} + STaylor1{N,T}(b .* a.coeffs) +end +function *(b::T, a::STaylor1{N,T}) where {N, T<:Number} + STaylor1{N,T}(b .* a.coeffs) +end diff --git a/test/onevariable.jl b/test/onevariable.jl index 0044534a..c3e52c20 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -534,80 +534,6 @@ function test_vs_Taylor1(x,y) flag end -@testset "Tests for STaylor1 expansions" begin - - @test STaylor1 <: AbstractSeries - @test STaylor1{1,Float64} <: AbstractSeries{Float64} - @test STaylor1([1.0, 2.0]) == STaylor1((1.0, 2.0)) - @test STaylor1(STaylor1((1.0, 2.0))) == STaylor1((1.0, 2.0)) - @test STaylor1(1.0, Val(2)) == STaylor1((1.0, 0.0, 0.0)) - - - @test +STaylor1([1.0, 2.0, 3.0]) == STaylor1([1.0, 2.0, 3.0]) - @test -STaylor1([1.0, 2.0, 3.0]) == -STaylor1([1.0, 2.0, 3.0]) - @test STaylor1([1.0, 2.0, 3.0]) + STaylor1([3.0, 2.0, 3.0]) == STaylor1([4.0, 4.0, 6.0]) - @test STaylor1([1.0, 2.0, 3.0]) - STaylor1([3.0, 2.0, 4.0]) == STaylor1([-2.0, 0.0, -1.0]) - @test STaylor1([1.0, 2.0, 3.0]) + 2.0 == STaylor1([3.0, 2.0, 3.0]) - @test STaylor1([1.0, 2.0, 3.0]) - 2.0 == STaylor1([-1.0, 2.0, 3.0]) - @test 2.0 + STaylor1([1.0, 2.0, 3.0]) == STaylor1([3.0, 2.0, 3.0]) - @test 2.0 - STaylor1([1.0, 2.0, 3.0]) == STaylor1([1.0, -2.0, -3.0]) - - @test zero(STaylor1([1.0, 2.0, 3.0])) == STaylor1([0.0, 0.0, 0.0]) - @test one(STaylor1([1.0, 2.0, 3.0])) == STaylor1([1.0, 0.0, 0.0]) - - @test isinf(STaylor1([Inf, 2.0, 3.0])) && ~isinf(STaylor1([0.0, 0.0, 0.0])) - @test isnan(STaylor1([NaN, 2.0, 3.0])) && ~isnan(STaylor1([1.0, 0.0, 0.0])) - @test iszero(STaylor1([0.0, 0.0, 0.0])) && ~iszero(STaylor1([0.0, 1.0, 0.0])) - - @test length(STaylor1([0.0, 0.0, 0.0])) == 3 - @test size(STaylor1([0.0, 0.0, 0.0])) == 3 - @test firstindex(STaylor1([0.0, 0.0, 0.0])) == 0 - @test lastindex(STaylor1([0.0, 0.0, 0.0])) == 2 - - st1 = STaylor1([1.0, 2.0, 3.0]) - @test st1(2.0) == 41.0 - @test st1() == 1.00 - st2 = typeof(st1)[st1; st1] - @test st2(2.0)[1] == st2(2.0)[2] == 41.0 - @test st2()[1] == st2()[2] == 1.0 - @test evaluate(st1,2.0) == 41.0 - @test evaluate(st1) == 1.00 - @test evaluate(st2,2.0)[1] == evaluate(st2,2.0)[2] == 41.0 - @test evaluate(st2)[1] == evaluate(st2)[2] == 1.0 - - # check that STaylor1 and Taylor yeild same result - t1 = STaylor1([1.1, 2.1, 3.1]) - t2 = Taylor1([1.1, 2.1, 3.1]) - for f in (exp, abs, TaylorSeries.square) - @test test_vs_Taylor1(f(t1), f(t2)) - end - t1a = STaylor1([2.1, 2.1, 3.1]) - t2a = Taylor1([2.1, 2.1, 3.1]) - @test isapprox((t1/t1a)[0], (t2/t2a)[0], atol=1E-10) - @test isapprox((t1/t1a)[1], (t2/t2a)[1], atol=1E-10) - @test isapprox((t1/t1a)[2], (t2/t2a)[2], atol=1E-10) - - @test isapprox((t1*t1a)[0], (t2*t2a)[0], atol=1E-10) - @test isapprox((t1*t1a)[1], (t2*t2a)[1], atol=1E-10) - @test isapprox((t1*t1a)[2], (t2*t2a)[2], atol=1E-10) - - a = STaylor1([0.0, 1.2, 2.3, 4.5, 0.0]) - @test findfirst(a) == 1 - @test findlast(a) == 3 - - a = STaylor1([5.0, 1.2, 2.3, 4.5, 0.0]) - @test isapprox(deg2rad(a)[0], 0.087266, atol=1E-5) - @test isapprox(deg2rad(a)[2], 0.040142, atol=1E-5) - @test isapprox(rad2deg(a)[0], 286.4788975, atol=1E-5) - @test isapprox(rad2deg(a)[2], 131.7802928, atol=1E-5) - @test real(a) == STaylor1([5.0, 1.2, 2.3, 4.5, 0.0]) - @test imag(a) == STaylor1([0.0, 0.0, 0.0, 0.0, 0.0]) - @test adjoint(a) == STaylor1([5.0, 1.2, 2.3, 4.5, 0.0]) - @test conj(a) == STaylor1([5.0, 1.2, 2.3, 4.5, 0.0]) - @test a == abs(a) - @test a == abs(-a) -end - @testset "Test `inv` for `Matrix{Taylor1{Float64}}``" begin t = Taylor1(5) a = Diagonal(rand(0:10,3)) + rand(3, 3)