50 Real numbers — the mathematics of approximation

50.1 Choosing a type, and comparing

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

/* absolute error: used for comparisons near 0 */
bool near_abs(double a, double b, double eps)
{
    return fabs(a - b) < eps;
}

/* relative error: as values grow the spacing grows, so the tolerance scales with size */
bool near_rel(double a, double b, double rel)
{
    double scale = fabs(a) > fabs(b) ? fabs(a) : fabs(b);
    return fabs(a - b) <= rel * scale;
}

int main(void)
{
    double sum = 0.1 + 0.2;

    printf("0.1 + 0.2 == 0.3 ?      %d\n", sum == 0.3);
    printf("compare with absolute eps? %d\n", near_abs(sum, 0.3, 1e-9));
    printf("%.20f\n", sum);
    printf("%.20f\n", 0.3);

    double big = 1e16;
    printf("1e16 + 1 == 1e16 ?      %d  (gap wider than 1)\n", big + 1 == big);
    printf("compare with relative eps? %d\n", near_rel(big + 1, big, 1e-12));
    return 0;
}

0.1 + 0.2 == 0.3 ?      0
compare with absolute eps? 1
0.30000000000000004441
0.29999999999999998890
1e16 + 1 == 1e16 ?      1  (gap wider than 1)
compare with relative eps? 1

The first three lines are the execution check of chapter 8′s first incident — == is false, and printed to twenty digits the two numbers diverge at the end. Use near_abs (absolute error) and it becomes true.

The latter part checks the third incident (absorption) — add 1 to ${10}^{16}$ and the value is unchanged, because at that size the gap between representable neighbours is already wider than 1 (chapter 8′s tick calculation). Here the right tool is not absolute error but near_rel (relative error).

Reduced to practical rules there are three. The default is double — as chapter 8 showed, the room in precision is of a different order, and float is chosen only where memory and bandwidth are tight. Do not use == — code asking whether two reals are equal is almost always suspect (there are exceptions, such as comparing against an integer or against 0, but they must be judged consciously). The tolerance comes from the problem — it is decided by looking at the nature of the calculation and the size of the values; there is no magic constant.

50.2 Special values — infinity and NaN

IEEE 754 (chapter 8) defines special values besides ordinary numbers. Infinity (positive and negative) comes out of overflow or division by zero, and NaN (Not a Number) means “not a number” — the result of an undefinable operation such as $\frac{0}{0}$ or the square root of a negative.

One thing has to be stated exactly here. Integer division by zero is outside the contract (chapter 28), and floating-point division by zero is commonly said to “be defined and give infinity”. That, however, is not a promise of the C language itself. The standard’s rule for division leaves the behaviour undefined when the second operand is zero — for reals as much as for integers. Infinity appears in implementations that support IEC 60559 (that is, IEEE 754) semantics. In such an implementation, dividing a finite non-zero value by zero gives a signed infinity and raises the divide-by-zero exception, while $\frac{0}{0}$ falls on the NaN side. Today’s mainstream compilers on x86-64 and AArch64 behave that way, but it is not something the standard forces on every implementation.

50.2.1 Opening the bits directly

Having seen the layout in chapter 8, we now print the bits of real values and check them. To move a representation we use memcpy rather than a union — by chapter 48′s rule that is the safest passage for “moving a value”, and compilers mostly make the copy disappear.

/* Opening up a real number — how sign, exponent and fraction really sit in it.
   Type punning goes through memcpy, not a union (the rule of chapters 37 and 48). */
#include <inttypes.h>
#include <math.h>
#include <stdio.h>
#include <string.h>

static uint64_t bits_of(double d)
{
    uint64_t u;
    memcpy(&u, &d, sizeof u);          /* the representation moved as it is */
    return u;
}

static uint32_t bits_of_f(float f)
{
    uint32_t u;
    memcpy(&u, &f, sizeof u);
    return u;
}

/* double: 1 sign + 11 exponent + 52 fraction */
static void dump(const char *label, double d)
{
    uint64_t u = bits_of(d);
    unsigned sign = (unsigned)(u >> 63);
    unsigned expo = (unsigned)((u >> 52) & 0x7FFu);
    uint64_t frac = u & 0xFFFFFFFFFFFFFu;

    printf("%-12s %016" PRIx64 "  sign %u  exp %4u(=%+5d)  frac %013" PRIx64,
           label, u, sign, expo,
           expo == 0 ? -1022 : (int)expo - 1023, frac);

    if (expo == 0x7FF)      printf("  <- %s", frac ? "NaN" : "infinity");
    else if (expo == 0)     printf("  <- %s", frac ? "subnormal" : "zero");
    printf("\n");
}

