## Setting up your IEEE environment

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`.

Function: void gsl_ieee_env_setup ()
  This function reads the environment variable `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.

To demonstrate the effects of different rounding modes consider the following program which computes e, the base of natural logarithms, by summing a rapidly-decreasing series,

    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 `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 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 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 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.