julia
julia copied to clipboard
JuliaLang
Faster and more accurate sinpi, cospi, sincospi for Float32, Float64
Instead of using sin_kernel and cos_kernel, use better kernels that need less precision.
I benchmark this as 14% faster for Float64, 7% faster for Float32.
Accuracy for Float64 is .64 ULPs compared to .75 for master.
Accuracy for Float32 is .5 ULP
Thanks Oscar. For my own curiosity, how did you generate these polynomial coefficients? If it's not too much code, it would be cool to record that for future usage with other special functions, so that you're not the only person who knows how to do this.
Test failures look like cospi_kernel(::BigFloat) is trying to be called.
The coefficient generation is slightly complex, but the TLDR is to use Remez to generate a minimax polynomial over the range of inputs, specialized to the even/oddness of the polynomial. I really need to put the general principals all into a blogpost at some point, but the TLDR for this is
using Remez
sinpi_coeffs = [(big(pi))^n/factorial(big(n)) for n in 1:2:40] #nonzero taylor coefficients in bigfloat
sinpi_kinda(x) = evalpoly(x, sinpi_coeffs)
bounds = (0, 0.25^2) #squared because we are getting coefficients for sqrt(sinpi(x))
n = 7 # degree 7 polynomial
exact_coeficients = ratfn_minimax(sinpi_kinda, bounds, n, 0)[1]
#now do lots of messing around with rounding to get things exact.
This has been surprisingly tricky to get right. At the moment, the only failing tests is the behavior as to whether 0.0 or -0.0 is returned. To me, neither is clearly right. Is there a good reason why the specific signs need to be tested? If not, I will probably loosen the tests slightly to check for ==0 rather than ===0.0
Is it worth trying to get this done? I suppose too late for 1.8 but maybe 1.9?
re: signed zeros, the IEEE754-2019 spec states
sinPi(+n) is +0 and sinPi(−n) is −0 for positive integers n cosPi(n + 1⁄2) is +0 for any integer n when n + 1⁄2 is representable
