floating-pointprecisionradixnumerical-analysis

How and why does this algorithm find the floating point base value of your computer?


The algorithm in pseudocode:

float a = 1
while (a + 1) - a == 1
    a = 2 * a

int i = 1
while (a + i) == a
    i = i + 1

return (a + i) - a

This will return the base your computer uses (most probably 2). Variable a must be a float or double for it to work.

I don't understand why and how exactly it works.


Solution

  • It's a search algorithm, searching for (nearly) the first place floating points can no longer represent all integers consecutevly. Just to make it easier to read, let's assume our base is 10. What does it mean (a+1) - a != 1? Mark

    a=s*10^e
    

    where s is the significand and e the exponent. Then:

    a+1=s*10^e + 1*10^0 = s*10^e + (1/10^e)*10^e=(s+1/10^e)*10^e
    

    Now 1/10^e = 1 * 10^-e is 0.0000...1 where there are e zeros. This is limited by the machines/language/type precision. When e is large enough, this will just be 0. So the first loop finds one of the first (ish) numbers this happens.

    Now the second loop finds a first integer that adding it to a is something the machine notices, the next representable value of a. Let's initially just guess that this integer is the base, 10. Then we have:

    a + i = s*10^e + 10 = 10*(s*10^(e-1) + 1)
    

    We know the RHS can be represented because e was the first exponent that adding 1 to was not representable (so e-1 is), and multiplying by the base (10) is just adding 1 to the exponent. Lets try a smaller integer:

    a+9=s*10^e + (9/10^e)*10^e = (s + 9/10^e)*10^e
    

    0.000....9 needs the same precision as 1/10^e = 0.000...1 to be different from 0, and hence will not change the value of a. This also makes clear another way to show i=10 is the first integer to be enough - we would have 10/10^e=1/10^(e-1), which again, by virtue of the first loop, is representable.

    Finally, it is enough just to print i.

    Note this works for any base, it's just easier to write in base 10 (where 1/10^e=0.000...1 in a familiar manner). Note also that we don't have to rely on a = 2*a, that is, that a is a power of the base as happens to be the case for you. It would have simplified the above a bit (s=1), but of course, that would be cheating (since we don't know the base a priori).