I would like to implement the sinc(x) = sin(x) / x function in GLSL. This is my first attempt:
float sinc(float x) {
if (abs(x) < 1.0e-5) {
return 1.0;
} else {
return sin(x) / x;
}
}
It seems well, but I would like to avoid branches, so I changed it to:
float sinc(float x) {
return mix(1.0, sin(x) / x, step(1.0e-5, abs(x)));
}
But this expression has a possible division by zero expression. Is it possible to avoid that while keeping the code branch free?
Yes. Your first attempt isn't that far off.
BTW Testing for equality with zero is marginally faster and works OK for true zero which is the only value that causes trouble for computing the division.
Sinc is an even function so that defining it for positive only values of x and avoiding exactly zero is sufficient. It is well behaved for any small non-zero x.
Although you might want to check that it is actually faster done this way.
float sinc(float x) {
x = abs(x)+1e-14; // add a tiny amount to avoid divide by zero
return sin(x) / x;
}
Does what you have asked for at float precision. You will also need to check in the disassembly that the implementation of abs(x) isn't hiding a conditional branch.
This adding a small constant trick is safe because for any reasonable value of x, 10^-14 is too small to alter the answer. And for very small x we can look at the polynomial expansion of sin(x)/x
sinc(x) = 1 - x^2/6
Once x^2/6 < machine_eps then sinc(x)==1 or sin(x) == x if you prefer.
For single precision this lower limit is at about x= 0.002 and for double precision it is around 2e-8. 10^-30 is a reasonable choice of fixup for double precision.
Depending on exactly what you are using it for and the range of x a lookup table might be faster still. Evaluating the function sin(x) is really quite slow and divisions have high latency too. I'd be surprised if you could measure the effect of the conditional branch here (especially if it is only taken when x==0.0). Better interpolating functions are available depending on what it is you want to do...