C / TYPES AND REPRESENTATION
Float precision, rounding and comparing with epsilon
Predict where a double loses precision, print the digits it really holds, and replace == with a tolerance that scales with the magnitude of the values.
What you will learn
- Print doubles with %.17g or %.17f to see the value actually stored
- Explain why 0.1 + 0.2 ends up one step above the double nearest 0.3
- Write fabs(a-b) <= ulps * DBL_EPSILON * larger magnitude instead of a == b
- Add an absolute floor for comparisons against zero, where relative slack vanishes
Understanding Float precision, rounding and comparing with epsilon
A double does not store 0.1. It stores a sign, a 53-bit binary significand and a power-of-two exponent, so the only values it holds exactly are integers scaled by powers of two. One tenth in binary is 0.0001100110011... repeating forever, just as one third has no finite decimal form, so the compiler substitutes the nearest representable double, 0.10000000000000000555. Every arithmetic result is rounded the same way, to the nearest double, with exact ties going to the candidate whose last significand bit is zero.
Picture the doubles as rungs on a ladder that grows coarser as you climb. Between 1.0 and 2.0 the rungs sit DBL_EPSILON = 2^-52 apart, about 2.22e-16; between 2^53 and 2^54 they sit 2.0 apart, which is why 1e16 + 1.0 == 1e16 is true; far below 1.0 they crowd together. An operation computes the exact answer and then slides to the nearest rung, so the error it introduces is roughly proportional to the size of the number rather than a fixed quantity. That is what DBL_DIG = 15 means: the guarantee covers significant digits, not decimal places.
Testing a == b asks whether two values landed on the same rung, a question about the rounding history of the calculation and almost never the question you meant. Use fabs(a - b) <= ulps * DBL_EPSILON * (larger of fabs(a) and fabs(b)), where the small integer factor covers however many rounded operations produced the values. A hardcoded absolute epsilon such as 1e-9 fails in both directions: it accepts 1e-12 and 5e-13 as equal, and it rejects two neighbouring doubles near 1e16. The one case a relative test cannot handle is comparison against zero, where the scale collapses to nothing; there you need an absolute floor drawn from the meaning of the numbers, such as one nanometre in a geometry program.
<stdio.h>
<math.h>
<float.h>
/* Tolerance measured in double-steps (ulps) at the magnitude of the operands. */
static int nearly_equal(double a, double b, double ulps)
{
double diff = fabs(a - b);
double scale = fabs(a) > fabs(b) ? fabs(a) : fabs(b);
return diff <= ulps * DBL_EPSILON * scale;
}
int main(void)
{
double sum = 0.1 + 0.2;
printf("0.1 + 0.2 = %.17f\n", sum);
printf("0.3 = %.17f\n", 0.3);
printf("sum == 0.3 : %d\n", sum == 0.3);
printf("|sum - 0.3| : %e\n", fabs(sum - 0.3));
printf("DBL_EPSILON : %e\n", DBL_EPSILON);
printf("nearly_equal: %d\n", nearly_equal(sum, 0.3, 4.0));
return 0;
}
Doubles round every literal and every operation to the nearest value on a scale whose spacing grows with magnitude, so equality has to become a tolerance that grows the same way.
Worked examples
Rounding error building up in a loop
Adding 0.1 ten times does not reach 1.0, and the default printf format hides it.
<stdio.h>
int main(void)
{
double sum = 0.0;
int i;
for (i = 0; i < 10; i++)
sum += 0.1;
printf("default printf : %f\n", sum);
printf("all the digits : %.17f\n", sum);
printf("1.0 - sum : %e\n", 1.0 - sum);
printf("sum == 1.0 : %d\n", sum == 1.0);
return 0;
}
Example explained
Line 1Each += rounds the running total again, so the errors do not simply pile up in one direction; the tenth partial sum lands on 1 - 2^-53, the largest double below 1.0.
Line 2%f shows six digits after the point, which rounds the total to 1.000000 and makes the failing comparison look impossible to explain.
Line 31.0 - sum is exact because subtracting two nearby doubles is itself representable, and the leftover 1.110223e-16 is exactly one step at this magnitude.
Line 4sum == 1.0 gives 0, so a loop written as while (sum != 1.0) would never terminate.
Why a fixed epsilon breaks at large magnitudes
Two adjacent doubles near 1e16 differ by 2.0, so a 1e-9 tolerance calls them different.
<stdio.h>
<math.h>
<float.h>
int main(void)
{
double a = 1e16;
double b = 1e16 + 2.0; /* the next double after 1e16 */
printf("a = %.1f\n", a);
printf("b = %.1f\n", b);
printf("b - a = %.1f\n", b - a);
printf("1e-9 test : %d\n", fabs(a - b) <= 1e-9);
printf("relative : %d\n", fabs(a - b) <= 4 * DBL_EPSILON * fabs(a));
printf("1e16 + 1.0 == 1e16 : %d\n", 1e16 + 1.0 == 1e16);
return 0;
}
Example explained
Line 11e16 is above 2^53, so the 53-bit significand has no bit left for odd values: doubles step by 2.0 there and nothing exists between a and b.
Line 2The 1e-9 test reports 0 because the tolerance is far below the smallest nonzero difference two distinct doubles can have at this magnitude.
Line 34 * DBL_EPSILON * fabs(a) is about 8.9, a few steps of slack in local units, so the relative test reports 1.
Line 41e16 + 1.0 sits exactly halfway between the two neighbours, and round-half-to-even picks the one with the even significand, handing back 1e16 unchanged.
float versus double, and comparing across them
A float keeps far fewer bits, so 0.1f widened to double is not the double nearest 0.1.
<stdio.h>
<float.h>
int main(void)
{
float f = 0.1f;
double d = 0.1;
printf("0.1f stored : %.20f\n", f);
printf("0.1 stored : %.20f\n", d);
printf("FLT_EPSILON : %e\n", FLT_EPSILON);
printf("f == 0.1 : %d\n", f == 0.1);
printf("reliable digits : float %d, double %d\n", FLT_DIG, DBL_DIG);
return 0;
}
Example explained
Line 1A float has 24 significand bits, so 0.1f is off by about 1.5e-9, while the 53 bits of a double put the error near 5.6e-18.
Line 2printf widens the float argument to double for you, and widening copies the coarse value; it cannot recover bits that were never stored.
Line 3f == 0.1 compares the widened float against the double nearest 0.1, two different values, so a literal that looks identical in the source fails.
Line 4FLT_DIG and DBL_DIG are the decimal digits that survive a decimal to binary to decimal round trip: 6 for float, 15 for double.
Important notes
DBL_EPSILON is the gap between 1.0 and the next double, not the smallest positive double (that is DBL_MIN) and not a tolerance by itself; it becomes one only after you multiply it by the magnitude being compared.
Results can shift if the compiler keeps intermediates wider than double (x87 on 32-bit x86, FLT_EVAL_METHOD 2) or reassociates under -ffast-math, and every comparison involving NaN is false, tolerance tests included.
Common mistakes
Testing total == 1.0 after adding 0.1 ten times: the value prints as 1.000000 under %f, so a branch that never fires looks impossible to explain.
Reusing one #define EPS 1e-9 everywhere: it calls 1e-12 and 5e-13 equal while calling two adjacent doubles near 1e16 different, so the same test is both too loose and too strict.
Comparing a float against a double literal, as in 0.1f == 0.1: the float widens to 0.10000000149011612 and the test is false even though the two literals look the same in the source.
Try it yourself
Change, predict, then run
Write nearly_equal(a, b, ulps) using a relative tolerance and print its verdict for the pairs (0.1 + 0.2, 0.3), (1e16, 1e16 + 2.0) and (0.0, 1e-20). Say which pair shows why a purely relative test still needs an absolute floor.
Open the C workspaceCheck your understanding
A tolerance of 1e-9 reports that 1e16 and 1e16 + 2.0 are different, even though no double exists between them. Why?
- 1e16 cannot be stored exactly, so the subtraction returns a garbage value.
- 1e-9 is smaller than DBL_EPSILON, so any comparison against it is false.
- Consecutive doubles near 1e16 are 2.0 apart, so a tolerance of 1e-9 sits far below the local spacing and can only ever accept bit-identical values.
- The subtraction underflows because the difference is too small for a double to represent.
Show answer
A double keeps 53 significand bits, so the gap between neighbours scales with magnitude: about 2.2e-16 near 1.0, and exactly 2.0 once you pass 2^53. Any tolerance smaller than that gap accepts nothing but identical values. The first option is tempting but wrong, since 1e16 is 2^16 times 5^16 and is stored exactly, as is 1e16 + 2.0, and their difference comes out as exactly 2.0. DBL_EPSILON is 2.2e-16, far smaller than 1e-9, so the second option has the comparison backwards.