mabruzzo
commented
2 years ago Here's the code for the branchless exponential function:
# in the file where you define your branchless exponential function,
# you want to make sure that the machine is little endian. The following
# code does that check for you...
const LITTLE_ENDIAN = (ENDIAN_BOM == 0x04030201)
@assert LITTLE_ENDIAN
# polynomial coefficients in correction were taken from
# https://gitlab.mpcdf.mpg.de/bbramas/inastemp/-/blob/master/Src/Common/InaFastExp.hpp
# (That above repository has an MIT license that needs to be noted somewhere if this is used,
# unless you derive them yourself)
# you can adjust the order of the polynomials to adjust accuracy (see the above link for
# more options).
#
# Some random spot-checks seemed to indicate that the following polynomial is accurate to
# better than 1e-5. But my checks were NOT rigorous
@inline function correction(x::F) where {F<:AbstractFloat}
coef_1, coef_2, coef_3, coef_4, coef_5 = (-3.70138142771437266806e-06,
3.07033820309224325662e-01,
-2.41638288055762540107e-01,
-5.16904731562965388814e-02,
-1.36976563343097993558e-02)
(((((coef_5 * x + coef_4) * x + coef_3) * x) + coef_2) * x) + coef_1
end
# this whole approach is based on
# https://www.researchgate.net/publication/272178514_Fast_Exponential_Computation_on_SIMD_Architectures
const COEFF_LOG2E = log2(ℯ)
const COEFF_A = 2^52
const COEFF_B = 2^52*1023
@inline my_ifelse(cond, a, b) = cond*a + (1-cond) * b
@inline function branchless_exp(x::Float64)
# inrange and my_ifelse are present purely for error handling.
# When x^2 > 709.7, they make the function returns -1
# Note: This error handling makes the function ~50% slower
# (but my benchmarks still suggest that it's ~3 times faster
# than exp)
inrange = (x^2 < 709.7^2)
x = my_ifelse(inrange, x, 0) * COEFF_LOG2E
fractional_part::Float64 = x - floor(x)
x -= correction(fractional_part)
casted_integer::Int64 = unsafe_trunc(Int64, COEFF_A * x + COEFF_B)
my_ifelse(inrange, reinterpret(Float64, casted_integer),-1.0)
end
# keep in mind: to use this function with ForwardDiff, you'll have to
# define a specialization of `branchless_exp` that takes a DualNumber.
# If you look at dual.jl, within ForwardDiff.jl, they do this for
# functions like `sincospi`When you broadcast over a call to this function, Julia has a really easy time auto-vectorizing. It would take a little effort to write a version of this function for 32-bit floats
ajwheeler
As promised, I'm posting my faster exponential function. I'm curious to hear if using this inside of
voigtbuys you any speed.Tonight, I also put together a branchless version of
voigt. Unfortunately, all of the extra work that is being done by the branchless voigt makes it slower than the unvectorized version. I suspect that if someone spent a lot of time refactoring, they might be able to get it to be slightly faster. But I don't think that's worth the time. As you had suggested a while back, I think a better use of time would be to replace the current approximation with a better approach.I did some googling and the consensus seems to be that we should leverage the relationship of the Voigt function with the real function. According to eqn 9 of Thompson (93),
V(u,a) = Re[exp(z²) erfc(z)], wherez = a + i*u. I suspect that nobody will have time to pursue this for a while, so I'll briefly document what I read. I've seen 2 ways to do this (there are probably more):I think that Turbospectrum's implementation is based on the same underlying approach as Kuntz (1997). (Interestingly, I think that turbospectrum and the article might be suggesting that for
V(a,v),z = a - iv...). Under this approach, I think that they directly approximate the whole expressionRe[exp(z²) erfc(z)], but it introduces at least some level of branching. I think this article talks a little about variants of this approach.At a glance, I suspect this might potentially be slightly faster to adopt a strategy like the one used in this article (if we target a particular accuracy). They seemed to use an approximation for
exp(z²) erfc(z) ≈ (Σⱼᵖ aⱼ zʲ)/(z⁽ᵖ⁺¹⁾ + Σⱼᵖ bⱼ zʲ)with coefficients from here, and directly computedexp(z²). The benefit of this approach is that there is basically no-branching. However, pluggingRe(exp(-1*(x + i*y)^2) * (a + i b))into wolfram alpha suggests that you also need to evaluate sin and cos (which I think boils down to 2 more polynomial evaluations)...If either of these approaches are adopted, I believe some care is required to avoid using Julia's complex datatype (since I think ForwardDiff.jl is incompatible with complex numbers)