高效计算 （a - K） / （a + K），精度更高

我真的没有答案（适当的浮点误差分析非常乏味），但有一些观察结果：

• 快速倒数指令（如 RCPSS）不如除法准确，因此如果使用这些指令，可能会降低准确性。
• m如果 a ∈ [0.5×K b， 2¹⁺ⁿ×K b），则精确计算，其中 K_b 是 K 低于 K 的 2 的幂（如果 K 是 2 的幂，则为 K 本身），n 是 K 有效数中尾随零的数目（即，如果 K 是 2 的幂，则 n=23）。
• 这类似于 Dekker （1971） 算法的简化形式：为了扩大范围（尤其是下限），您可能需要从中加入更多修正项（即存储为 2 秒的总和，或使用 ）。div2mfloatdouble

如果可以放宽 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+ka-k

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

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

测试代码

https://ideone.com/bHxAg8

如果 a 与 K 相比很大，则 （a-K）/（a+K） = 1 - 2K / （a + K） 将给出一个很好的近似值。如果 a 与 K 相比较小，则 2a / （a + K） - 1 将给出一个很好的近似值。如果 K/2 ≤ ≤ 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 个师，但如果我没有搞砸，应该接近一半。

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

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

如果 和 的指数相同（几乎），则不会丢失精度，并且如果指数之间的绝对差大于有效大小，则（如果具有更大的量级）或（如果具有更大的量级）。aK(a + K) == aa(a + K) == KK

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

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

第一种选择是为每个值预先计算一个查找表（使用“大实数”数学），这（通过一些技巧）最终对于 32 位浮点约为 2 GiB（对于 64 位浮点则完全疯狂）。当然，如果 的可能值范围小于“任何可以容纳在 32 位浮点数中的正值”，则查找表的大小将减小。aa

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

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

我主要想说的是（在我看来）你应该放弃并接受你无能为力来接近 0.5 ulp。不好意思。。:)

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

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

以相对较小的额外计算成本，（a+K）可以表示为双精度，即头尾对，这导致了我的原始代码的以下修改版本：floatp: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）。我们可以将其计算为双头尾对，这会导致以下代码变体：floatm: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²⁴] 提供几乎正确的四舍五入结果。不幸的是，与我问题中的代码相比，这需要花费十个额外的操作，这是一个高昂的代价，要从大约 1.625 ulps 获得最大误差，而朴素的计算下降到接近 0.5 ulp。a

就像我在问题中的原始代码一样，可以用 （a-K） 表示 （a+K），从而消除了 的尾部计算。此方法生成以下代码：pplo

/* 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，^{24， 4.23}*K] 的值，会产生几乎正确的四舍五入结果。总共有 14 或 15 个操作（取决于架构是支持完全预测还是仅支持条件移动），这需要比我的原始代码多 7 到 8 个操作。a

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

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