The 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 fpsetround, which should be used whenever they are available.
Unfortunately in the past there has been no universal API for
controlling their behavior -- each system has had its own way of
accessing them. For example, the Linux kernel provides the function
__setfpucw (set-fpu-control-word) to set IEEE modes, while
HP-UX and Solaris use the functions fpsetround and
fpsetmask. To help you write portable programs GSL allows you to
specify modes in a platform-independent way using the environment
variable GSL_IEEE_MODE. The library then takes care of all the
necessary machine-specific initializations for you when you call the
function gsl_ieee_env_setup.
GSL_IEEE_MODE and
attempts to set up the corresponding specified IEEE modes. The
environment variable should be a list of keywords, separated by
commas, like this,
GSL_IEEE_MODE = "keyword,keyword,..."
where keyword is one of the following mode-names,
single-precision
double-precision
extended-precision
round-to-nearest
round-down
round-up
round-to-zero
mask-all
mask-invalid
mask-denormalized
mask-division-by-zero
mask-overflow
mask-underflow
trap-inexact
trap-common
If GSL_IEEE_MODE is empty or undefined then the function returns
immediately and no attempt is made to change the system's IEEE
mode. When the modes from GSL_IEEE_MODE are turned on the
function prints a short message showing the new settings to remind you
that the results of the program will be affected.
If the requested modes are not supported by the platform being used then
the function calls the error handler and returns an error code of
GSL_EUNSUP.
The following combination of modes is convenient for many purposes,
GSL_IEEE_MODE="double-precision,"\
"mask-underflow,"\
"mask-denormalized"
This choice ignores any errors relating to small numbers (either denormalized, or underflowing to zero) but traps overflows, division by zero and invalid operations.
#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 round-to-nearest
mode. This is the IEEE default so it isn't really necessary to specify
it here,
GSL_IEEE_MODE="round-to-nearest" ./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.44089e-16 i=19 sum=2.718281828459045535 error=4.44089e-16
After nineteen terms the sum converges to within @c{$4 \times 10^{-16}$}
@math{4 \times 10^-16} of the correct value.
If we now change the rounding mode to
round-down the final result is less accurate,
GSL_IEEE_MODE="round-down" ./a.out i= 1 sum=1.000000000000000000 error=-1.71828 .... i=19 sum=2.718281828459041094 error=-3.9968e-15
The result is about
@math{4 \times 10^-15}
below the correct value, an order of magnitude worse than the result
obtained in the round-to-nearest mode.
If we change to rounding mode to round-up then the series no
longer converges (the reason is that when we add each term to the sum
the final result is always rounded up. This is guaranteed to increase the sum
by at least one tick on each iteration). To avoid this problem we would
need to use a safer converge criterion, such as while (fabs(sum -
oldsum) > epsilon), with a suitably chosen value of epsilon.
Finally we can see the effect of computing the sum using
single-precision rounding, in the default round-to-nearest
mode. In this case the program thinks it is still using double precision
numbers but the CPU rounds the result of each floating point operation
to single-precision accuracy. This simulates the effect of writing the
program using single-precision float variables instead of
double variables. The iteration stops after about half the number
of iterations and the final result is much less accurate,
GSL_IEEE_MODE="single-precision" ./a.out .... i=12 sum=2.718281984329223633 error=1.5587e-07
with an error of @math{O(10^-7)}, which corresponds to single precision accuracy (about 1 part in @math{10^7}). Continuing the iterations further does not decrease the error because all the subsequent results are rounded to the same value.