Arithmetic double expression and c/c++

Question

Why this code give me 2 different result on ris1 and ris2?

Why they are not equivalent???

#include <iostream>
#include <stdio.h>
#include <malloc.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>

using namespace std;

int main() {
    double x1 = 0.4628471891711442;
    double h = 0.51152525298500628;
    double lambda = 0.99999999999999989;
    double t = 0.10000000000000001;
    double ris1 = 0;
    double ris2 = 0;

    ris1 = x1 - (h * exp(-lambda * t));

    double sub = h * exp(-lambda * t);
    ris2 = x1 - sub;

    printf("\n\nris1 = %1.4e", ris1);
    printf("\n\nris2 = %1.4e", ris2);

    return 0;
}

Result:

ris1 = 8.9257e-017

ris2 = 1.1102e-016

I am working with mingw in eclipse. Here my mingw installed framework: http://ctrlv.in/638751

Is this all the code? I don't think it is, and the problem has to be somewhere else. Because this code as it is has no problems. Have you tried to test the code in an isolated compilation unit? Note that none of those two values is correct. The answer is `5.5511e-17`. — Iharob Al Asimi, Sep 19 '15 at 21:44
Could you tell us more about the environment in which you run your code? I am unable to reproduce your problem. In particular I would like you to post the version of your standard libraries, the compiler and all the options (especially optimization-related) used for compilation, as well as the data about processor which ran the code. — mszymborski, Sep 19 '15 at 21:47
And c++ although your code is c. There is no way to reproduce your problem. Try to run the program directly in a console window. — Iharob Al Asimi, Sep 19 '15 at 22:20
@mszymborski in the pasted link you can find my mingw configuration from the installation manager. — andrea, Sep 19 '15 at 22:43
Essentially, you are computing the difference between two floating point values that are very close in value (i.e. difference is near zero). That is one way to amplify the effects of rounding error, and increase sensitivity to order of operations. — Peter, Sep 19 '15 at 23:03
@Peter - it has nothing to do with it. Using a proper implementation (see my answer) the only possible difference could be casting a 80-bit value used in computations to a 64-bit value for storage and then promoting it back and using it in the subtraction. Other than that, the errors should propagate the same way in both expressions. — mszymborski, Sep 19 '15 at 23:58
@iharob: The code is C++, as OP included `` and uses `namespace` which do not exist in C. — too honest for this site, Sep 20 '15 at 00:07
@Olaf There is nothing in the code that makes it invalid c except for that and it's not used anywhere. I suspect that the OP does not know the difference. — Iharob Al Asimi, Sep 20 '15 at 00:12
@iharob: I did not say different;-) Question is how to treat such a situation (neither the first one, nor the last). As there is no clear problem statement and the error is likely due to rounding, I think it can be closed. — too honest for this site, Sep 20 '15 at 00:14

score 1 · Answer 1 · edited May 23 '17 at 11:51

1

Given your configuration, it's most likely a bug in implementation of MinGW's g++ 4.8.1-4 , as explained here. Please consider updating (or downgrading) your MinGW installation..

edited May 23 '17 at 11:51

Community

1
1

answered Sep 19 '15 at 23:53

mszymborski

1,615
1
20
28

Or even avoiding MinGW and using an operating system where the bug does not exist for g++. Or using MSVC c++ compiler. – Iharob Al Asimi Sep 20 '15 at 00:17

score 0 · Answer 2 · edited May 23 '17 at 11:58

The difference of x1 - sub is near the last few bits of double precision as the 2 variables are about the same value. @peter

Due to optimization and the degree of correctness of exp(), slightly variant answers are expected.

Another factor is C allows double math to be carried out as double or long double depending on FLT_EVAL_METHOD.

These 2 above may not be used equally with the different computation paths.

By forcing extra precision with long double, the answers are consistent

ris1 = 1.2357e-16
ris2 = 1.2357e-16

Arithmetic double expression and c/c++

2 Answers2