floating-point - 除以整数时最小化舍入误差

我试图通过将一个整数的乘积除以另一个整数的乘积来形成 double 浮点数(64位)。我希望这样做可以减少舍入误差。

我熟悉用于加法和减法的Kahan求和。哪些技术适用于除法？

分子是许多长值(成千上万)的乘积，分母也是如此。我也希望防止上溢和下溢。 (一个应用程序通过在足够数量的项之后停止来估计无限乘积。)

我尝试过的一件事是将容易分解的数字分解为因子(使用按已知的素数进行的除法运算最多可达一百万)并取消公因子，这虽然有帮助，但还不够。我的错误大约是1.0E-13。

我正在C#中工作，但是欢迎使用任何与IEEE标准浮点数兼容的代码。

研究:

我碰到了一篇很好的论文，讨论了+-x/的EFT(无错误转换)，霍纳法则(多项式)和平方根。标题是Philippe Langlois的“浮点4算术中的4精确4算法”。参见http://www.mathematik.hu-berlin.de/~gaggle/S09/AUTODIFF/projects/papers/langlois_4ccurate_4lgorithms_in_floating_point_4rithmetic.pdf

上面的内容将我指向Karp和Markstein(进行分组):https://cr.yp.to/bib/1997/karp.pdf

最佳答案

哪些技术适用于除法？

对于a/b除法，您可以评估残差(余数):

a = b*q + r

如果您有fusion-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不可用怎么办？

我们可以使用精确乘积à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的矫kill过正少了一点，但是使用一种为本部分工作提供本地fma的语言可能会更聪明...