I think the point is that if the power is a literal constant, as opposed to a variable, the compiler can simplify the calculation a lot. (It should suffice to write 1/7 rather than 0.14285714285714285, because the compiler can evaluate the former to the latter.)
However, I expect you could do even better by writing a special function to compute the seventh root of 1+x, similar to how we have special functions for cbrt and sqrt. Look at the code in cbrt.jl for the implementation of cbrt, for inspiration.
For example, when you compute (1+x)^{1/7}, is x always relatively small, e.g. |x| < 0.1? If so, a very simple approach is a Taylor expansion:
# compute (1+x)^(1/7) approximately, in double precision, by a
# Taylor series for small x
taylor7th(x) = evalpoly(x, (1, 1/7, -3/49, 13/343, -65/2401, 351/16807, -1989/117649, 81549/5764801, -489294/40353607, 2990130/282475249))
This matches (1+x)^(1/7) to machine precision for x=0.01, to 12 digits for x=0.1, and to 5 digits for x=0.5. However, it is about 5x faster.
If you need something accurate for larger x, look at cbrt.jl … basically you compute an approximation and then polish it with Newton steps, which converge very quickly.
@stevengj unfortunately I am not that good at this kind of bit manipulation etc. Is there a ressourcce which explains and goes through for example this cube root function?`
Here’s a version of the optimized algorithm for x^(1/7) It’s not perfect (I think I could make it a bit faster with some extra work), but it should work for all (finite) x and is accurate to a couple floating point epsilon:
function fancy7th(x::Float64)
# todo tune the magic constant
# initial guess based on fast inverse sqrt trick but adjusted to compute x^(1/8)
t = copysign(reinterpret(Float64, 0x37f306fe0a31b715 + reinterpret(UInt64,abs(x))>>3), x)
@fastmath for _ in 1:3
# newton's method for t^3 - x/t^4 = 0
t2 = t*t
t3 = t2*t
t4 = t2*t2
xot4 = x/t4
t = t - t*(t3 - xot4)/(4*t3 + 3*xot4)
end
t
end
It is a little slower than @stevengj’s solution, but it is more accurate and works over the full range.
Here is an (IMO better version). This one gets within 5*10^-11 over the range and is 5.6ns (it saves 1 iteration by making a slightly better initial guess)
function fancy7th(x::Float64)
# todo tune the magic constant
# initial guess based on fast inverse sqrt trick but adjusted to compute x^(1/7)
t = copysign(reinterpret(Float64, 0x36cd000000000000 + reinterpret(UInt64,abs(x))÷7), x)
@fastmath for _ in 1:2
# newton's method for t^3 - x/t^4 = 0
t2 = t*t
t3 = t2*t
t4 = t2*t2
xot4 = x/t4
t = t - t*(t3 - xot4)/(4*t3 + 3*xot4)
end
t
end
I doubt what I’m about to suggest is faster than the other things being suggested, but the following is definitely more accurate for small values of P * Cb. The original writing can suffer from catastrophic cancellation.
expm1(γ⁻¹ * log1p(P * Cb))
If you do something like Taylor series or like the fancy7th suggested above, it definitely makes sense to pull the ±1s you have into those (which would have the same accuracy benefit as expm1/log1p provide).
Thanks @mikmoore ! I’ve made sure to implement this. I’ve never seen these functions before in Julia haha, but from the ? I could see it is quite clever to e^(x-1) and then log1p since it all reduces to that expression you showed above.
It did not make code faster, perhaps a little bit slower - but good to be more precise, especially around the small numbers.
