Numerical Procedure generating different outputs in Windows10 vs Debian GNU/Linux 10-CodePudding

I was playing arround with the Maehly Procedure to polish the roots of a polinomial and stumbled upon something interessing: The exact samame code gave me two really different outputs depending on the machine it was compiled on.

The Code

#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]));
    }
}

The Windows output (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

The Linux Output(Debian GNU/Linux 10):

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

The x are the polished roots of the polinomial, start values are saved in the array roots[8].

Can you help me explain this behavior and,most important, help me understand how to avoid something similar in the future?

CodePudding user response：

We have to 2 issue, why different? Which is better - or may a 3rd way?

OP reports different FLT_EVAL_METHOD values of 2 and 0. This implies the Windows version is using long double math for intermediate calculations and the Linux one is just using double. Usually the FLT_EVAL_METHOD==2 is more correct.

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

FP have a weakness in subtracting values nearly the same. The cancelation of many common bits leads to the slight computational errors becoming more significant. Different compilations may have slightly different rounded result for various reasons, although I'd _expect the same. poly() does that with repeated calculations to find a power of x and then adding terms that cancel out.

Alternative code below does far fewer calculations are certainly provide a better result for both systems.

double poly(double x) {
  const double coeff[9] = {-61.688, //
      72.5235, 72.822, -108.519, -5.12949, //
      39.9139, -7.07373, -3.91823, 1.0};
  double result = 0.0;
  for (int i = 9; i-- > 0; ) {
    result = result * x   coeff[i];
  }
  return result;
}

On my system FLT_EVAL_METHOD is 0.
Result of just changing poly(). Similar issue applies to poly_der().

0: x:1.0072928773885645 poly:0.000000e 00 
1: x:-1.0004044550991309 poly:-7.105427e-15 
2: x:1.3770602924650401 poly:-2.842171e-14 
3: x:3.0001628929552013 poly:3.552714e-13 
4: x:-1.7318124315476964 poly:-2.842171e-14 
5: x:-2.5000428878301499 poly:6.821210e-13 
6: x:2.2341265341600485 poly:-7.105427e-15 
7: x:1.5318471775081424 poly:3.552714e-14

Fixing both poly() and 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;
}

Results with only double math rival the weaker computational code using long double math.

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