int main(void)
{
    puts("-- the bits of a double (1 sign + 11 exponent + 52 fraction) --");
    dump("1.0", 1.0);
    dump("-1.0", -1.0);
    dump("0.5", 0.5);
    dump("2.0", 2.0);
    dump("0.1", 0.1);
    dump("0.3", 0.3);
    dump("0.1+0.2", 0.1 + 0.2);
    dump("0.0", 0.0);
    dump("-0.0", -0.0);
    dump("inf", INFINITY);
    dump("NaN", NAN);

    puts("\n-- 0.1 + 0.2 and 0.3 differ in their bits --");
    printf("0.1+0.2 = %.20f\n", 0.1 + 0.2);
    printf("0.3     = %.20f\n", 0.3);
    printf("the bits: %016" PRIx64 " vs %016" PRIx64 "  (the last bit)\n",
           bits_of(0.1 + 0.2), bits_of(0.3));

    puts("\n-- move one ULP from 1.0 and the last digit of the fraction rises by 1 --");
    double one = 1.0;
    double next = nextafter(1.0, 2.0);
    printf("1.0        %016" PRIx64 "\n", bits_of(one));
    printf("next value %016" PRIx64 "  difference %.17g\n", bits_of(next), next - one);
    printf("equal to DBL_EPSILON? %s\n",
           (next - one) == 0x1p-52 ? "yes" : "no");

    puts("\n-- a float uses the same structure, more narrowly (1 + 8 + 23) --");
    float f = 0.1f;
    uint32_t fu = bits_of_f(f);
    printf("0.1f       %08" PRIx32 "  sign %u  exp %3u(=%+4d)  frac %06" PRIx32 "\n",
           fu, fu >> 31, (fu >> 23) & 0xFFu, (int)((fu >> 23) & 0xFFu) - 127,
           fu & 0x7FFFFFu);
    printf("printing (double)0.1f gives %.17g — the trace of narrowing to float\n",
           (double)f);

    puts("\n-- below the smallest normal number come the subnormals --");
    double small = 0x1p-1022;          /* the smallest normal number */
    dump("2^-1022", small);
    dump("half", small / 2);           /* subnormal */
    dump("2^-1074", 0x1p-1074);        /* the smallest subnormal */
    dump("half again", 0x1p-1074 / 2); /* it sinks to zero */
    return 0;
}

-- the bits of a double (1 sign + 11 exponent + 52 fraction) --
1.0          3ff0000000000000  sign 0  exp 1023(=   +0)  frac 0000000000000
-1.0         bff0000000000000  sign 1  exp 1023(=   +0)  frac 0000000000000
0.5          3fe0000000000000  sign 0  exp 1022(=   -1)  frac 0000000000000
2.0          4000000000000000  sign 0  exp 1024(=   +1)  frac 0000000000000
0.1          3fb999999999999a  sign 0  exp 1019(=   -4)  frac 999999999999a
0.3          3fd3333333333333  sign 0  exp 1021(=   -2)  frac 3333333333333
0.1+0.2      3fd3333333333334  sign 0  exp 1021(=   -2)  frac 3333333333334
0.0          0000000000000000  sign 0  exp    0(=-1022)  frac 0000000000000  <- zero
-0.0         8000000000000000  sign 1  exp    0(=-1022)  frac 0000000000000  <- zero
inf          7ff0000000000000  sign 0  exp 2047(=+1024)  frac 0000000000000  <- infinity
NaN          7ff8000000000000  sign 0  exp 2047(=+1024)  frac 8000000000000  <- NaN

-- 0.1 + 0.2 and 0.3 differ in their bits --
0.1+0.2 = 0.30000000000000004441
0.3     = 0.29999999999999998890
the bits: 3fd3333333333334 vs 3fd3333333333333  (the last bit)

-- move one ULP from 1.0 and the last digit of the fraction rises by 1 --
1.0        3ff0000000000000
next value 3ff0000000000001  difference 2.2204460492503131e-16
equal to DBL_EPSILON? yes

-- a float uses the same structure, more narrowly (1 + 8 + 23) --
0.1f       3dcccccd  sign 0  exp 123(=  -4)  frac 4ccccd
printing (double)0.1f gives 0.10000000149011612 — the trace of narrowing to float

-- below the smallest normal number come the subnormals --
2^-1022      0010000000000000  sign 0  exp    1(=-1022)  frac 0000000000000
half         0008000000000000  sign 0  exp    0(=-1022)  frac 8000000000000  <- subnormal
2^-1074      0000000000000001  sign 0  exp    0(=-1022)  frac 0000000000001  <- subnormal
half again   0000000000000000  sign 0  exp    0(=-1022)  frac 0000000000000  <- zero

Five things from the output are worth pointing at.

First, 1.0 is remarkably tidy. The exponent field holds 1023 (the bias itself, so the actual exponent is 0) and the fraction is all zeros — the hidden bit alone makes $1.0\times 2^0$. 2.0 raises the exponent by one, 0.5 lowers it by one, and flipping the sign bit gives -1.0.

Second, 0.1 shows the cut mark of an unending fraction. Its fraction ends in 999999999999a, and that final a is the trace of rounding. It is chapter 8′s mathematics box — “it does not come out even in binary” — laid bare in bits.

Third, 0.1 + 0.2 and 0.3 differ by one last bit. The two bit patterns end ...3334 and ...3333, exactly one apart. That is why the == comparison is false, and why this chapter talks about tolerances.

Fourth, the identity of one ULP becomes visible. Adding the integer 1 to the bits of 1.0 gives the very next real number, and the difference is DBL_EPSILON ($2^{\left\{-52\right\}}$). “The smallest distinguishable difference near 1.0” turns out to be a single bit.

Fifth, the subnormals appear at the floor. Halve the smallest normal number and the exponent cannot go lower, so zeros begin to fill the front of the fraction instead — that state, with the exponent field all zeros, is a subnormal. Precision is given up little by little on the way down to zero, and when the last bit disappears the value becomes zero. This design, fading out instead of falling abruptly to zero, is called gradual underflow.

NaN has one famous property — it is not even equal to itself. If x != x is true then x is NaN, and that is the classic idiom for detecting NaN (today one uses isnan()). Being a value that breaks the basic property of the relation “equality”, NaN mixed into sorting or searching algorithms produces strange results — which is why checking for NaN at the boundary is the practice when handling real-number data.

We can handle the world of approximation in C. The next chapter is this part’s central subject — how a program deals with the fact that a computation can fail: the story of errors and contracts.