From 02c6abbaf8dc882448cb815702ad3b4b69863d23 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Tue, 12 Dec 2023 20:09:38 -0600 Subject: [PATCH] Further additions and clean up --- ext/TaylorSeriesIAExt.jl | 610 ++++++++++++++++++++++++++++----------- 1 file changed, 446 insertions(+), 164 deletions(-) diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl index f5b3a39c..14ebba97 100644 --- a/ext/TaylorSeriesIAExt.jl +++ b/ext/TaylorSeriesIAExt.jl @@ -401,92 +401,460 @@ end # several math functions for T in (:Taylor1, :TaylorN) - @eval function log(a::$T{Interval{S}}) where {S<:Real} - TS._isthinzero(constant_term(a)) && throw(DomainError(a, - """The 0-th order coefficient must be non-zero in order to expand `log` around 0.""")) - a0 = intersect_interval(constant_term(a), interval(zero(S), S(Inf))) - order = a.order - aux = log(a0) - aa = one(aux) * a - aa[0] = one(aux) * a0 - c = $T( aux, order ) - for k in eachindex(a) - TS.log!(c, aa, k) + @eval begin + function log(a::$T{Interval{S}}) where {S<:Real} + TS._isthinzero(constant_term(a)) && throw(DomainError(a, + """The 0-th order coefficient must be non-zero in order to expand `log` around 0.""")) + a0 = intersect_interval(constant_term(a), interval(zero(S), S(Inf))) + order = a.order + aux = log(a0) + aa = one(aux) * a + aa[0] = one(aux) * a0 + c = $T( aux, order ) + for k in eachindex(a) + TS.log!(c, aa, k) + end + return c end - return c - end - @eval function asin(a::$T{Interval{S}}) where {S<:Real} - a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) - a0sqr = a0^2 - isequal_interval(a0sqr, one(a0)) && - throw(DomainError(a, - """Series expansion of asin(x) diverges at x = ±1.""")) - - order = a.order - aux = asin(a0) - aa = one(aux) * a - aa[0] = one(aux) * a0 - c = $T( aux, order ) - r = $T( sqrt(1 - a0sqr), order ) - for k in eachindex(a) - TS.asin!(c, aa, r, k) + function asin(a::$T{Interval{S}}) where {S<:Real} + a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) + a0sqr = a0^2 + isequal_interval(a0sqr, one(a0)) && + throw(DomainError(a, + """Series expansion of asin(x) diverges at x = ±1.""")) + + order = a.order + aux = asin(a0) + aa = one(aux) * a + aa[0] = one(aux) * a0 + c = $T( aux, order ) + r = $T( sqrt(1 - a0sqr), order ) + for k in eachindex(a) + TS.asin!(c, aa, r, k) + end + return c end - return c - end - @eval function acos(a::$T{Interval{S}}) where {S<:Real} - a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) - a0sqr = a0^2 - isequal_interval(a0sqr, one(a0)) && - throw(DomainError(a, - """Series expansion of acos(x) diverges at x = ±1.""")) - - order = a.order - aux = acos(a0) - aa = one(aux) * a - aa[0] = one(aux) * a0 - c = $T( aux, order ) - r = $T( sqrt(1 - a0sqr), order ) - for k in eachindex(a) - TS.acos!(c, aa, r, k) + function acos(a::$T{Interval{S}}) where {S<:Real} + a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) + a0sqr = a0^2 + isequal_interval(a0sqr, one(a0)) && + throw(DomainError(a, + """Series expansion of acos(x) diverges at x = ±1.""")) + + order = a.order + aux = acos(a0) + aa = one(aux) * a + aa[0] = one(aux) * a0 + c = $T( aux, order ) + r = $T( sqrt(1 - a0sqr), order ) + for k in eachindex(a) + TS.acos!(c, aa, r, k) + end + return c end - return c - end - @eval function acosh(a::$T{Interval{S}}) where {S<:Real} - a0 = intersect_interval(constant_term(a), interval(one(S), S(Inf))) - a0sqr = a0^2 - isequal_interval(a0sqr, one(a0)) && throw(DomainError(a, - """Series expansion of acosh(x) diverges at x = ±1.""")) - - order = a.order - aux = acosh(a0) - aa = one(aux) * a - aa[0] = one(aux) * a0 - c = $T( aux, order ) - r = $T( sqrt(a0sqr - 1), order ) - for k in eachindex(a) - TS.acosh!(c, aa, r, k) + function acosh(a::$T{Interval{S}}) where {S<:Real} + a0 = intersect_interval(constant_term(a), interval(one(S), S(Inf))) + a0sqr = a0^2 + isequal_interval(a0sqr, one(a0)) && throw(DomainError(a, + """Series expansion of acosh(x) diverges at x = ±1.""")) + + order = a.order + aux = acosh(a0) + aa = one(aux) * a + aa[0] = one(aux) * a0 + c = $T( aux, order ) + r = $T( sqrt(a0sqr - 1), order ) + for k in eachindex(a) + TS.acosh!(c, aa, r, k) + end + return c + end + + function atanh(a::$T{Interval{S}}) where {S<:Real} + order = a.order + a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) + aux = atanh(a0) + aa = one(aux) * a + aa[0] = one(aux) * a0 + c = $T( aux, order) + r = $T(one(aux) - a0^2, order) + TS._isthinzero(constant_term(r)) && throw(DomainError(a, + """Series expansion of atanh(x) diverges at x = ±1.""")) + + for k in eachindex(a) + TS.atanh!(c, aa, r, k) + end + return c + end + + # Some internal functions + @inline function exp!(c::$T{Interval{T}}, a::$T{Interval{T}}, k::Int) where {T<:Real} + if k == 0 + @inbounds c[0] = exp(constant_term(a)) + return nothing + end + if $T == Taylor1 + @inbounds c[k] = interval(k) * a[k] * c[0] + else + @inbounds mul!(c[k], interval(k) * a[k], c[0]) + end + @inbounds for i = 1:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * a[k-i] * c[i] + else + mul!(c[k], interval(k-i) * a[k-i], c[i]) + end + end + @inbounds c[k] = c[k] / interval(k) + return nothing + end + + @inline function expm1!(c::$T{Interval{T}}, a::$T{Interval{T}}, k::Int) where {T<:Real} + if k == 0 + @inbounds c[0] = expm1(constant_term(a)) + return nothing + end + c0 = c[0]+one(c[0]) + if $T == Taylor1 + @inbounds c[k] = interval(k) * a[k] * c0 + else + @inbounds mul!(c[k], interval(k) * a[k], c0) + end + @inbounds for i = 1:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * a[k-i] * c[i] + else + mul!(c[k], interval(k-i) * a[k-i], c[i]) + end + end + @inbounds c[k] = c[k] / interval(k) + return nothing + end + + @inline function log!(c::$T{Interval{T}}, a::$T{T}, k::Int) where {T<:Real} + if k == 0 + @inbounds c[0] = log(constant_term(a)) + return nothing + elseif k == 1 + @inbounds c[1] = a[1] / constant_term(a) + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k-1) * a[1] * c[k-1] + else + @inbounds mul!(c[k], interval(k-1)*a[1], c[k-1]) + end + @inbounds for i = 2:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * a[i] * c[k-i] + else + mul!(c[k], interval(k-i)*a[i], c[k-i]) + end + end + @inbounds c[k] = (a[k] - c[k]/interval(k)) / constant_term(a) + return nothing + end + + @inline function log1p!(c::$T{Interval{T}}, a::$T{Interval{T}}, k::Int) where {T<:Real} + a0 = constant_term(a) + a0p1 = a0+one(a0) + if k == 0 + @inbounds c[0] = log1p(a0) + return nothing + elseif k == 1 + @inbounds c[1] = a[1] / a0p1 + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k-1) * a[1] * c[k-1] + else + @inbounds mul!(c[k], interval(k-1)*a[1], c[k-1]) + end + @inbounds for i = 2:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * a[i] * c[k-i] + else + mul!(c[k], interval(k-i)*a[i], c[k-i]) + end + end + @inbounds c[k] = (a[k] - c[k]/interval(k)) / a0p1 + return nothing + end + + @inline function sincos!(s::$T{Interval{T}}, c::$T{Interval{T}}, a::$T{Interval{T}}, + k::Int) where {T<:Real} + if k == 0 + a0 = constant_term(a) + @inbounds s[0], c[0] = sincos( a0 ) + return nothing + end + + x = a[1] + if $T == Taylor1 + @inbounds s[k] = x * c[k-1] + @inbounds c[k] = -x * s[k-1] + else + mul!(s[k], x, c[k-1]) + mul!(c[k], -x, s[k-1]) + end + @inbounds for i = 2:k + x = interval(i) * a[i] + if $T == Taylor1 + s[k] += x * c[k-i] + c[k] -= x * s[k-i] + else + mul!(s[k], x, c[k-i]) + mul!(c[k], -x, s[k-i]) + end + end + + @inbounds s[k] = s[k] / interval(k) + @inbounds c[k] = c[k] / interval(k) + return nothing + end + + @inline function tan!(c::$T{Interval{T}}, a::$T{Interval{T}}, c2::$T{Interval{T}}, + k::Int) where {T<:Real} + if k == 0 + @inbounds aux = tan( constant_term(a) ) + @inbounds c[0] = aux + @inbounds c2[0] = aux^2 + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k) * a[k] * c2[0] + else + @inbounds mul!(c[k], interval(k) * a[k], c2[0]) + end + @inbounds for i = 1:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * a[k-i] * c2[i] + else + mul!(c[k], interval(k-i) * a[k-i], c2[i]) + end + end + @inbounds c[k] = a[k] + c[k]/interval(k) + sqr!(c2, c, k) + + return nothing + end + + @inline function asin!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, + k::Int) where {T<:Real} + if k == 0 + a0 = constant_term(a) + @inbounds c[0] = asin( a0 ) + @inbounds r[0] = sqrt( one(a0) - a0^2 ) + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k-1) * r[1] * c[k-1] + else + @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + end + @inbounds for i in 2:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * r[i] * c[k-i] + else + mul!(c[k], interval(k-i) * r[i], c[k-i]) + end + end + sqrt!(r, one(a0)-a^2, k) + @inbounds c[k] = (a[k] - c[k]/interval(k)) / constant_term(r) + return nothing end - return c - end - @eval function atanh(a::$T{Interval{S}}) where {S<:Real} - order = a.order - a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) - aux = atanh(a0) - aa = one(aux) * a - aa[0] = one(aux) * a0 - c = $T( aux, order) - r = $T(one(aux) - a0^2, order) - TS._isthinzero(constant_term(r)) && throw(DomainError(a, - """Series expansion of atanh(x) diverges at x = ±1.""")) - - for k in eachindex(a) - TS.atanh!(c, aa, r, k) + @inline function acos!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, + k::Int) where {T<:Real} + if k == 0 + a0 = constant_term(a) + @inbounds c[0] = acos( a0 ) + @inbounds r[0] = sqrt( one(a0) - a0^2 ) + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k-1) * r[1] * c[k-1] + else + @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + end + @inbounds for i in 2:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * r[i] * c[k-i] + else + mul!(c[k], interval(k-i) * r[i], c[k-i]) + end + end + sqrt!(r, one(a0)-a^2, k) + @inbounds c[k] = -(a[k] + c[k]/interval(k)) / constant_term(r) + return nothing + end + + @inline function atan!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, + k::Int) where {T<:Real} + if k == 0 + a0 = constant_term(a) + @inbounds c[0] = atan( a0 ) + @inbounds r[0] = one(a0) + a0^2 + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k-1) * r[1] * c[k-1] + else + @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + end + @inbounds for i in 2:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * r[i] * c[k-i] + else + mul!(c[k], interval(k-i) * r[i], c[k-i]) + end + end + @inbounds sqr!(r, a, k) + @inbounds c[k] = (a[k] - c[k]/interval(k)) / constant_term(r) + return nothing + end + + @inline function sinhcosh!(s::$T{Interval{T}}, c::$T{Interval{T}}, a::$T{Interval{T}}, + k::Int) where {T<:Real} + if k == 0 + @inbounds s[0] = sinh( constant_term(a) ) + @inbounds c[0] = cosh( constant_term(a) ) + return nothing + end + + x = a[1] + if $T == Taylor1 + @inbounds s[k] = x * c[k-1] + @inbounds c[k] = x * s[k-1] + else + @inbounds mul!(s[k], x, c[k-1]) + @inbounds mul!(c[k], x, s[k-1]) + end + @inbounds for i = 2:k + x = interval(i) * a[i] + if $T == Taylor1 + s[k] += x * c[k-i] + c[k] += x * s[k-i] + else + mul!(s[k], x, c[k-i]) + mul!(c[k], x, s[k-i]) + end + end + s[k] = s[k] / interval(k) + c[k] = c[k] / interval(k) + return nothing + end + + @inline function tanh!(c::$T{Interval{T}}, a::$T{Interval{T}}, c2::$T{Interval{T}}, + k::Int) where {T<:Real} + if k == 0 + @inbounds aux = tanh( constant_term(a) ) + @inbounds c[0] = aux + @inbounds c2[0] = aux^2 + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = k * a[k] * c2[0] + else + @inbounds mul!(c[k], k * a[k], c2[0]) + end + @inbounds for i = 1:k-1 + if $T == Taylor1 + c[k] += (k-i) * a[k-i] * c2[i] + else + mul!(c[k], (k-i) * a[k-i], c2[i]) + end + end + @inbounds c[k] = a[k] - c[k]/k + sqr!(c2, c, k) + + return nothing + end + + @inline function asinh!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, k::Int) where {T<:Real} + if k == 0 + a0 = constant_term(a) + @inbounds c[0] = asinh( a0 ) + @inbounds r[0] = sqrt( a0^2 + one(a0) ) + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k-1) * r[1] * c[k-1] + else + @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + end + @inbounds for i in 2:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * r[i] * c[k-i] + else + mul!(c[k], interval(k-i) * r[i], c[k-i]) + end + end + sqrt!(r, a^2+one(a0), k) + @inbounds c[k] = (a[k] - c[k]/interval(k)) / constant_term(r) + return nothing + end + + @inline function acosh!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, k::Int) where {T<:Real} + if k == 0 + a0 = constant_term(a) + @inbounds c[0] = acosh( a0 ) + @inbounds r[0] = sqrt( a0^2 - one(a0) ) + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k-1) * r[1] * c[k-1] + else + @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + end + @inbounds for i in 2:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * r[i] * c[k-i] + else + mul!(c[k], interval(k-i) * r[i], c[k-i]) + end + end + sqrt!(r, a^2-one(a0), k) + @inbounds c[k] = (a[k] - c[k]/interval(k)) / constant_term(r) + return nothing + end + + @inline function atanh!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, k::Int) where {T<:Real} + if k == 0 + a0 = constant_term(a) + @inbounds c[0] = atanh( a0 ) + @inbounds r[0] = one(a0) - a0^2 + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k-1) * r[1] * c[k-1] + else + @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + end + @inbounds for i in 2:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * r[i] * c[k-i] + else + mul!(c[k], interval(k-i) * r[i], c[k-i]) + end + end + @inbounds sqr!(r, a, k) + @inbounds c[k] = (a[k] + c[k]/interval(k)) / constant_term(r) + return nothing end - return c end end @@ -668,93 +1036,7 @@ function _normalize(a::TaylorN, I::AbstractVector{Interval{T}}, ::Val{false}) wh return a(x) end -# Printing-related methods -# function TS.pretty_print(a::Taylor1{Interval{T}}) where {T<:Real} -# z = zero(a[0]) -# var = TS._params_Taylor1_.var_name -# space = string(" ") -# bigO = TS.bigOnotation[end] ? -# string("+ 𝒪(", var, TS.superscriptify(a.order+1), ")") : -# string("") -# TS._isthinzero(a) && return string(space, z, space, bigO) -# strout::String = space -# ifirst = true -# for i in eachindex(a) -# monom::String = i==0 ? string("") : -# i==1 ? string(" ", var) : string(" ", var, TS.superscriptify(i)) -# @inbounds c = a[i] -# TS._isthinzero(c) && continue -# cadena = TS.numbr2str(c, ifirst) -# strout = string(strout, cadena, monom, space) -# ifirst = false -# end -# strout = strout * bigO -# return strout -# end - -# function TS.pretty_print(a::HomogeneousPolynomial{Interval{T}}) where -# {T<:Real} -# z = zero(a[1]) -# space = string(" ") -# TS._isthinzero(a) && return string(space, z) -# strout::String = TS.homogPol2str(a) -# return strout -# end - -# function TS.pretty_print(a::TaylorN{Interval{T}}) where {T<:Real} -# z = zero(a[0]) -# space = string("") -# bigO::String = TS.bigOnotation[end] ? -# string(" + 𝒪(‖x‖", TS.superscriptify(a.order+1), ")") : -# string("") -# TS._isthinzero(a) && return string(space, z, space, bigO) -# strout::String = space -# ifirst = true -# for ord in eachindex(a) -# pol = a[ord] -# TS._isthinzero(pol) && continue -# cadena::String = TS.homogPol2str( pol ) -# strsgn = (ifirst || ord == 0) ? -# string("") : string(" +") -# strout = string( strout, strsgn, cadena) -# ifirst = false -# end -# strout = strout * bigO -# strout -# end - - -# function TS.homogPol2str(a::HomogeneousPolynomial{Interval{T}}) where -# {T<:Real} -# numVars = get_numvars() -# order = a.order -# z = zero(a.coeffs[1]) -# space = string(" ") -# strout::String = space -# ifirst = true -# iIndices = zeros(Int, numVars) -# for pos = 1:TS.size_table[order+1] -# monom::String = string("") -# @inbounds iIndices[:] = TS.coeff_table[order+1][pos] -# for ivar = 1:numVars -# powivar = iIndices[ivar] -# if powivar == 1 -# monom = string(monom, TS.name_taylorNvar(ivar)) -# elseif powivar > 1 -# monom = string(monom, TS.name_taylorNvar(ivar), -# TS.superscriptify(powivar)) -# end -# end -# @inbounds c = a[pos] -# TS._isthinzero(c) && continue -# cadena = TS.numbr2str(c, ifirst) -# strout = string(strout, cadena, monom, space) -# ifirst = false -# end -# return strout[1:prevind(strout, end)] -# end - - +# Printing-related methods numbr2str function TS.numbr2str(zz::Interval, ifirst::Bool=false) TS._isthinzero(zz) && return string( zz ) plusmin = ifelse( ifirst, string(""), string("+ ") )