I haven't explained the reason of the error estimation in
Range#step for Float;

double n = (end - beg)/unit;
double err = (fabs(beg) + fabs(end) + fabs(end-beg)) / fabs(unit) * epsilon;

The reason is as follows. (including unicode characters)
This is based on the theory of the error propagation;
http://en.wikipedia.org/wiki/Propagation_of_uncertainty

If f(x,y,z) is given as a function of x, y, z,
¦¤f (the error of f) can be estimated as:

¦¤f^2 = |˘ßf/˘ßx|^2*¦¤x^2 + |˘ßf/˘ßy|^2*¦¤y^2 + |˘ßf/˘ßz|^2*¦¤z^2

This is a kind of `statistical' error.  Instead, `maximum' error
can be expressed as:

¦¤f = |˘ßf/˘ßx|*¦¤x + |˘ßf/˘ßy|*¦¤y + |˘ßf/˘ßz|*¦¤z

I considered the latter is enough for this case.
Now, the target function here is:

n = f(e,b,u) = (e-b)/u

The partial differentiations of f are:

˘ßf/˘ße = 1/u
˘ßf/˘ßb = -1/u
˘ßf/˘ßu = -(e-b)/u^2

The errors of floating point values are estimated as:

¦¤e = |e|*¦Ĺ
¦¤b = |b|*¦Ĺ
¦¤u = |u|*¦Ĺ

Finally, the error is derived as:

¦¤n = |˘ßn/˘ße|*¦¤e + |˘ßn/˘ßb|*¦¤b + |˘ßn/˘ßu|*¦¤u
= |1/u|*|e|*¦Ĺ + |1/u|*|b|*¦Ĺ + |(e-b)/u^2|*|u|*¦Ĺ
= (|e| + |b| + |e-b|)/|u|*¦Ĺ

Masahiro Tanaka