CODEWITHSUNDEEP
c++
user2509143
user2509143

Reputation: 31

g++ and double precision

Can anyone explain to me why executing the following code returns that (1./3.)*3. is equal to 1.? I was assuming that the rounding would prevent such an equality but apparently it doesn't work that way. Any explanation?

#include <iostream>

int main()
{
  double x = 1./3.;
  double y = 3.*x;
  std::cout.precision(30);
  std::cout << "x = " << x << std::endl;
  std::cout << "y = " << y << std::endl;
  if (y == 1.)
    std::cout << "Equality!" << std::endl;
  else
    std::cout << "Not equality." << std::endl;

  return 0;
}

It returns:

$ g++ test_dp.cpp -o test_dp.exe

$ ./test_dp.exe

x = 0.333333333333333314829616256247

y = 1

Equality!

Upvotes: 0

Views: 217

Answers (1)

TonyK
TonyK

Reputation: 17114

Assuming your double type is IEEE 754 compliant, the hexadecimal representation of 1.0/3.0 is 3FD5555555555555. This is not exactly 1/3. It breaks down as:

  • Sign = 0 (i.e. number is positive)
  • Biased exponent = 3FD (i.e. exponent = -2)
  • 53-bit Mantissa = 15555555555555 (the leading 1 is implied)

Multiplying 15555555555555 by 3 gives 3FFFFFFFFFFFFF. To squeeze this into 53 bits, we have to discard the least significant bit, which (in the default rounding mode for most environments) means that the mantissa is "rounded to even": trailing binary 01 goes to 0, and trailing binary 11 goes to 100 (which means that carries may be propagated). In this case, the carry propagates all the way up to the most significant bit, leaving 40000000000000 in the mantissa. We want bit 52 to be the most significant bit, so we shift this down to 10000000000000 and increase the exponent by 2. The end result is 3FF0000000000000, or exactly 1.0.

So your program is behaving correctly -- an inexact result multiplied by a small integer gives you an exact result. But one could take the view that this is more by luck than judgment :-)

Upvotes: 3

Related Questions