From 85a4cda784ddc759419090094afef3d2367dea00 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?= Date: Thu, 4 Jul 2024 14:43:57 +0200 Subject: [PATCH] WIP: avoid ifs using generated functions --- src/arithmetic.jl | 67 +++++++++++------- src/functions.jl | 175 +++++++++++++++++++++++++++------------------- src/power.jl | 54 +++++++------- 3 files changed, 172 insertions(+), 124 deletions(-) diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 2af3e8eb..64bb052c 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -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} @@ -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 diff --git a/src/functions.jl b/src/functions.jl index d4fc6001..5c550eaa 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -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 @@ -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} @@ -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) diff --git a/src/power.jl b/src/power.jl index 7e9788d0..a58b2695 100644 --- a/src/power.jl +++ b/src/power.jl @@ -414,38 +414,42 @@ c_k &= 2 \sum_{j=0}^{(k-2)/2} a_{k-j} a_j + (a_{k/2})^2, for T = (:Taylor1, :TaylorN) @eval begin - @inline function sqr!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number} - if k == 0 - sqr_orderzero!(c, a) - return nothing - end + @inline @generated function sqr!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number} + ex1 = quote + if k == 0 + sqr_orderzero!(c, a) + return nothing + end - # Sanity - zero!(c, k) + # Sanity + zero!(c, k) - # Recursion formula - kodd = k%2 - kend = (k - 2 + kodd) >> 1 + # Recursion formula + kodd = k%2 + kend = (k - 2 + kodd) >> 1 + end if $T == Taylor1 - @inbounds for i = 0:kend - c[k] += a[i] * a[k-i] + ex2 = quote + @inbounds for i = 0:kend + c[k] += a[i] * a[k-i] + end + @inbounds c[k] = 2 * c[k] + kodd == 1 && return nothing + @inbounds c[k] += a[k >> 1]^2 + return nothing end - @inbounds c[k] = 2 * c[k] else - @inbounds for i = 0:kend - mul!(c[k], a[i], a[k-i]) + ex2 = quote + @inbounds for i = 0:kend + mul!(c[k], a[i], a[k-i]) + end + @inbounds mul!(c, 2, c, k) + kodd == 1 && return nothing + @inbounds c[k] += a[k >> 1]^2 + return nothing end - @inbounds mul!(c, 2, c, k) end - kodd == 1 && return nothing - - if $T == Taylor1 - @inbounds c[k] += a[k >> 1]^2 - else - accsqr!(c[k], a[k >> 1]) - end - - return nothing + return Expr(:block, ex1, ex2) end end end