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

Question

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

为了便于探索，我正在试验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，如果 K 是 2 的幂，则上面的代码计算所有输入的商几乎正确四舍五入（最大误差非常接近 0.5 ulps），并且在中间结果中没有上溢或下溢。对于K 不是 2 的幂，这段代码仍然比基于除法的朴素算法更准确。就性能而言，此代码可以更快在平台上计算浮点倒数比浮点除法更快的幼稚方法。

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

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

---

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

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

Answer 1

我真的没有答案（正确的浮点错误分析非常乏味），但有一些观察：

• 快速倒数指令（例如 RCPSS）不如除法准确，因此如果使用这些指令，您可能会发现准确度有所下降。
• m 精确计算，如果 ∈ [0.5×K_b, 2¹⁺ⁿ×K_b)，其中 K_b 是2 低于 K（如果 K 是 2 的幂，则为 K 本身），n 是 K 的有效数字中尾随零的数量（即，如果 K 是 2 的幂，则 n=23）。
• 这类似于来自Dekker (1971) 的div2 算法的简化形式：要扩大范围（尤其是下限），您可能必须从这里合并更多的校正项（即存储m作为 2 个floats 的总和，或使用double）。

Answer 2

如果您可以放宽 API 以返回另一个对错误建模的变量，那么解决方案会变得更加简单：

float foo(float a, float k, float *res)
{
    float ret=(a-k)/(a+k);
    *res = fmaf(-ret,a+k,a-k)/(a+k);
    return ret;
}

该方案只处理除法的截断错误，不处理a+k和a-k的精度损失。

为了处理这些错误，我想我需要使用双精度，或者使用 bithack 来使用定点。

更新测试代码以人工生成非零最低有效位 在输入中

测试代码

https://ideone.com/bHxAg8

Answer 3

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

Answer 4

一种可能性是使用经典的 Dekker/Schewchuk 将 m 和 p 的误差跟踪到 m1 和 p1：

m=a-k;
k0=a-m;
a0=k0+m;
k1=k0-k;
a1=a-a0;
m1=a1+k1;

p=a+k;
k0=p-a;
a0=p-k0;
k1=k-k0;
a1=a-a0;
p1=a1+k1;

然后，纠正幼稚的划分：

q=m/p;
r0=fmaf(p,-q,m);
r1=fmaf(p1,-q,m1);
r=r0+r1;
q1=r/p;
q=q+q1;

这将花费你 2 个师，但如果我没有搞砸的话，应该会接近一半。

但是这些除法可以用 p 的倒数乘法代替，没有任何问题，因为第一个不正确舍入的除法将由余数 r 补偿，第二个不正确舍入的除法并不重要（校正 q1 的最后一位不会' t改变任何东西）。

Answer 5

问题在于(a + K) 中的添加。 (a + K) 中的任何精度损失都会被除法放大。问题不在于部门本身。

如果a 和K 的指数相同（几乎）没有精度损失，并且如果指数之间的绝对差大于有效数字大小，则(a + K) == a（如果a 具有更大的震级）或(a + K) == K（如果K有更大的震级）。

没有办法阻止这种情况。增加有效数字大小（例如，在 80x86 上使用 80 位“扩展双精度”）仅有助于稍微扩大“准确结果范围”。要了解原因，请考虑 smallest + largest（其中 smallest 是 32 位浮点数可以是最小的正非正规数）。在这种情况下（对于 32 位浮点数），您需要大约 260 位的有效位大小才能完全避免精度损失。做（例如）temp = 1/(a + K); result = a * temp - K / temp; 也无济于事，因为您仍然遇到完全相同的(a + K) 问题（但它可以避免(a - K) 中的类似问题）。你也不能这样做result = anything / p + anything_error/p_error，因为除法不是那样工作的。

对于可以适合 32 位浮点的 a 的所有可能正值，我只能想到 3 种替代方案来接近 0.5 ulps。没有一个可能是可以接受的。

第一个替代方案涉及为a 的每个值预先计算一个查找表（使用“大实数”数学），对于 32 位浮点（和对于 64 位浮点完全疯狂）。当然，如果a 的可能值范围小于“任何可以放入 32 位浮点数的正值”，则查找表的大小将会减小。

第二种选择是在运行时使用其他东西（“大实数”）进行计算（并转换为/从 32 位浮点数）。

第三种选择涉及“某物”（我不知道它叫什么，但它很贵）。将舍入模式设置为“舍入到正无穷大”并计算temp1 = (a + K); if(a < K) temp2 = (a - K);，然后切换到“舍入到负无穷大”并计算if(a >= K) temp2 = (a - K); lower_bound = temp2 / temp1;。接下来执行a_lower = a 并尽可能减少a_lower 并重复“lower_bound”计算，并继续执行此操作，直到获得lower_bound 的不同值，然后恢复为a_lower 的先前值。之后，您执行基本相同的操作（但舍入模式相反，并且递增而不是递减）以确定upper_bound 和a_upper（从a 的原始值开始）。最后，插值，如a_range = a_upper - a_lower; result = upper_bound * (a_upper - a) / a_range + lower_bound * (a - a_lower) / a_range;。请注意，您将需要计算初始上限和下限，如果它们相等，则跳过所有这些。另请注意，这都是“理论上的，完全未经测试”，我可能在某个地方搞砸了。

我的主要意思是（在我看来）你应该放弃并接受你无法做任何事情来接近 0.5 ulp。对不起.. :)

Answer 6

由于我的目标只是扩大获得准确结果的区间，而不是找到适用于所有可能的 a 值的解决方案，因此似乎对所有中间计算都使用 double-float 算法太贵了。

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

以相对较小的额外计算成本，(a+K) 可以表示为双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);

测试表明，这可以在 [K/2, 2²⁴*K) 中为 a 提供几乎正确的四舍五入结果，从而可以显着增加准确的区间上限取得了成果。

在下端加宽区间需要更准确的 (a-K) 表示。我们可以将其计算为双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);

详尽的测试表明这如何在区间 [K/2²⁴, K*2²⁴) 中为 a 提供几乎正确的舍入结果。不幸的是，与我的问题中的代码相比，这需要额外执行 10 次操作，这是为了将最大误差从大约 1.625 ulp 降低到接近 0.5 ulp 而付出的巨大代价。

正如我在问题中的原始代码一样，可以用 (a-K) 来表示 (a+K)，从而消除了 p、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²⁴ 中的 a 的值，会产生几乎正确的舍入结果， 4.23*K]。总共有 14 或 15 个操作（取决于架构是支持完全预测还是仅支持条件移动），这需要比我的原始代码多 7 到 8 个操作。

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

/* 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);