Types and how data lives in memory · about 50 min

Floating point: how real numbers are approximated

Why this matters

Measurements, physics, graphics, statistics, audio, machine learning: much of computing deals with numbers that have fractional parts. C represents them with floating-point types, which are fast and have an enormous range, but are approximations. A surprising number of bugs come from forgetting that: totals that are off by a cent, loops that never reach their end value, equality tests that fail for numbers that "are" equal. This lesson shows what is really stored, so the surprises become predictable.

float and double

C has three floating-point types:

type	size	significant decimal digits	range (approx.)	printf
`float`	4 bytes	about 6-7	10^-38 to 10^38	`%f`
`double`	8 bytes	about 15-16	10^-308 to 10^308	`%f`
`long double`	16 bytes on x86-64 Linux (80 bits used)	about 18-19	10^-4932 to 10^4932	`%Lf`

Use double by default. float saves memory in huge arrays (graphics, ML) but runs out of precision quickly. A literal with a decimal point or exponent, such as 3.14 or 2.5e-3, is a double; add f for a float literal: 3.14f. Both float and double print with %f, %e or %g, because a float passed to printf is automatically converted to double.

Scientific notation in binary

A floating-point number is stored like scientific notation, but in base 2. In decimal you might write 6.02 * 10^23: a sign, some significant digits (the significand or mantissa), and an exponent saying where the point goes. A double (IEEE 754 binary64, used by virtually every computer) stores:

  1 bit    11 bits          52 bits
  [sign]   [exponent]       [fraction]
    0      10000000000      1001001000011111101101010100010001000010110100011000

  value = (-1)^sign * 1.fraction (in binary) * 2^(exponent - 1023)

That example is pi, approximately. The leading 1. is implied (not stored), so a double has 53 significant bits. A float uses 1 + 8 + 23 bits, so 24 significant bits.

The key insight: the significand has a fixed number of bits, and the exponent says where the binary point goes. That means:

Precision is relative. Near 1.0, doubles are about 2 * 10^-16 apart. Near 10^16, they are 2 apart: 1e16 + 1 is stored as exactly 1e16.
Numbers that need infinitely many binary digits are rounded to 53 bits.

Why 0.1 is not 0.1

In decimal, 1/3 = 0.3333... never ends. In binary, one tenth has the same problem: 0.1 = 0.000110011001100110011... repeating forever. So the double nearest to 0.1 is not exactly 0.1:

#include <stdio.h>

int main(void)
{
    printf("%.20f\n", 0.1);
    double a = 0.1 + 0.2;
    printf("%f\n", a);
    printf("%.17g\n", a);
    printf("%d\n", a == 0.3);
    return 0;
}

0.10000000000000000555
0.300000
0.30000000000000004
0

%f rounds to six decimals, which hides the error. %.17g shows enough digits to identify the exact double. Each of 0.1 and 0.2 carries a tiny rounding error; their sum lands on a different double than the literal 0.3. Errors also accumulate:

double sum = 0.0;
/* add 0.1 ten times (you will learn loops soon) */

gives 0.99999999999999989, not 1. Fractions whose denominator is a power of two (0.5, 0.25, 0.375) are exact; most decimal fractions are not.

Precision of float versus double, computing 1/3:

float:  0.33333334326744079590
double: 0.33333333333333331483

The float is right to about 7 digits, the double to about 16.

Comparing floating-point values

Because of rounding, do not use == on computed floating-point values. Instead ask whether they are close enough. An absolute tolerance (fabs(a - b) < 1e-9) fails for very large numbers (where neighbouring doubles are further apart than that) and is meaningless for very small ones. A relative tolerance scales with the magnitude:

#include <math.h>
/* equal if they differ by at most one part in a billion */
fabs(a - b) <= 1e-9 * fmax(fabs(a), fabs(b))

The right tolerance depends on the problem; there is no universal value.

Special values: infinity and NaN

IEEE 754 reserves exponent patterns for special values. Dividing a nonzero double by zero gives infinity, and 0.0 / 0.0 gives NaN ("not a number"):

printf("%g %g %g\n", 1.0 / 0.0, -1.0 / 0.0, 0.0 / 0.0);

inf -inf -nan

NaN is contagious (any arithmetic with NaN gives NaN) and compares unequal to everything, including itself: if x is NaN, x == x is 0. <math.h> has isnan(x) and isinf(x) to test for them.

Warning: this is only for floating point. Integer division by zero is undefined behavior and on Linux usually kills the program with "Floating point exception" (the name is historical; it is an integer fault).

Overflowing a double gives infinity rather than undefined behavior; the largest double is about 1.8 * 10^308 (DBL_MAX in <float.h>, which also has DBL_EPSILON, FLT_DIG and friends).

The math library

