
Go to the first, previous, next, last section, table of contents.
IEEE floatingpoint arithmeticThis chapter describes functions for examining the representation of floating point numbers and controlling the floating point environment of your program. The functions described in this chapter are declared in the header file `gsl_ieee_utils.h'. Representation of floating point numbersThe IEEE Standard for Binary FloatingPoint Arithmetic defines binary formats for single and double precision numbers. Each number is composed of three parts: a sign bit (s), an exponent (E) and a fraction (f). The numerical value of the combination (s,E,f) is given by the following formula, (1)^s (1.fffff...) 2^E The sign bit is either zero or one. The exponent ranges from a minimum value E_min to a maximum value E_max depending on the precision. The exponent is converted to an unsigned number e, known as the biased exponent, for storage by adding a bias parameter, e = E + bias. The sequence fffff... represents the digits of the binary fraction f. The binary digits are stored in normalized form, by adjusting the exponent to give a leading digit of 1. Since the leading digit is always 1 for normalized numbers it is assumed implicitly and does not have to be stored. Numbers smaller than 2^(E_min) are be stored in denormalized form with a leading zero, (1)^s (0.fffff...) 2^(E_min) This allows gradual underflow down to 2^(E_min  p) for p bits of precision. A zero is encoded with the special exponent of 2^(E_min  1) and infinities with the exponent of 2^(E_max + 1). The format for single precision numbers uses 32 bits divided in the following way, seeeeeeeefffffffffffffffffffffff s = sign bit, 1 bit e = exponent, 8 bits (E_min=126, E_max=127, bias=127) f = fraction, 23 bits The format for double precision numbers uses 64 bits divided in the following way, seeeeeeeeeeeffffffffffffffffffffffffffffffffffffffffffffffffffff s = sign bit, 1 bit e = exponent, 11 bits (E_min=1022, E_max=1023, bias=1023) f = fraction, 52 bits It is often useful to be able to investigate the behavior of a calculation at the bitlevel and the library provides functions for printing the IEEE representations in a humanreadable form.
#include <stdio.h> #include <gsl/gsl_ieee_utils.h> int main (void) { float f = 1.0/3.0; double d = 1.0/3.0; double fd = f; /* promote from float to double */ printf(" f="); gsl_ieee_printf_float(&f); printf("\n"); printf("fd="); gsl_ieee_printf_double(&fd); printf("\n"); printf(" d="); gsl_ieee_printf_double(&d); printf("\n"); return 0; } The binary representation of 1/3 is 0.01010101... . The output below shows that the IEEE format normalizes this fraction to give a leading digit of 1, f= 1.01010101010101010101011*2^2 fd= 1.0101010101010101010101100000000000000000000000000000*2^2 d= 1.0101010101010101010101010101010101010101010101010101*2^2 The output also shows that a singleprecision number is promoted to doubleprecision by adding zeros in the binary representation.
Setting up your IEEE environmentThe IEEE standard defines several modes for controlling the behavior of floating point operations. These modes specify the important properties of computer arithmetic: the direction used for rounding (e.g. whether numbers should be rounded up, down or to the nearest number), the rounding precision and how the program should handle arithmetic exceptions, such as division by zero.
Many of these features can now be controlled via standard functions such
as
e = 1 + 1/2! + 1/3! + 1/4! + ... = 2.71828182846...
#include <math.h> #include <stdio.h> #include <gsl/gsl_ieee_utils.h> int main (void) { double x = 1, oldsum = 0, sum = 0; int i = 0; gsl_ieee_env_setup (); /* read GSL_IEEE_MODE */ do { i++; oldsum = sum; sum += x; x = x / i; printf("i=%2d sum=%.18f error=%g\n", i, sum, sum  M_E); if (i > 30) break; } while (sum != oldsum); return 0; }
Here are the results of running the program in GSL_IEEE_MODE="roundtonearest" ./a.out i= 1 sum=1.000000000000000000 error=1.71828 i= 2 sum=2.000000000000000000 error=0.718282 .... i=18 sum=2.718281828459045535 error=4.44089e16 i=19 sum=2.718281828459045535 error=4.44089e16
After nineteen terms the sum converges to within
4 \times 10^16 of the correct value.
If we now change the rounding mode to
GSL_IEEE_MODE="rounddown" ./a.out i= 1 sum=1.000000000000000000 error=1.71828 .... i=19 sum=2.718281828459041094 error=3.9968e15
The result is about
4 \times 10^15
below the correct value, an order of magnitude worse than the result
obtained in the
If we change to rounding mode to
Finally we can see the effect of computing the sum using
singleprecision rounding, in the default GSL_IEEE_MODE="singleprecision" ./a.out .... i=12 sum=2.718281984329223633 error=1.5587e07 with an error of O(10^7), which corresponds to single precision accuracy (about 1 part in 10^7). Continuing the iterations further does not decrease the error because all the subsequent results are rounded to the same value. References and Further ReadingThe reference for the IEEE standard is,
A more pedagogical introduction to the standard can be found in the paper "What Every Computer Scientist Should Know About FloatingPoint Arithmetic".
Go to the first, previous, next, last section, table of contents. 