How to calculate 32-bit floating-point epsilon?

Question

In book Game Engine Architecture : "..., let’s say we use a floating-point variable to track absolute game time in seconds. How long can we run our game before the magnitude of our clock variable gets so large that adding 1/30th of a second to it no longer changes its value? The answer is roughly 12.9 days." Why 12.9 days, how to calculate it ?

Answer 1

When the result of a floating point computation can't be represented exactly, it is rounded to the nearest value. So you want to find the smallest value x such that the increment f = 1/30 is less than half the width h between x and the next largest float, which means that x+f will round back to x.

Since the gap is the same for all elements in the same binade, we know that x must be the smallest element in its binade, which is a power of 2.

So if x = 2^k, then h = 2^k-23 since a float has a 24-bit significand. So we need to find the smallest integer k such that

2^k-23/2 > 1/30

which implies k > 19.09, hence k = 20, and x = 2²⁰ = 1048576 (seconds).

Note that x / (60 × 60 × 24) = 12.14 (days), which is a little bit less that what your answer proposes, but checks out empirically: in Julia

julia> x = 2f0^20
1.048576f6

julia> f = 1f0/30f0
0.033333335f0

julia> x+f == x
true

julia> p = prevfloat(x)
1.04857594f6

julia> p+f == p
false

UPDATE: Okay, so where did the 12.9 come from? The 12.14 is in game time, not actual time: these will have diverged due to the rounding error involved in floating point (especially near the end, when the rounding error is actually quite large relative to f). As far as I know, there's no way to calculate this directly, but it's actually fairly quick to iterate through 32-bit floats.

Again, in Julia:

julia> function timestuff(f)
           t = 0
           x = 0f0
           while true
               t += 1
               xp = x
               x += f
               if x == xp
                   return (t,x)
               end
           end
       end
timestuff (generic function with 1 method)

julia> t,x = timestuff(1f0/30f0)
(24986956,1.048576f6)

x matches our result we calculated earlier, and t is the clock time in 30ths of a second. Converting to days:

julia> t/(30*60*60*24)
9.640029320987654

which is even further away. So I don't know where the 12.9 came from...

UPDATE 2: My guess is that the 12.9 comes from the calculation

y = 4 × f / ε = 1118481.125 (seconds)

where ε is the standard machine epsilon (the gap between 1 and the next largest floating point number). Scaling this to days gives 12.945. This provides an upper bound on x, but it is not the correct answer as explained above.