Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
67 changes: 41 additions & 26 deletions src/arithmetic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -250,15 +250,18 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
end

## add! and subst! ##
function ($fc)(v::$T{T}, a::$T{T}, k::Int) where {T<:Number}
if $T == Taylor1
if $T == Taylor1
function ($fc)(v::$T{T}, a::$T{T}, k::Int) where {T<:Number}
@inbounds v[k] = ($f)(a[k])
else
return nothing
end
else
function ($fc)(v::$T{T}, a::$T{T}, k::Int) where {T<:Number}
@inbounds for l in eachindex(v[k])
v[k][l] = ($f)(a[k][l])
end
return nothing
end
return nothing
end

function ($fc)(v::$T{T}, a::T, k::Int) where {T<:Number}
Expand Down Expand Up @@ -507,52 +510,64 @@ end
# Internal multiplication functions
for T in (:Taylor1, :TaylorN)
# NOTE: For $T = TaylorN, `mul!` *accumulates* the result of a * b in c[k]
@eval @inline function mul!(c::$T{T}, a::$T{T}, b::$T{T}, k::Int) where {T<:Number}
@eval @inline @generated function mul!(c::$T{T}, a::$T{T}, b::$T{T}, k::Int) where {T<:Number}
if $T == Taylor1
@inbounds c[k] = a[0] * b[k]
@inbounds for i = 1:k
c[k] += a[i] * b[k-i]
return quote
@inbounds c[k] = a[0] * b[k]
@inbounds for i = 1:k
c[k] += a[i] * b[k-i]
end
end
else
@inbounds mul!(c[k], a[0], b[k])
@inbounds for i = 1:k
mul!(c[k], a[i], b[k-i])
return quote
@inbounds mul!(c[k], a[0], b[k])
@inbounds for i = 1:k
mul!(c[k], a[i], b[k-i])
end
end
end
return nothing
end

@eval @inline function mul_scalar!(c::$T{T}, scalar::NumberNotSeries, a::$T{T}, b::$T{T}, k::Int) where {T<:Number}
@eval @inline @generated function mul_scalar!(c::$T{T}, scalar::NumberNotSeries, a::$T{T}, b::$T{T}, k::Int) where {T<:Number}
if $T == Taylor1
@inbounds c[k] = scalar * a[0] * b[k]
@inbounds for i = 1:k
c[k] += scalar * a[i] * b[k-i]
return quote
@inbounds c[k] = scalar * a[0] * b[k]
@inbounds for i = 1:k
c[k] += scalar * a[i] * b[k-i]
end
end
else
@inbounds mul_scalar!(c[k], scalar, a[0], b[k])
@inbounds for i = 1:k
mul_scalar!(c[k], scalar, a[i], b[k-i])
return quote
@inbounds mul_scalar!(c[k], scalar, a[0], b[k])
@inbounds for i = 1:k
mul_scalar!(c[k], scalar, a[i], b[k-i])
end
end
end
return nothing
end

@eval @inline function mul!(v::$T, a::$T, b::NumberNotSeries, k::Int)
@eval @inline @generated function mul!(v::$T, a::$T, b::NumberNotSeries, k::Int)
if $T == Taylor1
@inbounds v[k] = a[k] * b
return :(@inbounds v[k] = a[k] * b)
else
@inbounds for i in eachindex(v[k])
v[k][i] = a[k][i] * b
return quote
@inbounds for i in eachindex(v[k])
v[k][i] = a[k][i] * b
end
end
end
return nothing
end
@eval @inline function mul!(v::$T, a::NumberNotSeries, b::$T, k::Int)
@eval @inline @generated function mul!(v::$T, a::NumberNotSeries, b::$T, k::Int)
if $T == Taylor1
@inbounds v[k] = a * b[k]
return :(@inbounds v[k] = a * b[k])
else
@inbounds for i in eachindex(v[k])
v[k][i] = a * b[k][i]
return quote
@inbounds for i in eachindex(v[k])
v[k][i] = a * b[k][i]
end
end
end
return nothing
Expand Down
175 changes: 102 additions & 73 deletions src/functions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -318,12 +318,14 @@ end

for T in (:Taylor1, :TaylorN)
@eval begin
@inline function identity!(c::$T{T}, a::$T{T}, k::Int) where {T<:NumberNotSeries}
@inline @generated function identity!(c::$T{T}, a::$T{T}, k::Int) where {T<:NumberNotSeries}
if $T == Taylor1
@inbounds c[k] = identity(a[k])
return :(@inbounds c[k] = identity(a[k]))
else
@inbounds for l in eachindex(c[k])
identity!(c[k], a[k], l)
return quote
@inbounds for l in eachindex(c[k])
identity!(c[k], a[k], l)
end
end
end
return nothing
Expand Down Expand Up @@ -364,69 +366,93 @@ for T in (:Taylor1, :TaylorN)

@inline abs2!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number} = sqr!(c, a, k)

