From 1d513985e5b8eced87669ed23ee3a79331d6c0e6 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Sat, 9 Feb 2019 23:01:52 -0600 Subject: [PATCH 1/2] Avoid duplicating order in product of `TaylorModelN`s --- src/arithmetic.jl | 57 ++++++++++++++++++++++++++++++----------------- test/TMN.jl | 9 +++----- 2 files changed, 40 insertions(+), 26 deletions(-) diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 31b9db49..3ea8dbab 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -208,30 +208,47 @@ end # Multiplication function *(a::TaylorModelN, b::TaylorModelN) @assert a.x0 == b.x0 && a.I == b.I - - # Polynomial product with extended order - order = max(get_order(a), get_order(b)) - @assert 2*order ≤ get_order() - aext = TaylorN(copy(a.pol.coeffs), 2*order) - bext = TaylorN(copy(b.pol.coeffs), 2*order) - res = aext * bext - - # Returned polynomial - bext = TaylorN( copy(res.coeffs[1:order+1]) ) - - # Bound for the neglected part of the product of polynomials - res[0:order] .= zero(eltype(res)) aux = a.I - a.x0 - Δnegl = res(aux) - # Remainder of the product - Δa = a.pol(aux) - Δb = b.pol(aux) - Δ = Δnegl + Δb * a.rem + Δa * b.rem + a.rem * b.rem + # Polynomial product with largest order from TaylorSeries._params_TaylorN_ + order_a = get_order(a) + order_b = get_order(b) + order_max = max(order_a, order_b) + order_prod = order_a + order_b + orderTS = get_order() + + # The returned polynomial has no neglected part + if order_prod ≤ orderTS + apol = TaylorN(a.pol.coeffs, order_prod) + bpol = TaylorN(b.pol.coeffs, order_prod) + res = apol * bpol + Δa = a.pol(aux) + Δb = b.pol(aux) + Δ = Δb * a.rem + Δa * b.rem + a.rem * b.rem + return TaylorModelN(res, Δ, a.x0, a.I) + end - return TaylorModelN(bext, Δ, a.x0, a.I) -end + # The returned polynomial has a neglected part + # Bound the neglected part of the product of polynomials + apol = TaylorN(a.pol.coeffs, orderTS) + bpol = TaylorN(b.pol.coeffs, orderTS) + res = apol * bpol + Δa = apol(aux) + Δb = bpol(aux) + Δ = Δb * a.rem + Δa * b.rem + a.rem * b.rem + + # We evaluate each term of the product of coefficients at `aux` + Δnegl = zero(Δ) + for order = order_max+1:order_prod + orderini = order - order_max + for inda = orderini:order_a + indb = order - inda + Δnegl += evaluate(apol[inda], aux) * evaluate(bpol[indb], aux) + end + end + return TaylorModelN(res, Δ+Δnegl, a.x0, a.I) +end # Multiplication by numbers function *(b::T, a::TaylorModelN) where {T<:NumberNotSeries} diff --git a/test/TMN.jl b/test/TMN.jl index 3ea3437d..c9b86371 100644 --- a/test/TMN.jl +++ b/test/TMN.jl @@ -20,8 +20,7 @@ function check_containment(ftest, xx::TaylorModelN{N,T,S}, tma::TaylorModelN{N,T end const _order = 2 -const _order_max = 2*(_order+1) -set_variables(Interval{Float64}, [:x, :y], order=_order_max) +set_variables(Interval{Float64}, [:x, :y], order=_order) @testset "Tests for TaylorModelN " begin b0 = Interval(0.0) × Interval(0.0) @@ -53,7 +52,7 @@ set_variables(Interval{Float64}, [:x, :y], order=_order_max) @test_throws BoundsError TaylorModelN(5, _order, b0, ib0) # wrong variable number # Tests for get_order and remainder - @test get_order() == 6 + @test get_order() == 2 @test get_order(xm) == 2 @test remainder(ym) == zi end @@ -79,8 +78,6 @@ set_variables(Interval{Float64}, [:x, :y], order=_order_max) b = b1[2] * a @test b == TaylorModelN( b1[2]*a.pol, Δ*b1[2], b1, ib1 ) @test b / b1[2] == a - @test_throws AssertionError TaylorModelN(TaylorN(1, order=_order_max), zi, b1, ib1) * - TaylorModelN(TaylorN(2, order=_order_max), zi, b1, ib1) remt = remainder(1/(1-TaylorModel1(_order, b1[1], ib1[1]))) @test remainder(1 / (1-xm)) == remt @@ -100,7 +97,7 @@ set_variables(Interval{Float64}, [:x, :y], order=_order_max) @test rpa(x->5*x, ym) == 5*ym remT = remainder(5*TaylorModel1(2, b1[1], ib1[1])^4) @test rpa(x->5*x^4, xm) == TaylorModelN(zero(xT), remT, b1, ib1) - remT = 5 * (ib1[1]-b1[1])^2 * (2*(ib1[2]-b1[2])+(ib1[2]-b1[2])^2) + remT = 5 * (ib1[1]-b1[1])^2 *(ib1[2]-b1[2]) * (2+(ib1[2]-b1[2])) @test rpa(x->5*x^2, xm*ym) == TaylorModelN( 5*xT^2, remT, b1, ib1) # Testing remainders of an RPA From 800c016a68a7a00427313a627bc011917751ade6 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Sun, 10 Feb 2019 12:14:08 -0600 Subject: [PATCH 2/2] Explicit evaluation of monomials, and revert a modified test --- src/arithmetic.jl | 45 +++++++++++++++++++++++++++++++++------------ test/TMN.jl | 3 ++- 2 files changed, 35 insertions(+), 13 deletions(-) diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 3ea8dbab..4787a178 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -210,14 +210,13 @@ function *(a::TaylorModelN, b::TaylorModelN) @assert a.x0 == b.x0 && a.I == b.I aux = a.I - a.x0 - # Polynomial product with largest order from TaylorSeries._params_TaylorN_ + # Some definitions related to the order of some polynomials + orderTS = get_order() # maximum order, fixed by `TaylorSeries._params_TaylorN_` order_a = get_order(a) order_b = get_order(b) - order_max = max(order_a, order_b) - order_prod = order_a + order_b - orderTS = get_order() + order_prod = 2 * max(order_a, order_b) - # The returned polynomial has no neglected part + # The returned polynomial has no terms beyond `orderTS` if order_prod ≤ orderTS apol = TaylorN(a.pol.coeffs, order_prod) bpol = TaylorN(b.pol.coeffs, order_prod) @@ -228,28 +227,50 @@ function *(a::TaylorModelN, b::TaylorModelN) return TaylorModelN(res, Δ, a.x0, a.I) end - # The returned polynomial has a neglected part - # Bound the neglected part of the product of polynomials + # The resulting product has a part whose order is larger than `orderTS`; + # a bound for this polynomial has to be included apol = TaylorN(a.pol.coeffs, orderTS) bpol = TaylorN(b.pol.coeffs, orderTS) - res = apol * bpol + res = apol * bpol # Trunctated polynomial to order `orderTS` Δa = apol(aux) Δb = bpol(aux) Δ = Δb * a.rem + Δa * b.rem + a.rem * b.rem + N = get_numvars() - # We evaluate each term of the product of coefficients at `aux` + # We evaluate *explicitely* each term of the product of coefficients at `aux` Δnegl = zero(Δ) - for order = order_max+1:order_prod - orderini = order - order_max + for order = orderTS+1:order_prod + orderini = order - orderTS for inda = orderini:order_a indb = order - inda - Δnegl += evaluate(apol[inda], aux) * evaluate(bpol[indb], aux) + # Δnegl += evaluate(apol[inda], aux) * evaluate(bpol[indb], aux) + + # Coefficient tables of corresponding `HomogeneousPolynomial`s + @inbounds cta = TaylorSeries.coeff_table[inda+1] + @inbounds ctb = TaylorSeries.coeff_table[indb+1] + + for i in eachindex(cta) + @inbounds a_coeff = apol.coeffs[inda+1][i] + iszero(a_coeff) && continue + for j in eachindex(ctb) + @inbounds b_coeff = bpol.coeffs[indb+1][j] + iszero(b_coeff) && continue + + # Evaluation of monomials + tmp = one(Δ) + for n in 1:N + @inbounds tmp *= aux[n]^(cta[i][n]+ctb[j][n]) + end + Δnegl += a_coeff * b_coeff * tmp + end + end end end return TaylorModelN(res, Δ+Δnegl, a.x0, a.I) end + # Multiplication by numbers function *(b::T, a::TaylorModelN) where {T<:NumberNotSeries} pol = a.pol * b diff --git a/test/TMN.jl b/test/TMN.jl index c9b86371..951e5791 100644 --- a/test/TMN.jl +++ b/test/TMN.jl @@ -97,7 +97,8 @@ set_variables(Interval{Float64}, [:x, :y], order=_order) @test rpa(x->5*x, ym) == 5*ym remT = remainder(5*TaylorModel1(2, b1[1], ib1[1])^4) @test rpa(x->5*x^4, xm) == TaylorModelN(zero(xT), remT, b1, ib1) - remT = 5 * (ib1[1]-b1[1])^2 *(ib1[2]-b1[2]) * (2+(ib1[2]-b1[2])) + # remT = 5 * (ib1[1]-b1[1])^2 *(ib1[2]-b1[2]) * (2+(ib1[2]-b1[2])) + remT = 5 * (ib1[1]-b1[1])^2 * (2*(ib1[2]-b1[2])+(ib1[2]-b1[2])^2) @test rpa(x->5*x^2, xm*ym) == TaylorModelN( 5*xT^2, remT, b1, ib1) # Testing remainders of an RPA