Why 10/3 it's exact in C?

Question

Take a look on this code. 10/3 return 3.3333332538604736328125000 and when I multiply by 3 in a calcutor i get 9.99, but if do the same by the code i get exactly 10.00. How it's posible ?

#include <stdlib.h>
#include <stdio.h>

int main() {

    float v = 10.f/3.f;
    float test = v*3.f;
    printf("10/3 => %25.25f \n (10/3)*3 => %25.25f\n",v,test);
    return 0;
}

This is the assembly code without printf, compiled using default gcc 7.2.1 parameters:

0000000000400497 <main>:
  400497:       55                      push   rbp
  400498:       48 89 e5                mov    rbp,rsp
  40049b:       f3 0f 10 05 b1 00 00    movss  xmm0,DWORD PTR [rip+0xb1]        # 400554 <_IO_stdin_used+0x4>
  4004a2:       00 
  4004a3:       f3 0f 11 45 fc          movss  DWORD PTR [rbp-0x4],xmm0
  4004a8:       f3 0f 10 4d fc          movss  xmm1,DWORD PTR [rbp-0x4]
  4004ad:       f3 0f 10 05 a3 00 00    movss  xmm0,DWORD PTR [rip+0xa3]        # 400558 <_IO_stdin_used+0x8>
  4004b4:       00 
  4004b5:       f3 0f 59 c1             mulss  xmm0,xmm1
  4004b9:       f3 0f 11 45 f8          movss  DWORD PTR [rbp-0x8],xmm0
  4004be:       b8 00 00 00 00          mov    eax,0x0
  4004c3:       5d                      pop    rbp
  4004c4:       c3                      ret    
  4004c5:       66 2e 0f 1f 84 00 00    nop    WORD PTR cs:[rax+rax*1+0x0]
  4004cc:       00 00 00 
  4004cf:       90                      nop

I think mulss is rounded by a CPU feature.

For note, 10/3 in the GNU BC program returns 3.3333333333333333333333 ( *3 => 9.9999) and in SciLab returns 3.3333333333333334813631 ( *3 => 10).

Comments are not for extended discussion; this conversation has been [moved to chat](http://chat.stackoverflow.com/rooms/160240/discussion-on-question-by-amanda-osvaldo-why-10-3-its-exact-in-c). — Andy, Nov 30 '17 at 22:44
The `mulss` result is rounded using the default IEEE754 rounding mode, which is round to nearest, with even as a tie-break. (Your program doesn't use `#pragma FENV_ACCESS ON` and set the FP rounding mode to something else, so the compiler produces a program that uses the default rounding mode.) — Peter Cordes, Dec 01 '17 at 01:30

dbush · Accepted Answer · 2017-12-01T01:22:15.783

You end up getting exactly 10 as a result because the representation happens to work out that way. I get the same on my implementation for both float and double.

Let's look at an example using double:

If we print out 10./3. in hexadecimal floating point notation using %a, we get this:

0x1.aaaaaaaaaaaabp+1

This matches up with the IEEE754 double representation 0x401aaaaaaaaaaaab.

The above number normalized is:

0x3.5555555555558

In binary:

11.0101010101010101010101010101010101010101010101011

To keep things simple, let's add three times instead of multiplying by 3:

     11.0101010101010101010101010101010101010101010101011
+    11.0101010101010101010101010101010101010101010101011
---------------------------------------------------------
    110.1010101010101010101010101010101010101010101010111
+    11.0101010101010101010101010101010101010101010101011
---------------------------------------------------------
   1010.0000000000000000000000000000000000000000000000000

Which is exactly 10.

EDIT:

Looks like I managed to botch the math on the last few digits. The actual sum:

     11.0101010101010101010101010101010101010101010101011
+    11.0101010101010101010101010101010101010101010101011
---------------------------------------------------------
    110.1010101010101010101010101010101010101010101010110
+    11.0101010101010101010101010101010101010101010101011
---------------------------------------------------------
   1010.0000000000000000000000000000000000000000000000001

So it's not exactly 10, but off by the least significant bit.

I noticed a similar difference when using float.

10.f/3.f printed with %a:

0x1.aaaaaap+1

Normalized:

0x3.555554

In binary:

11.0101010101010101010101

Then we add:

     11.0101010101010101010101
+    11.0101010101010101010101
------------------------------
    110.1010101010101010101010
+    11.0101010101010101010101
------------------------------
   1001.1111111111111111111111

Again, off by the least significant bit.

As for how the actual result is rounded, that I can't say.

Is that first addition right? The rightmost bits are ones, shouldn't `1 + 1` equal `10`, instead of `11` as shown here? And I think the addition should leave a stray `1` somewhere in the rightmost place. It'll probably get rounded off then, but... — ilkkachu, Nov 30 '17 at 22:55
@ilkkachu: FP multiplication doesn't round until the end, so any attempt to recreate it with addition needs to keep extra bits. Although I think this is supposed to be plain binary math, not exponent / mantissa (with the radix point in the right place to separate integer from fractional part). I agree this looks fishy, but the right fix is to keep extra precision. — Peter Cordes, Dec 01 '17 at 01:04
@PeterCordes, yeah, I meant rounding off at the end. The result has more significant bits in the integer part, so some bits should round off from the fractional part — ilkkachu, Dec 01 '17 at 08:54

score 2 · Answer 2 · answered Nov 30 '17 at 21:50

The reason for the discrepancy between what you see in C and what you see in SciLab is that you are using single-precision floating-point values (float) in C, whereas SciLab seems to use double-precision values (double) by default.

You can see the difference here (just remove the f suffix from your numbers and put double instead of float).

Why 10/3 it's exact in C?

2 Answers2