@inline function exp!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
if k == 0
@inbounds c[0] = exp(constant_term(a))
return nothing
end
zero!(c, k)
if $T == Taylor1
if $T == Taylor1
@inline function exp!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
if k == 0
@inbounds c[0] = exp(constant_term(a))
return nothing
end
zero!(c, k)
@inbounds for i = 0:k-1
c[k] += (k-i) * a[k-i] * c[i]
end
@inbounds div!(c, c, k, k)
else
return nothing
end
else
@inline function exp!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
if k == 0
@inbounds c[0] = exp(constant_term(a))
return nothing
end
zero!(c, k)
@inbounds for i = 0:k-1
mul_scalar!(c[k], k-i, a[k-i], c[i])
end
@inbounds div!(c[k], c[k], k)
return nothing
end
return nothing
end

@inline function expm1!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
if k == 0
@inbounds c[0] = expm1(constant_term(a))
return nothing
@inline @generated function expm1!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
ex1 = quote
if k == 0
@inbounds c[0] = expm1(constant_term(a))
return nothing
end
zero!(c, k)
c0 = c[0]+one(c[0])
end
zero!(c, k)
c0 = c[0]+one(c[0])
if $T == Taylor1
@inbounds c[k] = k * a[k] * c0
@inbounds for i = 1:k-1
c[k] += (k-i) * a[k-i] * c[i]
ex2 = quote
@inbounds c[k] = k * a[k] * c0
@inbounds for i = 1:k-1
c[k] += (k-i) * a[k-i] * c[i]
end
@inbounds div!(c, c, k, k)
end
@inbounds div!(c, c, k, k)
else
@inbounds mul_scalar!(c[k], k, a[k], c0)
@inbounds for i = 1:k-1
mul_scalar!(c[k], k-i, a[k-i], c[i])
ex2 = quote
@inbounds mul_scalar!(c[k], k, a[k], c0)
@inbounds for i = 1:k-1
mul_scalar!(c[k], k-i, a[k-i], c[i])
end
@inbounds div!(c[k], c[k], k)
end
@inbounds div!(c[k], c[k], k)
end
return nothing
end

@inline function log!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
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
ex3 = :(return nothing)
return Expr(:block, ex1, ex2, ex3)
end

@inline @generated function log!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
ex1 = quote
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
zero!(c, k)
end
zero!(c, k)
if $T == Taylor1
@inbounds for i = 1:k-1
c[k] += (k-i) * a[i] * c[k-i]
ex2 = quote
@inbounds for i = 1:k-1
c[k] += (k-i) * a[i] * c[k-i]
end
end
else
@inbounds for i = 1:k-1
mul_scalar!(c[k], k-i, a[i], c[k-i])
ex2 = quote
@inbounds for i = 1:k-1
mul_scalar!(c[k], k-i, a[i], c[k-i])
end
end
end
@inbounds c[k] = (a[k] - c[k]/k) / constant_term(a)
return nothing
ex3 = quote
@inbounds c[k] = (a[k] - c[k]/k) / constant_term(a)
return nothing
end
return Expr(:block, ex1, ex2, ex3)
end

@inline function log1p!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
Expand Down Expand Up @@ -456,41 +482,44 @@ for T in (:Taylor1, :TaylorN)
return nothing
end

@inline function sincos!(s::$T{T}, c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
if k == 0
a0 = constant_term(a)
if $T == Taylor1
@inbounds s[0], c[0] = sincos( a0 )
else
@inbounds s[0][1], c[0][1] = sincos( a0 )
end
return nothing
end
zero!(s, k)
zero!(c, k)

if $T == Taylor1
@inbounds for i = 1:k
x = i * a[i]
s[k] += x * c[k-i]
c[k] -= x * s[k-i]
@inline function sincos!(s::$T{T}, c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
if k == 0
a0 = constant_term(a)
@inbounds s[0], c[0] = sincos( a0 )
return nothing
end
zero!(s, k)
zero!(c, k)
@inbounds for i = 1:k
x = i * a[i]
s[k] += x * c[k-i]
c[k] -= x * s[k-i]
end
s[k] = s[k] / k
c[k] = c[k] / k
return nothing
end
else
@inbounds for i = 1:k
mul_scalar!(s[k], i, a[i], c[k-i])
mul_scalar!(c[k], -i, a[i], s[k-i])
@inline function sincos!(s::$T{T}, c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
if k == 0
a0 = constant_term(a)
@inbounds s[0][1], c[0][1] = sincos( a0 )
return nothing
end
zero!(s, k)
zero!(c, k)
@inbounds for i = 1:k
mul_scalar!(s[k], i, a[i], c[k-i])
mul_scalar!(c[k], -i, a[i], s[k-i])
end
@inbounds div!(s[k], s[k], k)
@inbounds div!(c[k], c[k], k)
return nothing
end
end

if $T == Taylor1
s[k] = s[k] / k
c[k] = c[k] / k
else
@inbounds div!(s[k], s[k], k)
@inbounds div!(c[k], c[k], k)
end
return nothing
end

@inline function sincospi!(s::$T{T}, c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
if k == 0
a0 = constant_term(a)
Expand Down
Loading