有效地计算(a-K)/(a + K)并提高精度

Question

在各种情况下,例如对于数学函数的参数减少,需要计算(a - K) / (a + K),其中a是正变量参数并且K是常数.在许多情况下,K是2的幂,这是与我的工作相关的用例.我正在寻找比直接划分更准确地计算这个商的有效方法.可以假设对融合乘法 - 加法(FMA)的硬件支持,因为此操作由所有主要的CPU和GPU架构提供,并且可以通过函数fma()和C/C++获得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);

对于a区间中的参数[K/2, 4.23*K],上面的代码计算所有输入几乎正确舍入的商(最大误差非常接近0.5 ulps),前提K是功率为2,并且中间结果中没有溢出或下溢.对于K不是2的幂,该代码仍然比基于除法的朴素算法更准确.在性能方面,这个代码可以比平台上的朴素方法更快,在这些平台上浮点倒数可以比浮点除法更快地计算.

我提出以下意见时K= 2 ^ñ:当绑定的工作间隔增加到上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的论证,改进除法的准确性并不比使用朴素方法好.

Answer 1

由于我的目标只是扩大获得准确结果的区间，而不是找到适用于 的所有可能值的解决方案a，因此对所有中间计算使用双重float算术似乎成本太高。

进一步思考这个问题，很明显，e在我的问题的代码中，除法余数的计算是获得更准确结果的关键部分。数学上，余数为 (aK) - q * (a+K)。在我的代码中，我简单地将m(aK) 表示为 ，并将 (a+k) 表示为m + 2*K，因为这在数值上比直接表示提供了更好的结果。

通过相对较小的额外计算成本， (a+K) 可以表示为 double- float，即头尾对p:plo，这导致我的原始代码的修改版本如下：

/* Compute q = (a - K) / (a + K) with improved accuracy. Variant 2 */
m = a - K;
p = a + K;
r = 1.0f / p;
q = m * r;
mx = fmaxf (a, K);
mn = fminf (a, K);
plo = (mx - p) + mn;
t = fmaf (q, -p, m);
e = fmaf (q, -plo, t);
q = fmaf (r, e, q);

测试表明，这可以a在 [K/2, 2 ²⁴ *K) 范围内提供几乎正确的舍入结果，从而可以大幅增加实现准确结果的区间上限。

扩大下端的区间需要更准确地表示 (aK)。我们可以将其计算为双头float尾对m:mlo，这导致以下代码变体：

/* Compute q = (a - K) / (a + K) with improved accuracy. Variant 3 */
m = a - K;
p = a + K;
r = 1.0f / p;
q = m * r;
plo = (a < K) ? ((K - p) + a) : ((a - p) + K);
mlo = (a < K) ? (a - (K + m)) : ((a - m) - K);
t = fmaf (q, -p, m);
e = fmaf (q, -plo, t);
e = e + mlo;
q = fmaf (r, e, q);

详尽的测试表明这如何a在区间 [K/2 ²⁴ , K*2 ²⁴ ) 内提供几乎正确的舍入结果。不幸的是，与我的问题中的代码相比，这需要十次额外的操作，这是为了将最大误差从 1.625 ulp 左右（通过简单计算降低到接近 0.5 ulp）而付出的高昂代价。

p正如我在问题中的原始代码中一样，可以用 (aK) 来表达 (a+K)，从而消除,尾部的计算plo。这种方法产生以下代码：

/* Compute q = (a - K) / (a + K) with improved accuracy. Variant 4 */
m = a - K;
p = a + K;
r = 1.0f / p;
q = m * r;
mlo = (a < K) ? (a - (K + m)) : ((a - m) - K);
t = fmaf (q, -2.0f*K, m);
t = fmaf (q, -m, t);
e = fmaf (q - 1.0f, -mlo, t);
q = fmaf (r, e, q);

如果主要焦点是减少间隔的下限，那么这将是有利的，这是我在问题中所解释的特别关注点。对单精度情况的详尽测试表明，当 K=2 时，对于区间 [K/2 ²⁴ , 4.23*K] 内的值，会生成^{n 个}几乎正确舍入的结果。总共有 14 或 15 个操作（取决于架构是否支持完整预测或仅支持条件移动），这比我的原始代码多需要七到八个操作。a

最后，可以直接基于原始变量进行残差计算，a以避免计算m和 时固有的误差p。这导致以下代码，对于 K = 2 ⁿ，计算a区间 [K/2 ²⁴ , K/3) 中几乎正确的舍入结果：

/* Compute q = (a - K) / (a + K) with improved accuracy. Variant 5 */
m = a - K;
p = a + K;
r = 1.0f / p;       
q = m * r;
t = fmaf (q + 1.0f, -K, a);
e = fmaf (q, -a, t);
q = fmaf (r, e, q);