如何在GLSL中无分支地实现sinc(x)？

Question

我想在GLSL中实现

sinc(x) = sin(x) / x

函数。这是我的第一次尝试：

float sinc(float x) {
    if (abs(x) < 1.0e-5) {
        return 1.0;
    } else {
        return sin(x) / x;
    }
}

看起来不错，但是我想避免分支，所以我将其更改为：

float sinc(float x) {
    return mix(1.0, sin(x) / x, step(1.0e-5, abs(x)));
}

但是这个表达式有一个可能被零除的表达式。是否可以在保持代码分支自由的同时避免这种情况？

Answer 1

是的。您的第一次尝试并不遥远。

顺便说一句，测试与零相等的速度稍微快一些，并且对于真正的零也可以正常工作，这是唯一会导致计算除法麻烦的值。

Sinc 是一个偶函数，因此将其定义为 x 的正值并避免恰好为零就足够了。对于任何小的非零 x 都表现良好。

尽管您可能想检查一下这样做实际上是否更快。

float sinc(float x) {
   x = abs(x)+1e-14;   // add a tiny amount to avoid divide by zero
   return sin(x) / x;
}

以浮点精度执行您所要求的操作。 您还需要在反汇编中检查 abs(x) 的实现是否隐藏条件分支。

添加一个小常数的技巧是安全的，因为对于任何合理的 x 值，10^-14 都太小而无法改变答案。对于非常小的 x，我们可以查看 sin(x)/x

的多项式展开式

sinc(x) = 1 - x^2/6

一次 x^2/6 < machine_eps then sinc(x)==1 or sin(x) == x if you prefer.

对于单精度，该下限约为 x= 0.002，对于双精度，该下限约为 2e-8。 10^-30 是双精度修正的合理选择。

具体取决于您使用它的用途和 x 的范围，查找表可能会更快。计算函数 sin(x) 确实非常慢，而且除法也有很高的延迟。如果您可以在这里测量条件分支的效果（特别是仅当 x==0.0 时才采用），我会感到惊讶。根据您想要做什么，可以使用更好的插值函数......

Answer 2

也许 ImageMagick“永远”使用的 SincFast 近似值会有用？ ImageMagick resize.c 源代码。例如：

const double xx = x*x;
/*
  Max. abs. rel. error 1.2e-12 < 1/2^39.
*/
const double c0 = 0.173611111110910715186413700076827593074e-2L;
const double c1 = -0.289105544717893415815859968653611245425e-3L;
const double c2 = 0.206952161241815727624413291940849294025e-4L;
const double c3 = -0.834446180169727178193268528095341741698e-6L;
const double c4 = 0.207010104171026718629622453275917944941e-7L;
const double c5 = -0.319724784938507108101517564300855542655e-9L;
const double c6 = 0.288101675249103266147006509214934493930e-11L;
const double c7 = -0.118218971804934245819960233886876537953e-13L;
const double p =
  c0+xx*(c1+xx*(c2+xx*(c3+xx*(c4+xx*(c5+xx*(c6+xx*c7))))));
const double d0 = 1.0L;
const double d1 = 0.547981619622284827495856984100563583948e-1L;
const double d2 = 0.134226268835357312626304688047086921806e-2L;
const double d3 = 0.178994697503371051002463656833597608689e-4L;
const double d4 = 0.114633394140438168641246022557689759090e-6L;
const double q = d0+xx*(d1+xx*(d2+xx*(d3+xx*d4)));
return((xx-1.0)*(xx-4.0)*(xx-9.0)*(xx-16.0)/q*p);