I cannot reproduce numeric_limits<double>::epsilon()

I have tried the following code:

#include <iostream>
#include <limits>

using namespace std;

int main()
  double x = 1;

  double s = 0;

  for(double y = x; y == x; y = x+(s+=x*1e-20)) {} 
  cout << "my calculation: " << s << endl; 

  cout << "epsilon(): " << numeric_limits<double>::epsilon() << endl;
  return 0;

Output is:

my calculation: 1.1108e-016
epsilon(): 2.22045e-016

Only half! Why?

Only half! Why?

That it is half is a clue, isn't it?
Answer: rounding error

Let f1 be a particular representable floating point number, and f2 be the next representable floating point number. If the IEEE-754 rounding method (round to nearest, ties to even) is used (as almost every current implementation does), then

f1 + delta yields f2 if delta == (f2-f1)/2

See: http://docs.sun.com/source/806-3568/ncg_goldberg.html#689

commented: excellent :) +34
commented: Good one +5