C++ compound assignment and type conversion issue

Question

I am calculating combination(15, 7) in C++.

1. I first used the following code and get the wrong answer due to a type promotion error.

   #include <iostream>
   
   int main()
   {
       int a = 15;
       double ans = 1;
       for(int i = 1; i <= 7; i++)
           ans *= (a + 1 - i) / i;
       std::cout << (int) ans;
   
       return 0;
   }
   

   Output: 2520

2. So I changed ans *= (a + 1 - i) / i; to ans *= (double)(a + 1 - i) / i; and still get the wrong answer.

   #include <iostream>
   
   int main()
   {
       int a = 15;
       double ans = 1;
       for(int i = 1; i <= 7; i++)
           ans *= (double) (a + 1 - i) / i;
       std::cout << (int) ans;
   
       return 0;
   }
   

   Output: 6434

3. Finally, I tried ans = ans * (a + 1 - i) / i, which gives the right answer.

   #include <iostream>
   
   int main()
   {
       int a = 15;
       double ans = 1;
       for(int i = 1; i <= 7; i++)
           ans = ans * (a + 1 - i) / i;
       std::cout << (int) ans;
   
       return 0;
   }
   

   Output: 6435

Could someone tell me why the second one did not work?

Answer 1

If you print out ans without casting it to (int) you'll see the second result is 6434.9999999999990905052982270717620849609375. That's pretty darn close to the right answer of 6535, so it's clearly not a type promotion error any more.

No, this is classic floating point inaccuracy. When you write ans *= (double) (a + 1 - i) / i you are doing the equivalent of:

ans = ans * ((double) (a + 1 - i) / i);

Compare this to the third version:

ans = ans * (a + 1 - i) / i;

The former performs division first followed by multiplication. The latter operates left to right and so the multiplication precedes the division. This change in order of operations causes the results of the two to be slightly different. Floating point calculations are extremely sensitive to order of operations.

Quick fix: Don't truncate the result; round it.

Better fix: Don't use floating point for integral arithmetic. Save the divisions until after all the multiplications are done. Use long, long long, or even a big number library.

Answer 2

First one did not work because you have integer division there.

Difference btw second one and third one is this:

ans = ans * (double(a + 1 - i) / i); // second is equal to this

vs:

ans = (ans * (a + 1 - i)) / i; // third is equal to this

so difference is in order of multiplication and division. If you round double to integer instead of simply dropping fractional part you will get the same result.

std::cout << int( ans + 0.5 ) << std::endl;