我的主要目标是最小化ulp错误,次要目标是最小化错误舍入结果的数量,在改进代码最多比下面发布的版本最低的约束下.任何提高准确度的改进,比如0.2 ulp,都是值得欢迎的.添加几个FMA(融合乘法 – 加法)会很好,另一方面我希望有人能够找到一个采用快速rsqrtf()(倒数平方根)的解决方案.
生成的C99代码应该适用于矢量化,可能是通过一些简单的直接转换.所有中间计算必须以函数参数和结果的精度发生,因为任何切换到更高精度可能会产生严重的负面性能影响.代码必须在IEEE-754非正常支持和FTZ(刷新到零)模式下正常工作.
到目前为止,我已经确定了以下两个候选实现.请注意,可以通过单次调用log1pf()轻松地将代码转换为无分支可向量化版本,但在此阶段我还没有这样做,以避免不必要的混淆.
/* for a >= 0,asinh(a) = log (a + sqrt (a*a+1))
= log1p (a + (sqrt (a*a+1) - 1))
= log1p (a + sqrt1pm1 (a*a))
= log1p (a + (a*a / (1 + sqrt(a*a + 1))))
= log1p (a + a * (a / (1 + sqrt(a*a + 1))))
= log1p (fma (a / (1 + sqrt(a*a + 1)),a,a)
= log1p (fma (1 / (1/a + sqrt(1/a*a + 1)),a)
*/
float my_asinhf (float a)
{
float fa,t;
fa = fabsf (a);
#if !USE_RECIPROCAL
if (fa >= 0x1.0p64f) { // prevent overflow in intermediate computation
t = log1pf (fa) + 0x1.62e430p-1f; // log(2)
} else {
t = fmaf (fa / (1.0f + sqrtf (fmaf (fa,fa,1.0f))),fa);
t = log1pf (t);
}
#else // USE_RECIPROCAL
if (fa > 0x1.0p126f) { // prevent underflow in intermediate computation
t = log1pf (fa) + 0x1.62e430p-1f; // log(2)
} else {
t = 1.0f / fa;
t = fmaf (1.0f / (t + sqrtf (fmaf (t,t,fa);
t = log1pf (t);
}
#endif // USE_RECIPROCAL
return copysignf (t,a); // restore sign
}
使用特定的log1pf()实现,该实现精确到< 0.6 ulps,我在所有232个可能的IEEE-754单精度输入上进行详尽测试时,我观察到以下错误统计数据.当USE_RECIPROCAL = 0时,最大错误为1.49486 ulp,并且有353,587,822个错误的舍入结果. USE_RECIPROCAL = 1时,最大误差为1.50805 ulp,并且只有77,569,390个错误的舍入结果. 在性能方面,如果倒数和完全除法占用大致相同的时间量,则变体USE_RECIPROCAL = 0将更快,但如果可获得非常快的互惠支持,则变体USE_RECIPROCAL = 1可能更快. 答案可以假设所有基本算术,包括FMA(融合乘法 – 加法)都是根据IEEE-754舍入到最接近或偶数模式正确舍入的.此外,可以使用更快,几乎正确舍入的版本的倒数和rsqrtf(),其中“几乎正确舍入”意味着最大的ulp误差将限制为0.53 ulps以及绝大多数结果,例如> 95％,正确圆润.可以使用定向舍入的基本算术,而无需额外的性能成本.
