K）/（a + K），具有改进的准确度
在各种情况下，例如为了减少数学函数的参数，需要计算(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 ，上面的代码计算所有输入的商几乎正确舍入（最大误差非常接近0.5 ulps），前提是K是2的幂，并且存在中间结果没有溢出或下溢。 对于K不是2的幂，这个代码比基于划分的朴素算法更精确。 就性能而言，此代码可以比平台上的朴素方法更快，因为平台中浮点倒数的计算速度可能比浮点数快。
当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，改进分割的准确性并不比使用朴素方法好得多。
如果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=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的最后一位不会改变任何东西）。
我没有真正的答案（适当的浮点错误分析非常繁琐），但有一些观察结果：
m准确计算，如果在＆; [0.5×Kb，21 + n×Kb），其中Kb是低于K的2的幂（或者K本身如果K是2的幂），并且n是K的有效数中的尾随零的数量（即，if K是2的幂，则n = 23）。 div2算法的简化形式类似：为了扩大范围（特别是下限），您可能必须包含更多的校正项（即存储m作为2 float s，或使用double ）。