除以整数时最小化舍入误差
Minimize rounding error when dividing by integers
我试图通过取一个整数乘积除以另一个整数乘积的比值来形成双精度浮点数（64 位）。我希望以减少舍入误差的方式这样做。
我熟悉加减法的卡汉求和。哪些技术适用于除法？
分子是许多长值（数万）的乘积，分母也是如此。我也希望防止溢出和下溢。 （一个应用程序通过在足够数量的项后停止来估计无限乘积。）
我尝试过的一件事是对容易因式分解的数字进行因式分解（使用最多一百万的已知质数进行试除）并取消公因数，这有帮助，但还不够。我的错误大约是 1.0E-13。
我正在使用 C#，但欢迎使用 IEEE 标准浮点数的任何代码。
研究：
我看到一篇讨论 + - x / 的 EFT（无误差变换）、霍纳法则（多项式）和平方根的好论文。标题是 Philippe Langlois 的“4ccurate 4lgorithms in Floating Point 4rithmetic”。参见 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
您要查找的 Kahan 求和的乘法等效项是“双双乘法”。在这里，如果您的整数可以表示为 double 值，那么 crlibm 中的函数 Mul122 就足够了。
#define Mul122(resh,resl,a,bh,bl) \
{ \
double _t1, _t2, _t3, _t4; \
\
Mul12(&_t1,&_t2,(a),(bh)); \
_t3 = (a) * (bl); \
_t4 = _t2 + _t3; \
Add12((*(resh)),(*(resl)),_t1,_t4); \
}
bh 和 bl 是 运行 乘积，以两个 double 值之和的形式存储，精度更高。 a 是下一个整数（我们假设它被精确地转换为 double）。 resh 和 resl 收到下一个 运行 产品，其中考虑了因素 a。
为了避免下溢和上溢，您可以将指数外部化为您希望宽度的整数。这是通过定期将 frexp 函数应用于 运行 乘积的高部分，然后通过将两个分量除以相同的 2 次幂来归一化 运行 乘积（跟踪总运行 乘积除以 2 的幂可以在旁边用所需宽度的整数变量完成。
应用的频率 frexp 取决于您乘以整数的界限。如果整数低于 253，这将有助于将它们精确表示为 double 值，您可以在必须标准化 运行 之前进行大约 19 次乘法运算乘积，因为双精度指数上升到 1023。
计算完分子和分母对应的乘积后，丢弃低分量，除以高分量。这只会引入大约 1ULP 的误差。您的目标不是小于双精度 ULP 的误差，是吗？
不要忘记你在旁边留下的分子和分母的2的幂！将它们相减并使用 ldexp 函数将差值应用于商。
什么技术适用于除法？
对于除法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 不可用怎么办？
我们可以使用 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 的语言可能会更聪明...
除法不会遭受与加法和减法相同的灾难性抵消效应，并且使用 IEEE 浮点数是正确舍入的，因此应该有大约 1/2 ulps (~2e-16) 的相对误差。任何大于该值的错误很可能是中间产品的结果，因此需要小心这些。
Dekker (1971) 有一些算法可以扩展初等数学运算的精度：正如另一个答案所指出的，如果您可以访问 fma 运算，则可以简化这些算法。
如果您可以访问 FMA（融合乘加），则其他答案很好，但 C# 不使用它。我继续寻找快速解决方案，但我找到了一个准确的解决方案。
第一步：分别收集分子和分母。
第 2 步：去掉符号并计算负数的乘数以了解答案的符号。
第 3 步：遍历所有数字，计算每个数字的自然对数。
第四步：分别累加分子和分母对数的补偿和。 （使用 Kahan 求和。）
第 5 步：取两个总和之间的差值并计算指数。
第 6 步：恢复标志。
我针对分子中的 100,000 个随机整数和分母中的相同数字进行了测试，但两组都以不同的随机顺序洗牌。如果我使用常规乘法和除法的简单方法，我的累积误差约为 2x10^-15。使用我的补偿日志方法，错误为零。 （也许我很幸运？）我将对更难的案例进行更多测试。尽管如此，通过补偿日志的总和，我在最后舍入之前得到了几乎两倍的精度。
我很惊讶它的效果如此之好。显然执行200,000次对数并不理想。
理论笔记：
累积舍入误差就像随机游走。在 N 次计算之后，您可以预期出现 sqrt(N)*ULP/2 的错误。如果 ULP/2 是 5.0E-18 并且 N 是 200,000，那么你会得到 2.2E-15，这接近于我用天真的方法得到的结果。
我试图通过取一个整数乘积除以另一个整数乘积的比值来形成双精度浮点数（64 位）。我希望以减少舍入误差的方式这样做。
我熟悉加减法的卡汉求和。哪些技术适用于除法？
分子是许多长值（数万）的乘积，分母也是如此。我也希望防止溢出和下溢。 （一个应用程序通过在足够数量的项后停止来估计无限乘积。）
我尝试过的一件事是对容易因式分解的数字进行因式分解（使用最多一百万的已知质数进行试除）并取消公因数，这有帮助，但还不够。我的错误大约是 1.0E-13。
我正在使用 C#，但欢迎使用 IEEE 标准浮点数的任何代码。
研究：
我看到一篇讨论 + - x / 的 EFT（无误差变换）、霍纳法则（多项式）和平方根的好论文。标题是 Philippe Langlois 的“4ccurate 4lgorithms in Floating Point 4rithmetic”。参见 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
您要查找的 Kahan 求和的乘法等效项是“双双乘法”。在这里，如果您的整数可以表示为 double 值，那么 crlibm 中的函数 Mul122 就足够了。
#define Mul122(resh,resl,a,bh,bl) \
{ \
double _t1, _t2, _t3, _t4; \
\
Mul12(&_t1,&_t2,(a),(bh)); \
_t3 = (a) * (bl); \
_t4 = _t2 + _t3; \
Add12((*(resh)),(*(resl)),_t1,_t4); \
}
bh 和 bl 是 运行 乘积，以两个 double 值之和的形式存储，精度更高。 a 是下一个整数（我们假设它被精确地转换为 double）。 resh 和 resl 收到下一个 运行 产品，其中考虑了因素 a。
为了避免下溢和上溢，您可以将指数外部化为您希望宽度的整数。这是通过定期将 frexp 函数应用于 运行 乘积的高部分，然后通过将两个分量除以相同的 2 次幂来归一化 运行 乘积（跟踪总运行 乘积除以 2 的幂可以在旁边用所需宽度的整数变量完成。
应用的频率 frexp 取决于您乘以整数的界限。如果整数低于 253，这将有助于将它们精确表示为 double 值，您可以在必须标准化 运行 之前进行大约 19 次乘法运算乘积，因为双精度指数上升到 1023。
计算完分子和分母对应的乘积后，丢弃低分量，除以高分量。这只会引入大约 1ULP 的误差。您的目标不是小于双精度 ULP 的误差，是吗？
不要忘记你在旁边留下的分子和分母的2的幂！将它们相减并使用 ldexp 函数将差值应用于商。
什么技术适用于除法？
对于除法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 不可用怎么办？
我们可以使用 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 的语言可能会更聪明...
除法不会遭受与加法和减法相同的灾难性抵消效应，并且使用 IEEE 浮点数是正确舍入的，因此应该有大约 1/2 ulps (~2e-16) 的相对误差。任何大于该值的错误很可能是中间产品的结果，因此需要小心这些。
Dekker (1971) 有一些算法可以扩展初等数学运算的精度：正如另一个答案所指出的，如果您可以访问 fma 运算，则可以简化这些算法。
如果您可以访问 FMA（融合乘加），则其他答案很好，但 C# 不使用它。我继续寻找快速解决方案，但我找到了一个准确的解决方案。
第一步：分别收集分子和分母。
第 2 步：去掉符号并计算负数的乘数以了解答案的符号。
第 3 步：遍历所有数字，计算每个数字的自然对数。
第四步：分别累加分子和分母对数的补偿和。 （使用 Kahan 求和。）
第 5 步：取两个总和之间的差值并计算指数。
第 6 步：恢复标志。
我针对分子中的 100,000 个随机整数和分母中的相同数字进行了测试，但两组都以不同的随机顺序洗牌。如果我使用常规乘法和除法的简单方法，我的累积误差约为 2x10^-15。使用我的补偿日志方法，错误为零。 （也许我很幸运？）我将对更难的案例进行更多测试。尽管如此，通过补偿日志的总和，我在最后舍入之前得到了几乎两倍的精度。
我很惊讶它的效果如此之好。显然执行200,000次对数并不理想。
理论笔记：
累积舍入误差就像随机游走。在 N 次计算之后，您可以预期出现 sqrt(N)*ULP/2 的错误。如果 ULP/2 是 5.0E-18 并且 N 是 200,000，那么你会得到 2.2E-15，这接近于我用天真的方法得到的结果。