How to add two large double precision numbers in c++

Question

I have the following piece of code

#include <iostream>
#include <iomanip>

int main()

{
    double  x = 7033753.49999141693115234375;
    double  y = 7033753.499991415999829769134521484375;
    double z = (x+ y)/2.0;

    std::cout  << "y is " << std::setprecision(40) << y << "\n";
    std::cout  << "x is " <<  std::setprecision(40) << x << "\n";
    std::cout  << "z is " << std::setprecision(40) << z << "\n";

    return 0;
}

When the above code is run I get,

y is 7033753.499991415999829769134521484375
x is 7033753.49999141693115234375
z is 7033753.49999141693115234375

When I do the same in Wolfram Alpha the value of z is completely different

 z = 7033753.4999914164654910564422607421875 #Wolfram answer

I am familiar with floating point precision and that large numbers away from zero can not be exactly represented. Is that what is happening here? Is there anyway in c++ where I can get the same answer as Wolfram without any performance penalty?

Answer 1

large numbers away from zero can not be exactly represented. Is that what is happening here?

Yes.

Note that there are also infinitely many rational numbers that cannot be represented near zero as well. But the distance between representable values does grow exponentially in larger value ranges.

Is there anyway in c++ where I can get the same answer as Wolfram ...

You can potentially get the same answer by using long double. My system produces exactly the same result as Wolfram. Note that precision of long double varies between systems even among systems that conform to IEEE 754 standard.

More generally though, if you need results that are accurate to many significant digits, then don't use finite precision math.

... without any performance penalty?

No. Precision comes with a cost.