Comment by moefh

Nice write-up.

Let me offer a nitpck: in the "Gradual underflow" section it says this about subnormal numbers:

    Bonus: we have now acquired extra armour against a division by zero:

    if ( x != y ) z = 1.0 / ( x - y );

But that's not that useful: just because you're not dividing by zero doesn't mean the result won't overflow to infinity, which is what you get when you do divide by zero.

Think about it this way: the smallest subnormal double is on the order of 10^-324, but the largest double is on the order of 10^308. If `x - y` is smaller than 10^-308, `1.0 / (x - y)` will be larger than 10^308, which can't be represented and must overflow to infinity.

This C program demonstrates this:

    #include <stdio.h>
    #include <float.h>
    #include <math.h>

    // return a (subnormal) number that results in zero when divided by 2:
    double calc_tiny(void)
    {
        double x = DBL_MIN; // the normal number closest to 0.0
        while (1) {
            double next = x / 2.0;
            if (next == 0.0) {
                return x;
            }
            x = next;
        }
    }

    int main(void)
    {
        double y = calc_tiny();
        double x = 2.0 * y;
        if (x != y) {
            double z = 1.0 / (x - y);
            printf("division is safe: %g\n", z);
        } else {
            printf("refusing to divide by zero\n");
        }
    }

(It will print something like "division is safe: inf", or however else your libc prints infinity)