(a - K) / (a + K)
,其中 a
是一个正的变量参数,而 K
是一个常数。 在许多情况下,K
是二的幂,这是我的工作相关用例。我正在寻找比直接除法更有效的方法来更准确地计算这个商。可以假设硬件支持融合乘加(FMA),因为此操作由所有主要的 CPU 和 GPU 架构提供,在 C/C++ 中可通过函数 fma()
和 fmaf()
获得。为了便于探索,我正在尝试使用
float
算法。由于我计划将方法移植到 double
算法中,因此不能使用高于参数和结果本机精度的操作。目前为止,我最好的解决方案是: /* Compute q = (a - K) / (a + K) with improved accuracy. Variant 1 */
m = a - K;
p = a + K;
r = 1.0f / p;
q = m * r;
t = fmaf (q, -2.0f*K, m);
e = fmaf (q, -m, t);
q = fmaf (r, e, q);
对于区间
[K/2, 4.23*K]
内的参数a
,上述代码计算出所有输入的商几乎正确舍入(最大误差极接近0.5 ulps),前提是K
为2的幂,并且中间结果没有溢出或下溢。对于非2的幂的K
,此代码仍比基于除法的朴素算法更精确。在性能方面,在浮点倒数可以比浮点除法更快地计算的平台上,此代码可能比朴素方法更快。当
K
=2n时,我做出以下观察:当工作区间的上限增加到8*K
、16*K
等时,最大误差逐渐增加,并开始从下方缓慢逼近朴素计算的最大误差。不幸的是,对于区间的下限却不是这样。如果下限降至0.25*K
,改进后的方法的最大误差等于朴素方法的最大误差。是否有一种方法来计算
q=(a-K)/(a+K)
,以实现与朴素方法和上述代码序列相比,更小的最大误差(用ulp来衡量)并且适用于更广的区间,特别是对于下限小于0.5*K
的区间?效率很重要,但比上述代码中使用的操作稍多一些可能是可以容忍的。
在下面的一个答案中,指出我可以通过将商作为两个操作数的未求值和返回,即作为头尾对
q:qlo
来提高精度,类似于众所周知的双float
和双double
格式。在我的上述代码中,这意味着将最后一行改为qlo=r*e
。这种方法确实有用,我已经考虑过在
pow()
中使用它来扩展精度对数。但它并不能从根本上帮助扩大增强计算提供更准确商的区间。在我正在研究的一个特定情况中,我想使用K=2
(单精度)或K=4
(双精度)来保持主要近似区间较窄,而a
的区间大约为[0,28]。我面临的实际问题是,对于小于0.25*K
的参数,改进后的除法的精度与朴素方法的精度没有太大差别。
(a / (a + k)) - (k / (a + k))
? - Brett Halea
接近K
时会发生减法抵消。 - njuffa