什么技术适用于除法？
对于除法a/b，可以求出余数（余数）：
a = b*q + r
如果你有 fused-multiply-add，这个余数 r 很容易获得
q = a/b ;
r = fma(b,q,-a) ;
同样的 fma 技巧可以应用于乘法：
y = a*b ;
r = fma(a,b,-y) ; // the result is y+r
如果你在乘积 (a0+ra) / (b0+rb) 之后得到两个近似操作数，那么你对 (a0+ra) = q*(b0+rb) + r 感兴趣。
您可以先评估：
q0 = a0/b0 ;
r0 = fma(b0,q0,-a0);
然后将余数近似为：
r = fma(q0,rb,r0-ra);
然后将商更正为：
q = q0 + r/b0;
编辑：如果 fma 不可用怎么办？
我们可以使用 à la Dekker 的精确乘积来模拟 fma，它被分解为 2 个浮点的精确总和，然后使用 Boldo-Melquiond roundToOdd 技巧确保将 3 个浮点的总和精确舍入。
但这将是矫枉过正。我们只使用 fma 来评估残差，因此我们通常使 c 非常接近 -ab。在这种情况下，ab+c 是精确的，我们只有 2 个浮点数来求和，而不是 3 个。
无论如何，我们只是粗略估计了一堆操作的残差，所以这个残差的最后一点不会那么重要。
所以fma可以这样写：
/* extract the high 26 bits of significand */
double upperHalf( double x ) {
double secator = 134217729.0; /* 1<<27+1 */
double p = x * secator; /* simplified... normally we should check if overflow and scale down */
return p + (x - p);
}
/* emulate a fused multiply add: roundToNearestFloat(a*b+c)
Beware: use only when -c is an approximation of a*b
otherwise there is NO guaranty of correct rounding */
double emulated_fma(a,b,c) {
double aup = upperHalf(a);
double alo = a-aup;
double bup = upperHalf(b);
double blo = b-bup;
/* compute exact product of a and b
which is the exact sum of ab and a residual error resab */
double high = aup*bup;
double mid = aup*blo + alo*bup;
double low = alo*blo;
double ab = high + mid;
double resab = (high - ab) + mid + low;
double fma = ab + c; /* expected to be exact, so don't bother with residual error */
return resab + fma;
}
嗯，比一般的模拟 fma 少了一点矫枉过正，但使用为这部分工作提供原生 fma 的语言可能更聪明...