有效地计算（a-K）/（a + K）并提高精度
在各种情况下，例如对于数学函数的参数减少，需要计算(a - K) / (a + K) ，其中a是正变量参数而K是常数。 在许多情况下， K是2的幂，这是与我的工作相关的用例。 我正在寻找比直接划分更准确地计算这个商的有效方法。 可以假设对融合乘法 – 加法（FMA）的硬件支持，因为此操作由此时所有主要CPU和GPU架构提供，并且可通过函数fma()和fmaf()以C / C ++forms提供。
为了便于探索，我正在尝试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不是2的幂，该代码仍然比基于除法的朴素算法更准确。 在性能方面，这个代码可以比平台上的朴素方法更快 ，在这些平台上，浮点倒数可以比浮点除法更快地计算。
当K = 2 n时，我做了以下观察：当工作区间的上限增加到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的论证，改进除法的准确性并不比使用朴素方法好。
如果a与K相比较大，则（aK）/（a + K）= 1 – 2K /（a + K）将给出良好的近似。 如果a与K相比较小，则2a /（a + K）-1将给出良好的近似值。 如果K /2≤a≤2K，则aK是精确的运算，因此进行除法将得到一个合适的结果。
一种可能性是使用经典的Dekker / Schewchuk跟踪m和p的误差为m1和p1：
m=ak; k0=am; a0=k0+m; k1=k0-k; a1=a-a0; m1=a1+k1; p=a+k; k0=pa; 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的最后几位不会改变任何东西） ）。
我真的没有答案（适当的浮点误差分析非常繁琐），但有一些观察结果：
- 快速倒数指令（如RCPSS ）不如除法精确，因此如果使用这些指令，您可能会看到准确度降低。
- 如果a∈[0.5×K b ，2 1 + n ×K b ），则精确计算
m，其中K b是低于K的2的幂（或者K本身，如果K是2的幂），并且n是数字在K的有效数中的尾随零（即，如果K是2的幂，则n = 23）。 - 这类似于Dekker（1971）的
div2算法的简化forms：为了扩展范围（特别是下限），你可能必须从中加入更多的修正项（即存储m作为2float的总和） s，或使用double）。
由于我的目标只是扩大实现精确结果的时间间隔，而不是找到适用于a的所有可能值的解决方案，因此对所有中间计算使用双float算法似乎成本太高。
更多地思考这个问题，很明显，在我的问题的代码中，除法的其余部分的计算是实现更准确结果的关键部分。 在数学上，余数是（aK） – q *（a + K）。 在我的代码中，我只是使用m来表示（aK）并将（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]提供了几乎正确的舍入结果，允许显着增加达到精确结果的间隔的上限。
加宽下端的间隔需要更准确地表示（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);
详尽的测试表明，这可以在[K / 2 24 ，K * 2 24 ]区间内提供几乎正确的舍入结果。 不幸的是，与我的问题中的代码相比，这需要花费10个额外的操作，这是一个陡峭的代价，从最小的1.625 ulps获得最大误差，天真计算下降到接近0.5 ulp。
正如在我的问题的原始代码中，可以用（aK）表示（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);
如果主要焦点是减小间隔的下限，则这certificate是有利的，这是我在问题中解释的特别关注点。 对单精度情况的详尽测试表明，当K = 2时，对于区间[K / 2 24,4.23 * K]中的a值，产生几乎正确的圆形结果。 总共有14或15次操作（取决于架构是支持完全预测还是只是条件移动），这需要比原始代码多7到8个操作。
最后，可以将残差计算直接基于原始变量a以避免计算m和p固有的误差。 这导致以下代码，对于K = 2 n ，在[K / 2 24 ，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);
如果您可以放松API以返回另一个模拟错误的变量，那么解决方案变得更加简单：
float foo(float a, float k, float *res) { float ret=(ak)/(a+k); *res = fmaf(-ret,a+k,ak)/(a+k); return ret; }
该解决方案仅处理除法的截断误差，但不处理a+k和ak的精度损失。
为了处理这些错误，我想我需要使用双精度或bithack来使用定点。
更新测试代码以在输入中人为地生成非零最低有效位
测试代码
问题是(a + K)的加法。 (a + K)任何精度损失都会被除法放大。 问题不在于分裂本身。
如果a和K的指数相同（几乎）没有精度丢失，并且如果指数之间的绝对差值大于有效数量大小那么(a + K) == a （如果a具有更大的量值）或(a + K) == K （如果K具有更大的幅度）。
没有办法阻止这一点。 增加有效数字大小（例如，在80×86上使用80位“扩展双精度”）仅有助于稍微扩大“精确结果范围”。 要理解为什么，请考虑smallest + largest （其中smallest是最小的正非正规，32位浮点数可以）。 在这种情况下（对于32位浮点数），您需要大约260位的有效位大小才能完全避免精度损失。 做（例如） temp = 1/(a + K); result = a * temp - K / temp; temp = 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-完全疯狂）位浮点）。 当然，如果a的可能值的范围小于“可以适合32位浮点数的任何正值”，则查找表的大小将减小。
第二种方法是在运行时使用其他东西（“大实数”）进行计算（并转换为/从32位浮点转换）。
第三种选择涉及“某事”（我不知道它叫什么，但它很昂贵）。 将舍入模式设置为“舍入到正无穷大”并计算temp1 = (a + K); if(a < K) temp2 = (a - K); temp1 = (a + K); if(a < K) temp2 = (a - K); 然后切换到“舍入到负无穷大”并计算if(a >= K) temp2 = (a - K); lower_bound = temp2 / temp1; 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; a_range = a_upper - a_lower; result = upper_bound * (a_upper - a) / a_range + lower_bound * (a - a_lower) / a_range; 。 请注意，如果它们相等，您将需要计算初始上限和下限并跳过所有这些。 另外要注意的是，这一切都“理论上完全没有经过考验”，我可能会把它搞砸到某个地方。
主要是我所说的（在我看来）你应该放弃并接受你无法做到接近0.5 ulp。 对不起.. 🙂