This floating-point precission error surprised me

> // so you would expect Zero

> // or 1e-1000 ....

The smallest number that can be represented with a normal double precision number is approximately 1e-308 (with subnormal numbers you can get 1e-324), so 1e-1000 is completely unrealistic.

But each calculation is done with a precision of approximately 15 digits, not 308 digits. In particular, if you see 0.87 the actual number is something like

  0.8700000000000002365426121879219
  
  0.87000000000000094383991287

  0.87000000000000026190262019827

Moreover, each calculation is rounded with approximately a precision of 15 digits. So the most optimistic possibility is to get a result that has approximately 15 accurate digits, not 300 or 1000 accurate digits

But in this example you are repeating the addition N = 4000000 times. With some hand waving and unsupported assumptions, you can estimate that if each sum has a relative error of 1e-15 and you sum 4e6 of them, the final relative error will be 4e-9. [You shouldn't add relative errors ...]

The final value of RxN in the example is 3480000.00030716, so the relative error is 0.00030716/3480000 = 8.8e-11 that is 100 smaller that my easy hand waving argument. So it's time to use wild guess ...

a) You are lucky and the error is smaller. ???

b) The rounding is sometimes upwards and sometimes downwards (look the details in the standard) and they somehow cancel. [If they were random the error would grow as sqrt(N) instead of N.] ???

c) In the first 1000000 iterations the value of RxN is smaller, so the same relative error gives a smaller absolute error. So we can ignore the error in the first millions of iterations. ???

I'm not sure if a), b) or c) is the better explanation, but an error of 1/1e11 is not surprising for me.