c - 高效计算 (a - K)/(a + K) 并提高准确性

Question

在各种情况下，例如对于数学函数的参数约简，需要计算 (a - K)/(a + K)，其中 a 是一个正变量参数，K 是一个常量。在许多情况下，K 是 2 的幂，这是与我的工作相关的用例。我正在寻找比直接除法更准确地计算这个商的有效方法。可以假定对融合乘加 (FMA) 的硬件支持，因为目前所有主要 CPU 和 GPU 架构都提供此操作，并且可通过函数 fma() 在 C/C++ 中使用和 fmaf()。

为了便于探索，我正在试验 float 算术。由于我还计划将该方法移植到 double 算术，因此不得使用比参数和结果的 native 精度更高的操作。到目前为止我最好的解决方案是:

 /* 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 的幂，并且中间结果没有上溢或下溢。对于 K 不是 2 的幂，此代码仍然比基于除法的朴素算法更准确。就性能而言，此代码可以更快比平台上的朴素方法更快，在平台上计算浮点倒数比浮点除法更快。

当 K = 2ⁿ 时，我做了以下观察:当工作区间的上限增加到 8*K 时，16*K, ... 最大误差逐渐增加，并开始从下方慢慢逼近朴素计算的最大误差。不幸的是，区间的下限似乎并非如此。如果下界下降到0.25*K，则上述改进方法的最大误差等于朴素方法的最大误差。

有没有一种计算 q = (a - K)/(a + K) 的方法可以实现比朴素方法更小的最大误差(以 ulp 衡量与数学结果相比)和上面的代码序列，在更宽的区间内，特别是对于下限小于 0.5*K 的区间？ 效率很重要，但比实际操作多一些在上面的代码中使用可能是可以容忍的。

---

在下面的一个回答中，有人指出我可以通过将商返回为两个操作数的未计算总和来提高准确性，即作为头尾对 q:qlo，即类似于众所周知的 double float 和 double double 格式。在我上面的代码中，这意味着将最后一行更改为 qlo = r * e。

这种方法当然很有用，我已经考虑过将其用于 pow() 中的扩展精度对数。但它并不能从根本上帮助扩大增强计算提供更准确商数的区间。在我正在查看的特定情况下，我想使用 K=2(对于单精度)或 K=4(对于 double )来保持主要近似值区间变窄，a 的区间大致为 [0,28]。我面临的实际问题是，对于 < 0.25*K 的参数，改进除法的准确性并不比使用朴素方法好很多。

Answer 1

如果 a 比 K 大，则 (a-K)/(a+K) = 1 - 2K/(a + K) 将给出一个很好的近似值。如果 a 与 K 相比较小，则 2a/(a + K) - 1 将给出一个很好的近似值。如果 K/2 ≤ a ≤ 2K，则 a-K 是精确运算，因此进行除法会得到不错的结果。