<math.h> declares the usual functions, all taking and returning double: sqrt, pow, exp, log, log10, sin, cos, tan, atan2, fabs, floor, ceil, round, fmod, fmax, fmin, and more. Remember to link with -lm (the track's build command does).

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

int main(void)
{
    printf("%.2f\n", 1000 * pow(1.05, 10));
    printf("%.4f\n", sqrt(2.0));
    printf("%.1f %.1f %.1f\n", floor(2.7), ceil(2.1), round(2.5));
    return 0;
}

1628.89
1.4142
2.0 3.0 3.0

Common mistakes

Using == on computed doubles. Use a tolerance, as above.
Using floating point for money. $0.10 cannot be stored exactly. Store whole cents in an integer type (long long cents) and format on output.
Expecting M_PI to exist. It is not part of standard C, so with -std=c17 you get Define your own (const double pi = 3.14159265358979323846;) or use acos(-1.0).
```
error: 'M_PI' undeclared (first use in this function)
```
Integer division where you wanted a fraction. 1 / 3 is 0, because both operands are integers; write 1.0 / 3. The next module covers this properly.
Printing a double with %d (or an int with %f): garbage output and a -Wformat warning.
Using float without a reason. Its 7 digits run out fast: a float cannot even count past 16,777,216 by ones.

Recap

float (4 bytes, about 7 digits) and double (8 bytes, about 16 digits) store a sign, an exponent and a significand in base 2 (IEEE 754). Prefer double.
Most decimal fractions, including 0.1, are rounded; errors accumulate.
Precision is relative: the gap between neighbouring values grows with their size.
Compare with a relative tolerance, not ==.
Floating-point division by zero gives inf or nan; integer division by zero is undefined.
<math.h> (link with -lm) has the math functions; <float.h> has the limits.

Going further: the classic article "What Every Computer Scientist Should Know About Floating-Point Arithmetic" (David Goldberg) goes deep. For a quick experiment, print 0.1 + 0.2 with %a, which shows the exact binary value in hexadecimal notation.

Exercises

Exercise 1 · quiz

Does 0.1 + 0.2 equal 0.3?

What does this print?

double a = 0.1 + 0.2;
printf("%d\n", a == 0.3);

(== compares two values and gives 1 if they are equal, 0 if not.)

1
0
0.3
It does not compile

Answer

0.1, 0.2 and 0.3 cannot be represented exactly in binary. The nearest doubles to 0.1 and 0.2 add up to 0.30000000000000004, while the nearest double to 0.3 is 0.29999999999999998890. They differ in the last bit, so == says no. Never compare computed floating-point values with ==; compare with a tolerance.

Exercise 2 · quiz

Big integers in a float

A float has a 24-bit significand. What does this print?

float f = 16777217;     /* 2 to the 24th, plus 1 */
printf("%.1f\n", f);

16777217.0
16777216.0
16777218.0
inf

Answer

16777216.0

A float can represent every integer up to 2^24 = 16777216 exactly, but above that the gap between neighbouring floats is 2, so 16777217 is rounded to the nearest representable value, 16777216. For double the limit is 2^53. This is why money, counters and IDs belong in integer types.

Exercise 3 · write code

Compound interest

1000 dollars are invested at 5 percent interest per year, compounded yearly. After n years the balance is 1000 * 1.05 to the power n. Use pow from <math.h> to print the balance after 10 years and after 30 years, rounded to cents, exactly:

after 10 years: 1628.89
after 30 years: 4321.94

File: main.c

Starting point

#include <stdio.h>

int main(void)
{
    double principal = 1000.0;
    double rate = 0.05;
    /* print the two lines using pow(base, exponent) */
    return 0;
}

Example

expected output

after 10 years: 1628.89
after 30 years: 4321.94

Hint 1

Include <math.h>. pow(x, y) returns x to the power y as a double.

Hint 2

%.2f prints with two decimals, rounding correctly.

Hint 3

printf("after 10 years: %.2f\n", principal * pow(1.0 + rate, 10));

Reference solution

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

int main(void)
{
    double principal = 1000.0;
    double rate = 0.05;
    printf("after 10 years: %.2f\n", principal * pow(1.0 + rate, 10));
    printf("after 30 years: %.2f\n", principal * pow(1.0 + rate, 30));
    return 0;
}

Exercise 4 · write code

Compare with a tolerance

Write a function

int nearly_equal(double a, double b)

that returns 1 when a and b are equal within a relative tolerance of 1e-9, and 0 otherwise. "Relative" means the allowed difference scales with the size of the numbers:

|a - b| <= 1e-9 * max(|a|, |b|)

Use fabs (absolute value of a double) and fmax (larger of two doubles) from <math.h>. A comparison such as x <= y is itself an expression whose value is 1 (true) or 0 (false), so the whole function can be a single return statement.

The check calls it with pairs such as 0.1 + 0.2 and 0.3 (equal), 1e20 and 1e20 + 1e5 (equal: the difference is tiny relative to 1e20), 1.0 and 1.001 (not equal), and 0.0 and 0.0 (equal).

File: main.c

Starting point

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

int nearly_equal(double a, double b)
{
    return a == b;
}

int main(void)
{
    printf("%d\n", nearly_equal(0.1 + 0.2, 0.3));
    return 0;
}

Hint 1

The difference is fabs(a - b).

Hint 2

The allowed difference is 1e-9 * fmax(fabs(a), fabs(b)).

Hint 3

Return the comparison directly: return fabs(a - b) <= 1e-9 * fmax(fabs(a), fabs(b));

Reference solution

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

int nearly_equal(double a, double b)
{
    return fabs(a - b) <= 1e-9 * fmax(fabs(a), fabs(b));
}

int main(void)
{
    printf("%d\n", nearly_equal(0.1 + 0.2, 0.3));
    return 0;
}