高效计算 (a - K) / (a + K) 并提高准确性
Efficiently computing (a - K) / (a + K) with improved accuracy
在各种情况下，例如对于数学函数的参数约简，需要计算 (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的幂，中间结果没有上溢或下溢。对于 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 的参数，改进除法的准确性并不比使用朴素方法好得多。
我真的没有答案（适当的浮点错误分析非常乏味）但有几点观察：
- 快速倒数指令（例如 RCPSS）不如除法准确，因此如果使用它们，您可能会发现精度会降低。
m 如果一个 ∈ [0.5×Kb, 21+n×Kb), 其中Kb 是低于 K 的 2 的幂（如果 K 是 2 的幂，则为 K 本身），n 是 K 的尾随零的数量（即，如果 K 是 2 的幂, 那么 n=23).- 这类似于 Dekker (1971) 中
div2 算法的简化形式：要扩大范围（尤其是下限），您可能需要从中合并更多校正项（即将 m 存储为 2 float 的总和，或使用 double）。
如果您可以将 API 放宽为 return 另一个对错误建模的变量，那么解决方案就会变得简单得多：
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 来使用定点。
测试代码更新为人工生成非零最低有效位
在输入
测试代码
如果 a 比 K 大，则 (a-K)/(a+K) = 1 - 2K / (a + K) 将给出一个很好的近似值。如果 a 与 K 相比较小，则 2a / (a + K) - 1 将给出一个很好的近似值。如果 K/2 ≤ a ≤ 2K，那么 a-K 是精确运算，所以做除法会得到一个不错的结果。
一种可能性是用经典的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 个师，但如果我没有搞砸的话应该接近一半 ulp。
但是这些除法可以用 p 的倒数乘法代替，没有任何问题，因为第一个不正确的舍入除法将由余数 r 补偿，而第二个不正确的舍入除法并不重要（修正 q1 的最后一位不会改变任何东西）。
问题是 (a + K) 中的加法。 (a + K) 中的任何精度损失都会被除法放大。问题不在于部门本身。
如果 a 和 K 的指数相同（几乎），则不会丢失精度，并且如果指数之间的绝对差大于尾数大小，则 (a + K) == a（如果a有更大的量级）或(a + K) == K（如果K有更大的量级）。
没有办法阻止这种情况。增加有效数字大小（例如，在 80x86 上使用 80 位 "extended double"）只会稍微扩大 "accurate result range"。要理解原因，请考虑 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 的备选方案。 None 很可能会被接受table.
第一个替代方案涉及为 a 的每个值预计算查找 table（使用 "big real number" 数学），最终（使用一些技巧）约为 2 GiB对于 32 位浮点（对于 64 位浮点完全疯狂）。当然，如果 a 的可能值范围小于 "any positive value that can fit in a 32-bit float"，则查找的大小 table 将减少。
第二种选择是在运行时使用其他东西（"big real number"）进行计算（并转换to/from 32位浮点数）。
第三种选择涉及，"something"（我不知道它叫什么，但它很贵）。将舍入模式设置为"round to positive infinity"并计算temp1 = (a + K); if(a < K) temp2 = (a - K);然后切换到"round to negative infinity"并计算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;。请注意，您需要计算初始上限和下限，如果它们相等，则跳过所有这些。还要注意，这就是全部 "in theory, completely untested"，我可能在某个地方搞砸了。
我主要想说的是（在我看来）你应该放弃并接受你无法做任何事情来接近 0.5 ulp。对不起..:)
由于我的目标只是扩大获得准确结果的时间间隔，而不是找到适用于 a 所有可能值的解决方案，因此使用双 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、224*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/224、K*2 区间内为 a 提供几乎正确的四舍五入结果24)。不幸的是，与我的问题中的代码相比，这是以十次额外的操作为代价的，这是将最大误差从大约 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=2n 时，在 [K/224, 4.23*K]。总共有 14 或 15 个操作（取决于架构是否支持完整谓词或仅支持条件移动），这比我的原始代码需要多七到八个操作。
最后，残差计算可以直接基于原始变量a以避免m和p计算中固有的错误。这导致以下代码，对于 K = 2n，计算区间 [K/224 中 a 几乎正确的舍入结果, 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);
在各种情况下，例如对于数学函数的参数约简，需要计算 (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的幂，中间结果没有上溢或下溢。对于 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 的参数，改进除法的准确性并不比使用朴素方法好得多。
我真的没有答案（适当的浮点错误分析非常乏味）但有几点观察：
- 快速倒数指令（例如 RCPSS）不如除法准确，因此如果使用它们，您可能会发现精度会降低。
m如果一个 ∈ [0.5×Kb, 21+n×Kb), 其中Kb 是低于 K 的 2 的幂（如果 K 是 2 的幂，则为 K 本身），n 是 K 的尾随零的数量（即，如果 K 是 2 的幂, 那么 n=23).- 这类似于 Dekker (1971) 中
div2算法的简化形式：要扩大范围（尤其是下限），您可能需要从中合并更多校正项（即将m存储为 2float的总和，或使用double）。
如果您可以将 API 放宽为 return 另一个对错误建模的变量，那么解决方案就会变得简单得多：
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 来使用定点。
测试代码更新为人工生成非零最低有效位 在输入
测试代码
如果 a 比 K 大，则 (a-K)/(a+K) = 1 - 2K / (a + K) 将给出一个很好的近似值。如果 a 与 K 相比较小，则 2a / (a + K) - 1 将给出一个很好的近似值。如果 K/2 ≤ a ≤ 2K，那么 a-K 是精确运算，所以做除法会得到一个不错的结果。
一种可能性是用经典的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 个师，但如果我没有搞砸的话应该接近一半 ulp。
但是这些除法可以用 p 的倒数乘法代替，没有任何问题，因为第一个不正确的舍入除法将由余数 r 补偿，而第二个不正确的舍入除法并不重要（修正 q1 的最后一位不会改变任何东西）。
问题是 (a + K) 中的加法。 (a + K) 中的任何精度损失都会被除法放大。问题不在于部门本身。
如果 a 和 K 的指数相同（几乎），则不会丢失精度，并且如果指数之间的绝对差大于尾数大小，则 (a + K) == a（如果a有更大的量级）或(a + K) == K（如果K有更大的量级）。
没有办法阻止这种情况。增加有效数字大小（例如，在 80x86 上使用 80 位 "extended double"）只会稍微扩大 "accurate result range"。要理解原因，请考虑 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 的备选方案。 None 很可能会被接受table.
第一个替代方案涉及为 a 的每个值预计算查找 table（使用 "big real number" 数学），最终（使用一些技巧）约为 2 GiB对于 32 位浮点（对于 64 位浮点完全疯狂）。当然，如果 a 的可能值范围小于 "any positive value that can fit in a 32-bit float"，则查找的大小 table 将减少。
第二种选择是在运行时使用其他东西（"big real number"）进行计算（并转换to/from 32位浮点数）。
第三种选择涉及，"something"（我不知道它叫什么，但它很贵）。将舍入模式设置为"round to positive infinity"并计算temp1 = (a + K); if(a < K) temp2 = (a - K);然后切换到"round to negative infinity"并计算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;。请注意，您需要计算初始上限和下限，如果它们相等，则跳过所有这些。还要注意，这就是全部 "in theory, completely untested"，我可能在某个地方搞砸了。
我主要想说的是（在我看来）你应该放弃并接受你无法做任何事情来接近 0.5 ulp。对不起..:)
由于我的目标只是扩大获得准确结果的时间间隔，而不是找到适用于 a 所有可能值的解决方案，因此使用双 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、224*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/224、K*2 区间内为 a 提供几乎正确的四舍五入结果24)。不幸的是，与我的问题中的代码相比，这是以十次额外的操作为代价的，这是将最大误差从大约 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=2n 时，在 [K/224, 4.23*K]。总共有 14 或 15 个操作（取决于架构是否支持完整谓词或仅支持条件移动），这比我的原始代码需要多七到八个操作。
最后，残差计算可以直接基于原始变量a以避免m和p计算中固有的错误。这导致以下代码，对于 K = 2n，计算区间 [K/224 中 a 几乎正确的舍入结果, 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);