Arithmetic and Data Types - Accuracy of floating point arithmetic

Chapter chap4 section 5

As a final point concerning floating point numbers, it is worth remembering that not all numbers can be stored exactly using floating point representation, as the following program shows.

main()
{
	double	z,y;
	z = 1.0/3.0;	/* one third */
	y = 1.0 - z - z - z;	/* should be zero */
	printf("%20.18le\n",y);
}
producing the output
1.110223024625156540e-16
This, clearly, isn't the expected zero. You will probably get different results on different computers. Of course, it is no great surprise that 1/3 cannot be represented exactly, in conventional decimal notation 0.3333.... is only ever an approximation. However it is less widely realised that fractions that can be expressed exactly in decimal notation (such as 1/10) cannot be stored exactly using computer floating point formats. The reason is that computers use an underlying binary notation rather than a decimal notation and only fractions with powers of 2 in the denominator (such 376/1024) can be represented exactly.

The following example illustrates the effects of precision. The program repeatedly adds smaller and smaller numbers to a variable originally holding 50000.0.

First the example using double variables

main()
{
        double  x = 50000.0;
        double  y = 0.01;
        x += y; /* adds 0.01 */
        printf("x = %30.15lf\n",x);
        y /= 100;
        x += y; /* adds 0.0001 */
        printf("x = %30.15lf\n",x);
        y /= 100;
        x += y; /* adds 0.000001 */
        printf("x = %30.15lf\n",x);
        y /= 100;
        x += y; /* adds 0.00000001 */
        printf("x = %30.15lf\n",x);
        y /= 100;
        x += y; /* adds 0.0000000001 */
        printf("x = %30.15lf\n",x);
        y /= 100;
        x += y; /* adds 0.000000000001 */
        printf("x = %30.15lf\n",x);
}
producing the output
x =          50000.010000000002037
x =          50000.010099999999511
x =          50000.010100999999850
x =          50000.010101009997015
x =          50000.010101010098879
x =          50000.010101010098879
The results aren't quite what a mathematician would expect, we've got 15 signifcant digits of accuracy guaranteed by our compiler but we've displayed the results to 20 significant figures, the last 5 should, of course, be regarded with suspicion. Notice that the final addition has made no difference to the stored number.

Modifying the program by changing the variables to float and the output conversion from "lf" to "f" and recompiling gave the following results.

x =          50000.011718750000000
x =          50000.011718750000000
x =          50000.011718750000000
x =          50000.011718750000000
x =          50000.011718750000000
x =          50000.011718750000000
Detailed comment hardly seems necessary. However this simple example does underline the fact that care needs to be taken with any program that performs floating point calculations, for fuller details a textbook on numerical analysis should be consulted.