在Windows10和Debian /Linux10中生成不同输出的数值过程

Question

我在玩Maehly的程序，以擦亮一个多项式的根部，偶然发现了一些东西:精确的样本代码给了我两个不同的输出，取决于它是在哪台机器上编译的。

代码

#include <stdio.h>

#define MAX_ITERATION 1000

double poly(double x){
    double coeff[9]={-61.688, 72.5235, 72.822, -108.519, -5.12949, 39.9139,-7.07373, -3.91823, 1.0};
    double result=coeff[0];
    double buffer;
    
    for(int i=1; i<9;i++){
        buffer=coeff[i];
        for(int j=1;j<=i;j++){
            buffer*=x;
        }
        result+=buffer;
    }
    return result;

}
double poly_der(double x){
    double coeff[8]={ 72.5235, 72.822, -108.519, -5.12949, 39.9139,-7.07373, -3.91823, 1.0};
    double result=coeff[0];
    double buffer;
    
    for(int i=1; i<8;i++){
        buffer=coeff[i]*(i+1);
        for(int j=1;j<=i;j++){
            buffer*=x;
        }
        result+=buffer;
    }
    return result;
}

int main(){
    double roots[8]={0.9, -1.1, 1.4, 1.4, -2.0, -2.0, 2.2, 2.2};
    double factor;
    double pol_eval;
    //Implement Maehly-procedure

    for(int i=0; i<MAX_ITERATION;i++){ 
        for(int k=0;k<8;k++){
            factor=0;
            for(int j=0;j<k;j++){
                factor+=1/(roots[k]-roots[j]);
            }
            pol_eval=poly(roots[k]);
            roots[k]-=pol_eval/(poly_der(roots[k])-(pol_eval*factor));

        }
    }
   

    for(int i=0;i<8;i++){
        printf("\n%d: x:%0.16f poly:%e \n",i,roots[i],poly(roots[i]));
    }
}

Windows输出(Windows10):

0: x:1.0072928773885637 poly:-8.437695e-015 

1: x:-1.0004044550991309 poly:-2.375877e-014 

2: x:1.3770602924650244 poly:-3.552714e-015  

3: x:-2.5000428878301499 poly:0.000000e+000  

4: x:-1.7318124315476966 poly:-1.136868e-013

5: x:3.0001628929552053 poly:9.094947e-013

6: x:2.2341265341600458 poly:-2.273737e-013

7: x:3.0001628929552049 poly:0.000000e+000

Linux输出(Debian/Linux10)：

0: x:1.0072928773885637 poly:-8.437695e-15

1: x:-1.0004044550991309 poly:-2.375877e-14

2: x:1.3770602924650244 poly:-3.552714e-15

3: x:-2.5000428878301499 poly:0.000000e+00

4: x:-1.7318124315476959 poly:2.842171e-14

5: x:3.0001628929552093 poly:-1.818989e-12

6: x:2.2341265341600458 poly:-2.273737e-13

7: x:1.5318471775081237 poly:0.000000e+00

X是polinomial的抛光根，开始值保存在数组roots[8]中。

你能帮我解释一下这种行为吗?最重要的是，你能帮我理解如何避免将来发生类似的事情吗？

Answer 1

我们有两个问题，为什么不同？哪条路更好--还是第三条路呢？

OP reports不同的FLT_EVAL_METHOD值为2和0。这意味着Windows版本正在使用long double数学进行中间计算，而Linux版本则只是使用double。通常情况下，FLT_EVAL_METHOD==2更正确。

#include <float.h>
printf("%d\n", FLT_EVAL_METHOD);

FP在减去值方面有一个弱点，几乎相同。许多常见比特的抵消导致微小的计算错误变得更加重要。由于各种原因，不同的编译可能有略微不同的舍入结果，尽管我的_expect也是如此。poly()通过反复计算找到x的幂，然后添加抵消项。

在我的系统中，FLT_EVAL_METHOD是0。

修复poly()和poly_dev()

static const double coeff[9] = {-61.688, //
    72.5235, 72.822, -108.519, -5.12949, //
    39.9139, -7.07373, -3.91823, 1.0};

double poly(double x) {
  double result = 0.0;
  for (int i = 9; i-- > 0; ) {
    result = result * x + coeff[i];
  }
  return result;
}

double poly_der(double x) {
  double result = 0.0;
  for (int i = 9; i-- > 1; ) {
    result = result * x + coeff[i]*i;
  }
  return result;
}

结果只有double数学可以与使用long double数学的较弱的计算代码相媲美。

printf("%d: x:% 0.16f poly:% e \n", i, roots[i], poly(roots[i]));
0: x: 1.0072928773885645 poly: 0.000000e+00 
1: x:-1.0004044550991309 poly:-7.105427e-15 
2: x: 1.3770602924650324 poly:-4.973799e-14 
3: x:-2.5000428878301495 poly:-6.394885e-13 
4: x:-1.7318124315476964 poly:-2.842171e-14 
5: x: 1.5318471775081377 poly: 0.000000e+00 
6: x: 2.2341265341600520 poly: 0.000000e+00 
7: x: 3.0001628929552009 poly:-2.771117e-13