Also: Something you could do if you’re happy with a maximum absolute relative error of 1.2e-12 (best of the three ImageMagick approximations; following trick also works with the less accurate/cheaper two) is exploit periodicity, that is:
If x is outside of [-4,4], Use that
sin(pi * x) / (pi * x) = sin(pi * y) / (pi * y) * (y / x) = Sinc(y) * (y/x) where y = x mod 4.
This also tells you how to get derivatives analytically.
Relative error should not be affected by mod and multiplication by y/x.
Whether this is faster than using sin(pi * x) / (pi * x) outside of [-4, 4] is of course an open question. My guess is that it probably is.
Note that this ensures that if x is an integer other than 0, the approximate Sinc is zero, because all my approximations are exactly 0 at -4, -3, -2, -1, 1, 2, 3, 4.
Also: Or course this trick works with any approximation of the normalized Sinc which covers at least [0, 1] (using that Sinc is even and sin is odd and 2 pi periodic).
BETTER:
Suppose that you have a good normalized Sinc approximation over [0, 2]. Call it S. (My approximations are good over [-4, 4] so [0,2] is covered.)
Set y = abs(x)
if y <= 2 (or 4, with my approximation):
return S(y)
else:
z = y * 0.5 (or y * 0.25 with my approximation)
z = 2. * (z - floor(z)) (or 4. * (z - floor(z)) with my approximation; this may help with derivative discontinuities; this can be checked.)
return S(z) * (z/y)
Use the floor function (or whatever gives you z -floor(x)) which is fastest on non-negative numbers of the relevant floating point type.
Note you can make this branchless with conditional move, e.g. let z = y and “factor” = 1. in the first branch so that you only have one call to S.
EVEN BETTER BUT NOT DONE YET:
With a tiny bit more work you can make this work with a good approximation on [0,1] (by reusing the * (z/y) trick on [1,2] for example). I could produce a superior even polynomial (even because you want the derivative to be exactly 0 at 0) or Padé-ish approximant on that interval (or on [0,2]; obviously the shorter the interval the more accurate one can make it with the same number of flops) but unfortunately I’d need to be paid to do the work.
P.S. Different context, but this illustrates how little difference there is between using ImageMagick’s SincFast and the exact Sinc. repeatedly rotating a thin color line with ImageMagick
P.S. You may find Section 20 and Appendix C interesting Chantal Racette’s Master’s Thesis
Both of us are more than happy to give permission to use our code modifications with suitable attribution. It should be very easy to modify it to get a high accuracy Sinc approximation on [0,1] or [0,2] (in relative error if you prefer, and you should prefer). I can provide guidance. I believe I did further work and improved this. I’d have to dig (long time ago).
1 Like