This is part of the GNU C Intro and Reference Manual and covered by its license.

Floating Point in Depth
***********************

Floating-Point Representations
==============================

Storing numbers as floating point allows representation of numbers with fractional values, in a range larger than that of hardware integers.  A floating-point number consists of a sign bit, a significand (also called the mantissa), and a power of a fixed base.  GNU C uses the floating-point representations specified by the IEEE 754-2008 Standard for Floating-Point Arithmetic.

The IEEE 754-2008 specification defines basic binary floating-point formats of five different sizes: 16-bit, 32-bit, 64-bit, 128-bit, and 256-bit.  The formats of 32, 64, and 128 bits are used for the standard C types float, double, and long double. GNU C supports the 16-bit floating point type _Float16 on some platforms, but does not support the 256-bit floating point type.

Each of the formats encodes the floating-point number as a sign bit. After this comes an exponent that specifies a power of 2 (with a fixed offset).  Then comes the significand.

The first bit of the significand, before the binary point, is always 1, so there is no need to store it in memory.  It is called the hidden bit because it doesn't appear in the floating-point number as used in the computer itself.

All of those floating-point formats are sign-magnitude representations, so +0 and -0 are different values.

Besides the IEEE 754 format 128-bit float, GNU C also offers a format consisting of a pair of 64-bit floating point numbers.  This lacks the full exponent range of the IEEE 128-bit format, but is useful when the underlying hardware platform does not support that.

Floating-Point Type Specifications
==================================

The standard library header file float.h defines a number of constants that describe the platform's implementation of floating-point types float, double and long double.  They include:

FLT_MIN
DBL_MIN
LDBL_MIN
     Defines the minimum normalized positive floating-point values that can be represented with the type.

FLT_HAS_SUBNORM
DBL_HAS_SUBNORM
LDBL_HAS_SUBNORM
     Defines if the floating-point type supports subnormal (or ``denormalized'') numbers or not (see subnormal numbers).

FLT_TRUE_MIN
DBL_TRUE_MIN
LDBL_TRUE_MIN
     Defines the minimum positive values (including subnormal values) that can be represented with the type.

FLT_MAX
DBL_MAX
LDBL_MAX
     Defines the largest values that can be represented with the type.

FLT_DECIMAL_DIG
DBL_DECIMAL_DIG
LDBL_DECIMAL_DIG
     Defines the number of decimal digits n such that any floating-point number that can be represented in the type can be rounded to a floating-point number with n decimal digits, and back again, without losing any precision of the value.

Special Floating-Point Values
=============================

IEEE floating point provides for special values that are not ordinary numbers.

infinities
     +Infinity and -Infinity are two different infinite values, one positive and one negative.  These result from operations such as 1 / 0, Infinity + Infinity, Infinity * Infinity, and Infinity + FINITE, and also from a result that is finite, but larger than the most positive possible value or smaller than the most negative possible value.

     See Handling Infinity, for more about working with infinities.

NaNs (not a number)
     There are two special values, called Not-a-Number (NaN): a quiet NaN (QNaN), and a signaling NaN (SNaN).

     A QNaN is produced by operations for which the value is undefined in real arithmetic, such as 0 / 0, sqrt (-1), Infinity - Infinity, and any basic operation in which an operand is a QNaN.

     The signaling NaN is intended for initializing otherwise-unassigned storage, and the goal is that unlike a QNaN, an SNaN does cause an interrupt that can be caught by a software handler, diagnosed, and reported.  In practice, little use has been made of signaling NaNs, because the most common CPUs in desktop and portable computers fail to implement the full IEEE 754 Standard, and supply only one kind of NaN, the quiet one.  Also, programming-language standards have taken decades to catch up to the IEEE 754 standard, and implementations of those language standards make an additional delay before programmers become willing to use these features.

     To enable support for signaling NaNs, use the GCC command-line option -fsignaling-nans, but this is an experimental feature and may not work as expected in every situation.

     A NaN has a sign bit, but its value means nothing.

     See Handling NaN, for more about working with NaNs.

subnormal numbers
     It can happen that a computed floating-point value is too small to represent, such as when two tiny numbers are multiplied.  The result is then said to underflow.  The traditional behavior before the IEEE 754 Standard was to use zero as the result, and possibly to report the underflow in some sort of program output.

     The IEEE 754 Standard is vague about whether rounding happens before detection of floating underflow and overflow, or after, and CPU designers may choose either.

     However, the Standard does something unusual compared to earlier designs, and that is that when the result is smaller than the smallest normalized representable value (i.e., one in which the leading significand bit is 1), the normalization requirement is relaxed, leading zero bits are permitted, and precision is gradually lost until there are no more bits in the significand.  That phenomenon is called gradual underflow, and it serves important numerical purposes, although it does reduce the precision of the final result.  Some floating-point designs allow you to choose at compile time, or even at run time, whether underflows are gradual, or are flushed abruptly to zero.  Numbers that have entered the region of gradual underflow are called subnormal.

     You can use the library functions fesetround and fegetround to set and get the rounding mode.  Rounding modes are defined (if supported by the platform) in fenv.h as: FE_UPWARD to round toward positive infinity; FE_DOWNWARD to round toward negative infinity; FE_TOWARDZERO to round toward zero; and FE_TONEAREST to round to the nearest representable value, the default mode.  It is best to use FE_TONEAREST except when there is a special need for some other mode.

Invalid Optimizations
=====================

Signed zeros, Infinity, and NaN invalidate some optimizations by programmers and compilers that might otherwise have seemed obvious:

  * x + 0 and x - 0 are not the same as x when x is zero, because the result depends on the rounding rule. See Rounding, for more about rounding rules.

  * x * 0.0 is not the same as 0.0 when x is Infinity, a NaN, or negative zero.

  * x / x is not the same as 1.0 when x is Infinity, a NaN, or zero.

  * (x - y) is not the same as -(y - x) because when the operands are finite and equal, one evaluates to +0 and the other to -0.

  * x - x is not the same as 0.0 when X is Infinity or a NaN.

  * x == x and x != x are not equivalent to 1 and 0 when X is a NaN.

  * x < y and isless (x, y) are not equivalent, because the first sets a sticky exception flag (see Exception Flags) when an operand is a NaN, whereas the second does not affect that flag.  The same holds for the other isxxx functions that are companions to relational operators.  See FP Comparison Functions (The GNU C Library Reference Manual).

The -funsafe-math-optimizations option enables these optimizations.

Floating Arithmetic Exception Flags
===================================

Sticky exception flags record the occurrence of particular conditions: once set, they remain set until the program explicitly clears them.

The conditions include invalid operand, division-by_zero, inexact result (i.e., one that required rounding), underflow, and overflow.  Some extended floating-point designs offer several additional exception flags.  The functions feclearexcept, feraiseexcept, fetestexcept, fegetexceptflags, and fesetexceptflags provide a standardized interface to those flags.  See Status bit operations (The GNU C Library Reference Manual).

One important use of those flags is to do a computation that is normally expected to be exact in floating-point arithmetic, but occasionally might not be, in which case, corrective action is needed.  You can clear the inexact result flag with a call to feclearexcept (FE_INEXACT), do the computation, and then test the flag with fetestexcept (FE_INEXACT); the result of that call is 0 if the flag is not set (there was no rounding), and 1 when there was rounding (which, we presume, implies the program has to correct for that).

Exact Floating-Point Arithmetic
===============================

As long as the numbers are exactly representable (fractions whose denominator is a power of 2), and intermediate results do not require rounding, then floating-point arithmetic is exact.  It is easy to predict how many digits are needed for the results of arithmetic operations:

  * addition and subtraction of two N-digit values with the same exponent require at most N + 1 digits, but when the exponents differ, many more digits may be needed;

  * multiplication of two N-digit values requires exactly 2 N digits;

  * although integer division produces a quotient and a remainder of no more than N-digits, floating-point remainder and square root may require an unbounded number of digits, and the quotient can need many more digits than can be stored.

Whenever a result requires more than N digits, rounding is needed.

Rounding
========

When floating-point arithmetic produces a result that can't fit exactly in the significand of the type that's in use, it has to round the value.  The basic arithmetic operations---addition, subtraction, multiplication, division, and square root---always produce a result that is equivalent to the exact, possibly infinite-precision result rounded to storage precision according to the current rounding rule.

Rounding sets the FE_INEXACT exception flag (see Exception Flags).  This enables programs to determine that rounding has occurred.

Rounding consists of adjusting the exponent to bring the significand back to the required base-point alignment, then applying the current rounding rule to squeeze the significand into the fixed available size.

The current rule is selected at run time from four options.  Here they are:

  * round-to-nearest, with ties rounded to an even integer;

  * round-up, towards +Infinity;

  * round-down, towards -Infinity;

  * round-towards-zero.

Under those four rounding rules, a decimal value -1.2345 that is to be rounded to a four-digit result would become -1.234, -1.234, -1.235, and -1.234, respectively.

The default rounding rule is round-to-nearest, because that has the least bias, and produces the lowest average error.  When the true result lies exactly halfway between two representable machine numbers, the result is rounded to the one that ends with an even digit.

The round-towards-zero rule was common on many early computer designs, because it is the easiest to implement: it just requires silent truncation of all extra bits.

The two other rules, round-up and round-down, are essential for implementing interval arithmetic, whereby each arithmetic operation produces lower and upper bounds that are guaranteed to enclose the exact result.

See Rounding Control, for details on getting and setting the current rounding mode.

Rounding Issues
===============

The default IEEE 754 rounding mode minimizes errors, and most normal computations should not suffer any serious accumulation of errors from rounding.

Of course, you can contrive examples where that is not so.  Here is one: iterate a square root, then attempt to recover the original value by repeated squaring.

     #include <stdio.h>
     #include <math.h>

     int main (void)
     {
       double x = 100.0;
       double y;
       int n, k;
       for (n = 10; n <= 100; n += 10)
         {
           y = x;
           for (k = 0; k < n; ++k) y = sqrt (y);
           for (k = 0; k < n; ++k) y *= y;
           printf ("n = %3d; x = %.0f\ty = %.6f\n", n, x, y);
         }
       return 0;
     }

Here is the output:

     n =  10; x = 100        y = 100.000000
     n =  20; x = 100        y = 100.000000
     n =  30; x = 100        y = 99.999977
     n =  40; x = 100        y = 99.981025
     n =  50; x = 100        y = 90.017127
     n =  60; x = 100        y = 1.000000
     n =  70; x = 100        y = 1.000000
     n =  80; x = 100        y = 1.000000
     n =  90; x = 100        y = 1.000000
     n = 100; x = 100        y = 1.000000

After 50 iterations, y has barely one correct digit, and soon after, there are no correct digits.

Significance Loss
=================

A much more serious source of error in floating-point computation is significance loss from subtraction of nearly equal values.  This means that the number of bits in the significand of the result is fewer than the size of the value would permit.  If the values being subtracted are close enough, but still not equal, a single subtraction can wipe out all correct digits, possibly contaminating all future computations.

Floating-point calculations can sometimes be carefully designed so that significance loss is not possible, such as summing a series where all terms have the same sign.  For example, the Taylor series expansions of the trigonometric and hyperbolic sines have terms of identical magnitude, of the general form X**(2*N + 1) / (2*N + 1)!.  However, those in the trigonometric sine series alternate in sign, while those in the hyperbolic sine series are all positive.  Here is the output of two small programs that sum K terms of the series for sin (X), and compare the computed sums with known-to-be-accurate library functions:

     x = 10      k = 51
     s (x)   = -0.544_021_110_889_270
     sin (x) = -0.544_021_110_889_370

     x = 20      k = 81
     s (x)   = 0.912_945_250_749_573
     sin (x) = 0.912_945_250_727_628

     x = 30      k = 109
     s (x)   = -0.987_813_746_058_855
     sin (x) = -0.988_031_624_092_862

     x = 40      k = 137
     s (x)   = 0.617_400_430_980_474
     sin (x) = 0.745_113_160_479_349

     x = 50      k = 159
     s (x)   = 57_105.187_673_745_720_532
     sin (x) = -0.262_374_853_703_929

     // sinh(x) series summation with positive signs
     // with k terms needed to converge to machine precision

     x = 10      k = 47
     t (x)    = 1.101_323_287_470_340e+04
     sinh (x) = 1.101_323_287_470_339e+04

     x = 20      k = 69
     t (x)    = 2.425_825_977_048_951e+08
     sinh (x) = 2.425_825_977_048_951e+08

     x = 30      k = 87
     t (x)    = 5.343_237_290_762_229e+12
     sinh (x) = 5.343_237_290_762_231e+12

     x = 40      k = 105
     t (x)    = 1.176_926_334_185_100e+17
     sinh (x) = 1.176_926_334_185_100e+17

     x = 50      k = 121
     t (x)    = 2.592_352_764_293_534e+21
     sinh (x) = 2.592_352_764_293_536e+21

We have added underscores to the numbers to enhance readability.

The sinh (X) series with positive terms can be summed to high accuracy.  By contrast, the series for sin (X) suffers increasing significance loss, so that when X = 30 only two correct digits remain.  Soon after, all digits are wrong, and the answers are complete nonsense.

An important skill in numerical programming is to recognize when significance loss is likely to contaminate a computation, and revise the algorithm to reduce this problem.  Sometimes, the only practical way to do so is to compute in higher intermediate precision, which is why the extended types like long double are important.

Fused Multiply-Add
==================

In 1990, when IBM introduced the POWER architecture, the CPU provided a previously unknown instruction, the fused multiply-add (FMA).  It computes the value x * y + z with an exact double-length product, followed by an addition with a single rounding.  Numerical computation often needs pairs of multiply and add operations, for which the FMA is well-suited.

On the POWER architecture, there are two dedicated registers that hold permanent values of 0.0 and 1.0, and the normal multiply and add instructions are just wrappers around the FMA that compute x * y + 0.0 and x * 1.0 + z, respectively.

In the early days, it appeared that the main benefit of the FMA was getting two floating-point operations for the price of one, almost doubling the performance of some algorithms.  However, numerical analysts have since shown numerous uses of the FMA for significantly enhancing accuracy.  We discuss one of the most important ones in the next section.

A few other architectures have since included the FMA, and most provide variants for the related operations x * y - z (FMS), -x * y + z (FNMA), and -x * y - z (FNMS).

The functions fmaf, fma, and fmal implement fused multiply-add for the float, double, and long double data types.  Correct implementation of the FMA in software is difficult, and some systems that appear to provide those functions do not satisfy the single-rounding requirement.  That situation should change as more programmers use the FMA operation, and more CPUs provide FMA in hardware.

Use the -ffp-contract=fast option to allow generation of FMA instructions, or -ffp-contract=off to disallow it.

Error Recovery
==============

When two numbers are combined by one of the four basic operations, the result often requires rounding to storage precision.  For accurate computation, one would like to be able to recover that rounding error.  With historical floating-point designs, it was difficult to do so portably, but now that IEEE 754 arithmetic is almost universal, the job is much easier.

For addition with the default round-to-nearest rounding mode, we can determine the error in a sum like this:

     volatile double err, sum, tmp, x, y;

     if (fabs (x) >= fabs (y))
       {
         sum = x + y;
         tmp = sum - x;
         err = y - tmp;
       }
     else /* fabs (x) < fabs (y) */
       {
         sum = x + y;
         tmp = sum - y;
         err = x - tmp;
       }

Now, x + y is exactly represented by sum + err. This basic operation, which has come to be called twosum in the numerical-analysis literature, is the first key to tracking, and accounting for, rounding error.

To determine the error in subtraction, just swap the + and - operators.

We used the volatile qualifier (see volatile) in the declaration of the variables, which forces the compiler to store and retrieve them from memory, and prevents the compiler from optimizing err = y - ((x + y) - x) into err = 0.

For multiplication, we can compute the rounding error without magnitude tests with the FMA operation (see Fused Multiply-Add), like this:

     volatile double err, prod, x, y;
     prod = x * y;                /* rounded product */
     err  = fma (x, y, -prod);    /* exact product = PROD + ERR */

For addition, subtraction, and multiplication, we can represent the exact result with the notional sum of two values.  However, the exact result of division, remainder, or square root potentially requires an infinite number of digits, so we can at best approximate it. Nevertheless, we can compute an error term that is close to the true error: it is just that error value, rounded to machine precision.

For division, you can approximate x / y with quo + err like this:

     volatile double err, quo, x, y;
     quo = x / y;
     err = fma (-quo, y, x) / y;

For square root, we can approximate sqrt (x) with root + err like this:

     volatile double err, root, x;
     root = sqrt (x);
     err = fma (-root, root, x) / (root + root);

With the reliable and predictable floating-point design provided by IEEE 754 arithmetic, we now have the tools we need to track errors in the five basic floating-point operations, and we can effectively simulate computing in twice working precision, which is sometimes sufficient to remove almost all traces of arithmetic errors.

Exact Floating-Point Constants
==============================

One of the frustrations that numerical programmers have suffered with since the dawn of digital computers is the inability to precisely specify numbers in their programs.  On the early decimal machines, that was not an issue: you could write a constant 1e-30 and be confident of that exact value being used in floating-point operations.  However, when the hardware works in a base other than 10, then human-specified numbers have to be converted to that base, and then converted back again at output time.  The two base conversions are rarely exact, and unwanted rounding errors are introduced.

As computers usually represent numbers in a base other than 10, numbers often must be converted to and from different bases, and rounding errors can occur during conversion.  This problem is solved in C using hexademical floating-point constants.  For example, +0x1.fffffcp-1 is the number that is the IEEE 754 32-bit value closest to, but below, 1.0. The significand is represented as a hexadecimal fraction, and the power of two is written in decimal following the exponent letter p (the traditional exponent letter e is not possible, because it is a hexadecimal digit).

In printf and scanf and related functions, you can use the %a and %A format specifiers for writing and reading hexadecimal floating-point values.  %a writes them with lower case letters and %A writes them with upper case letters.  For instance, this code reproduces our sample number:

     printf ("%a\n", 1.0 - pow (2.0, -23));
         -| 0x1.fffffcp-1

The strtod family was similarly extended to recognize numbers in that new format.

If you want to ensure exact data representation for transfer of floating-point numbers between C programs on different computers, then hexadecimal constants are an optimum choice.

Handling Infinity
=================

As we noted earlier, the IEEE 754 model of computing is not to stop the program when exceptional conditions occur.  It takes note of exceptional values or conditions by setting sticky exception flags, or by producing results with the special values Infinity and QNaN.  In this section, we discuss Infinity; see Handling NaN for the other.

In GNU C, you can create a value of negative Infinity in software like this:

     double x;

     x = -1.0 / 0.0;

GNU C supplies the __builtin_inf, __builtin_inff, and __builtin_infl macros, and the GNU C Library provides the INFINITY macro, all of which are compile-time constants for positive infinity.

GNU C also provides a standard function to test for an Infinity: isinf (x) returns 1 if the argument is a signed infinity, and 0 if not.

Infinities can be compared, and all Infinities of the same sign are equal: there is no notion in IEEE 754 arithmetic of different kinds of Infinities, as there are in some areas of mathematics.  Positive Infinity is larger than any finite value, and negative Infinity is smaller than any finite value.

Infinities propagate in addition, subtraction, multiplication, and square root, but in division, they disappear, because of the rule that finite / Infinity is 0.0.  Thus, an overflow in an intermediate computation that produces an Infinity is likely to be noticed later in the final results.  The programmer can then decide whether the overflow is expected, and acceptable, or whether the code possibly has a bug, or needs to be run in higher precision, or redesigned to avoid the production of the Infinity.

Handling NaN
============

NaNs are not numbers: they represent values from computations that produce undefined results.  They have a distinctive property that makes them unlike any other floating-point value: they are unequal to everything, including themselves!  Thus, you can write a test for a NaN like this:

     if (x != x)
       printf ("x is a NaN\n");

This test works in GNU C, but some compilers might evaluate that test expression as false without properly checking for the NaN value. A more portable way to test for NaN is to use the isnan function declared in math.h:

     if (isnan (x))
       printf ("x is a NaN\n");

See Floating Point Classes (The GNU C Library Reference Manual).

One important use of NaNs is marking of missing data.  For example, in statistics, such data must be omitted from computations.  Use of any particular finite value for missing data would eventually collide with real data, whereas such data could never be a NaN, so it is an ideal marker.  Functions that deal with collections of data that may have holes can be written to test for, and ignore, NaN values.

It is easy to generate a NaN in computations: evaluating 0.0 / 0.0 is the commonest way, but Infinity - Infinity, Infinity / Infinity, and sqrt (-1.0) also work. Functions that receive out-of-bounds arguments can choose to return a stored NaN value, such as with the NAN macro defined in math.h, but that does not set the invalid operand exception flag, and that can fool some programs.

Like Infinity, NaNs propagate in computations, but they are even stickier, because they never disappear in division. Thus, once a NaN appears in a chain of numerical operations, it is almost certain to pop out into the final results.  The programmer has to decide whether that is expected, or whether there is a coding or algorithmic error that needs repair.

In general, when function gets a NaN argument, it usually returns a NaN.  However, there are some exceptions in the math-library functions that you need to be aware of, because they violate the NaNs-always-propagate rule:

  * pow (x, 0.0) always returns 1.0, even if x is 0.0, Infinity, or a NaN.

  * pow (1, y) always returns 1, even if y is a NaN.

  * hypot (INFINITY, y) and hypot (-INFINITY, y) both always return INFINITY, even if y is a Nan.

  * If just one of the arguments to fmax (x, y) or fmin (x, y) is a NaN, it returns the other argument.  If both arguments are NaNs, it returns a NaN, but there is no requirement about where it comes from: it could be x, or y, or some other quiet NaN.

NaNs are also used for the return values of math-library functions where the result is not representable in real arithmetic, or is mathematically undefined or uncertain, such as sqrt (-1.0) and sin (Infinity).  However, note that a result that is merely too big to represent should always produce an Infinity, such as with exp (1000.0) (too big) and exp (Infinity) (truly infinite).

Signed Zeros
============

The sign of zero is significant, and important, because it records the creation of a value that is too small to represent, but came from either the negative axis, or from the positive axis.  Such fine distinctions are essential for proper handling of branch cuts in complex arithmetic (see Complex Arithmetic).

The key point about signed zeros is that in comparisons, their sign does not matter: 0.0 == -0.0 must always evaluate to 1 (true).  However, they are not the same number, and -0.0 in C code stands for a negative zero.

Scaling by Powers of the Base
=============================

We have discussed rounding errors several times in this chapter, but it is important to remember that when results require no more bits than the exponent and significand bits can represent, those results are exact.

One particularly useful exact operation is scaling by a power of the base.  While one, in principle, could do that with code like this:

     y = x * pow (2.0, (double)k);   /* Undesirable scaling: avoid! */

that is not advisable, because it relies on the quality of the math-library power function, and that happens to be one of the most difficult functions in the C math library to make accurate. What is likely to happen on many systems is that the returned value from pow will be close to a power of two, but slightly different, so the subsequent multiplication introduces rounding error.

The correct, and fastest, way to do the scaling is either via the traditional C library function, or with its C99 equivalent:

     y = ldexp (x, k);            /* Traditional pre-C99 style. */
     y = scalbn (x, k);           /* C99 style. */

Both functions return x * 2**k. See Normalization Functions (The GNU C Library Reference Manual).

Rounding Control
================

Here we describe how to specify the rounding mode at run time.  System header file fenv.h provides the prototypes for these functions. See Rounding (The GNU C Library Reference Manual).

That header file also provides constant names for the four rounding modes: FE_DOWNWARD, FE_TONEAREST, FE_TOWARDZERO, and FE_UPWARD.

The function fegetround examines and returns the current rounding mode.  On a platform with IEEE 754 floating point, the value will always equal one of those four constants. On other platforms, it may return a negative value. The function fesetround sets the current rounding mode.

Changing the rounding mode can be slow, so it is useful to minimize the number of changes.  For interval arithmetic, we seem to need three changes for each operation, but we really only need two, because we can write code like this example for interval addition of two reals:

     {
       struct interval_double
         {
           double hi, lo;
         } v;
       extern volatile double x, y;
       int rule;

       rule = fegetround ();

       if (fesetround (FE_UPWARD) == 0)
         {
           v.hi = x + y;
           v.lo = -(-x - y);
         }
       else
         fatal ("ERROR: failed to change rounding rule");

       if (fesetround (rule) != 0)
         fatal ("ERROR: failed to restore rounding rule");
     }

The volatile qualifier (see volatile) is essential on x86 platforms to prevent an optimizing compiler from producing the same value for both bounds.

Machine Epsilon
===============

In any floating-point system, three attributes are particularly important to know: base (the number that the exponent specifies a power of), precision (number of digits in the significand), and range (difference between most positive and most negative values).  The allocation of bits between exponent and significand decides the answers to those questions.

A measure of the precision is the answer to the question: what is the smallest number that can be added to 1.0 such that the sum differs from 1.0?  That number is called the machine epsilon.

We could define the needed machine-epsilon constants for float, double, and long double like this:

     static const float  epsf = 0x1p-23;  /* about 1.192e-07 */
     static const double eps  = 0x1p-52;  /* about 2.220e-16 */
     static const long double epsl = 0x1p-63;  /* about 1.084e-19 */

Instead of the hexadecimal constants, we could also have used the Standard C macros, FLT_EPSILON, DBL_EPSILON, and LDBL_EPSILON.

It is useful to be able to compute the machine epsilons at run time, and we can easily generalize the operation by replacing the constant 1.0 with a user-supplied value:

     double
     macheps (double x)
     { /* Return machine epsilon for X,  */
        /* such that X + macheps (X) > X.  */
       static const double base = 2.0;
       double eps;

       if (isnan (x))
           eps = x;
       else
         {
           eps = (x == 0.0) ? 1.0 : x;

           while ((x + eps / base) != x)
               eps /= base;          /* Always exact!  */
         }

       return (eps);
     }

If we call that function with arguments from 0 to 10, as well as Infinity and NaN, and print the returned values in hexadecimal, we get output like this:

     macheps (  0) = 0x1.0000000000000p-1074
     macheps (  1) = 0x1.0000000000000p-52
     macheps (  2) = 0x1.0000000000000p-51
     macheps (  3) = 0x1.8000000000000p-52
     macheps (  4) = 0x1.0000000000000p-50
     macheps (  5) = 0x1.4000000000000p-51
     macheps (  6) = 0x1.8000000000000p-51
     macheps (  7) = 0x1.c000000000000p-51
     macheps (  8) = 0x1.0000000000000p-49
     macheps (  9) = 0x1.2000000000000p-50
     macheps ( 10) = 0x1.4000000000000p-50
     macheps (Inf) = infinity
     macheps (NaN) = nan

Notice that macheps has a special test for a NaN to prevent an infinite loop.

Our code made another test for a zero argument to avoid getting a zero return.  The returned value in that case is the smallest representable floating-point number, here the subnormal value 2**(-1074), which is about 4.941e-324.

No special test is needed for an Infinity, because the eps-reduction loop then terminates at the first iteration.

Our macheps function here assumes binary floating point; some architectures may differ.

The C library includes some related functions that can also be used to determine machine epsilons at run time:

     #include <math.h>           /* Include for these prototypes. */

     double      nextafter  (double x, double y);
     float       nextafterf (float x, float y);
     long double nextafterl (long double x, long double y);

These return the machine number nearest X in the direction of Y.  For example, nextafter (1.0, 2.0) produces the same result as 1.0 + macheps (1.0) and 1.0 + DBL_EPSILON. See FP Bit Twiddling (The GNU C Library Reference Manual).

It is important to know that the machine epsilon is not symmetric about all numbers.  At the boundaries where normalization changes the exponent, the epsilon below X is smaller than that just above X by a factor 1 / base.  For example, macheps (1.0) returns +0x1p-52, whereas macheps (-1.0) returns +0x1p-53.  Some authors distinguish those cases by calling them the positive and negative, or big and small, machine epsilons.  You can produce their values like this:

     eps_neg = 1.0 - nextafter (1.0, -1.0);
     eps_pos = nextafter (1.0, +2.0) - 1.0;

If X is a variable, such that you do not know its value at compile time, then you can substitute literal Y values with either -inf() or +inf(), like this:

     eps_neg = x - nextafter (x, -inf ());
     eps_pos = nextafter (x, +inf() - x);

In such cases, if X is Infinity, then the nextafter functions return Y if X equals Y.  Our two assignments then produce +0x1.fffffffffffffp+1023 (that is a hexadecimal floating point constant and its value is around 1.798e+308; see Floating Constants) for EPS_NEG, and Infinity for EPS_POS.  Thus, the call nextafter (INFINITY, -INFINITY) can be used to find the largest representable finite number, and with the call nextafter (0.0, 1.0), the smallest representable number (here, 0x1p-1074 (about 4.491e-324), a number that we saw before as the output from macheps (0.0)).

Complex Arithmetic
==================

We've already looked at defining and referring to complex numbers (see Complex Data Types).  What is important to discuss here are some issues that are unlikely to be obvious to programmers without extensive experience in both numerical computing, and in complex arithmetic in mathematics.

The first important point is that, unlike real arithmetic, in complex arithmetic, the danger of significance loss is pervasive, and affects every one of the basic operations, and almost all of the math-library functions.  To understand why, recall the rules for complex multiplication and division:

     a = u + I*v              /* First operand. */
     b = x + I*y              /* Second operand. */

     prod = a * b
          = (u + I*v) * (x + I*y)
          = (u * x - v * y) + I*(v * x + u * y)

     quo  = a / b
          = (u + I*v) / (x + I*y)
          = [(u + I*v) * (x - I*y)] / [(x + I*y) * (x - I*y)]
          = [(u * x + v * y) + I*(v * x - u * y)] / (x**2 + y**2)

There are four critical observations about those formulas:

  * the multiplications on the right-hand side introduce the possibility of premature underflow or overflow;

  * the products must be accurate to twice working precision;

  * there is always one subtraction on the right-hand sides that is subject to catastrophic significance loss; and

  * complex multiplication has up to six rounding errors, and complex division has ten rounding errors.

Another point that needs careful study is the fact that many functions in complex arithmetic have branch cuts.  You can view a function with a complex argument, f (z), as f (x + I*y), and thus, it defines a relation between a point (x, y) on the complex plane with an elevation value on a surface.  A branch cut looks like a tear in that surface, so approaching the cut from one side produces a particular value, and from the other side, a quite different value.  Great care is needed to handle branch cuts properly, and even small numerical errors can push a result from one side to the other, radically changing the returned value.  As we reported earlier, correct handling of the sign of zero is critically important for computing near branch cuts.

The best advice that we can give to programmers who need complex arithmetic is to always use the highest precision available, and then to carefully check the results of test calculations to gauge the likely accuracy of the computed results.  It is easy to supply test values of real and imaginary parts where all five basic operations in complex arithmetic, and almost all of the complex math functions, lose all significance, and fail to produce even a single correct digit.

Even though complex arithmetic makes some programming tasks easier, it may be numerically preferable to rework the algorithm so that it can be carried out in real arithmetic.  That is commonly possible in matrix algebra.

GNU C can perform code optimization on complex number multiplication and division if certain boundary checks will not be needed.  The command-line option -fcx-limited-range tells the compiler that a range reduction step is not needed when performing complex division, and that there is no need to check if a complex multiplication or division results in the value Nan + I*NaN.  By default these checks are enabled.  You can explicitly enable them with the -fno-cx-limited-range option.

Round-Trip Base Conversion
==========================

Most numeric programs involve converting between base-2 floating-point numbers, as represented by the computer, and base-10 floating-point numbers, as entered and handled by the programmer.  What might not be obvious is the number of base-2 bits vs. base-10 digits required for each representation.  Consider the following tables showing the number of decimal digits representable in a given number of bits, and vice versa:

binary in     24   53   64   113   237
decimal out   9    17   21   36    73

decimal in    7    16   34    70
binary out    25   55   114   234

We can compute the table numbers with these two functions:

     int
     matula(int nbits)
     {   /* Return output decimal digits needed for nbits-bits input. */
         return ((int)ceil((double)nbits / log2(10.0) + 1.0));
     }

     int
     goldberg(int ndec)
     {   /* Return output bits needed for ndec-digits input. */
         return ((int)ceil((double)ndec / log10(2.0) + 1.0));
     }

One significant observation from those numbers is that we cannot achieve correct round-trip conversion between the decimal and binary formats in the same storage size!  For example, we need 25 bits to represent a 7-digit value from the 32-bit decimal format, but the binary format only has 24 available.  Similar observations hold for each of the other conversion pairs.

The general input/output base-conversion problem is astonishingly complicated, and solutions were not generally known until the publication of two papers in 1990 that are listed later near the end of this chapter.  For the 128-bit formats, the worst case needs more than 11,500 decimal digits of precision to guarantee correct rounding in a binary-to-decimal conversion!

For further details see the references for Bennett Goldberg and David Matula.

Further Reading
===============

The subject of floating-point arithmetic is much more complex than many programmers seem to think, and few books on programming languages spend much time in that area.  In this chapter, we have tried to expose the reader to some of the key ideas, and to warn of easily overlooked pitfalls that can soon lead to nonsensical results.  There are a few good references that we recommend for further reading, and for finding other important material about computer arithmetic:

We include URLs for these references when we were given them, when they are morally legitimate to recommend; we have omitted the URLs that are paywalled or that require running nonfree JavaScript code in order to function.  We hope you can find morally legitimate sites where you can access these works.

  * sort-key: Abbott Paul H. Abbott and 15 others, Architecture and software support in IBM S/390 Parallel Enterprise Servers for IEEE Floating-Point arithmetic, IBM Journal of Research and Development 43(5/6) 723--760 (1999), This article gives a good description of IBM's algorithm for exact decimal-to-binary conversion, complementing earlier ones by Clinger and others.

  * sort-key: Beebe Nelson H. F. Beebe, The Mathematical-Function Computation Handbook: Programming Using the MathCW Portable Software Library, Springer (2017). This book describes portable implementations of a large superset of the mathematical functions available in many programming languages, extended to a future 256-bit format (70 decimal digits), for both binary and decimal floating point.  It includes a substantial portion of the functions described in the famous NIST Handbook of Mathematical Functions, Cambridge (2018), ISBN 0-521-19225-0. See https://www.math.utah.edu/pub/mathcw/ for compilers and libraries.

  * sort-key: Clinger-1990 William D. Clinger, How to Read Floating Point Numbers Accurately,  ACM SIGPLAN Notices 25(6) 92--101 (June 1990), https://doi.org/10.1145/93548.93557. See also the papers by Steele & White.

  * sort-key: Goldberg-1967 I. Bennett Goldberg, 27  Bits Are Not Enough For 8-Digit Accuracy, Communications of the ACM 10(2) 105--106 (February 1967), https://doi.acm.org/10.1145/363067.363112.  This paper, and its companions by David Matula, address the base-conversion problem, and show that the naive formulas are wrong by one or two digits.

  * sort-key: Goldberg-1991 David Goldberg, What Every Computer Scientist Should Know About Floating-Point Arithmetic, ACM Computing Surveys 23(1) 5--58 (March 1991), corrigendum 23(3) 413 (September 1991), https://doi.org/10.1145/103162.103163. This paper has been widely distributed, and reissued in vendor programming-language documentation.  It is well worth reading, and then rereading from time to time.

  * sort-key: Juffa Norbert Juffa and Nelson H. F. Beebe, A Bibliography of Publications on Floating-Point Arithmetic, https://www.math.utah.edu/pub/tex/bib/fparith.bib. This is the largest known bibliography of publications about floating-point, and also integer, arithmetic.  It is actively maintained, and in mid 2019, contains more than 6400 references to original research papers, reports, theses, books, and Web sites on the subject matter.  It can be used to locate the latest research in the field, and the historical coverage dates back to a 1726 paper on signed-digit arithmetic, and an 1837 paper by Charles Babbage, the intellectual father of mechanical computers.  The entries for the Abbott, Clinger, and Steele & White papers cited earlier contain pointers to several other important related papers on the base-conversion problem.

  * sort-key: Kahan William Kahan, Branch Cuts for Complex Elementary Functions, or Much Ado About Nothing's Sign Bit, (1987), https://people.freebsd.org/~das/kahan86branch.pdf. This Web document about the fine points of complex arithmetic also appears in the volume edited by A. Iserles and M. J. D. Powell, The State of the Art in Numerical Analysis: Proceedings of the Joint IMA/SIAM Conference on the State of the Art in Numerical Analysis held at the University of Birmingham, 14--18 April 1986, Oxford University Press (1987), ISBN 0-19-853614-3 (xiv + 719 pages).  Its author is the famous chief architect of the IEEE 754 arithmetic system, and one of the world's greatest experts in the field of floating-point arithmetic.  An entire generation of his students at the University of California, Berkeley, have gone on to careers in academic and industry, spreading the knowledge of how to do floating-point arithmetic right.

  * sort-key: Knuth Donald E. Knuth, A Simple Program Whose Proof Isn't, in Beauty is our business: a birthday salute to Edsger W. Dijkstra, W. H. J. Feijen, A. J. M. van Gasteren, D. Gries, and J. Misra (eds.), Springer (1990), ISBN 1-4612-8792-8,  This book chapter supplies a correctness proof of the decimal to binary, and binary to decimal, conversions in fixed-point arithmetic in the TeX typesetting system.  The proof evaded its author for a dozen years.

  * sort-key: Matula-1968a David W. Matula, In-and-out conversions, Communications of the ACM 11(1) 57--50 (January 1968), https://doi.org/10.1145/362851.362887.

  * sort-key: Matula-1968b David W. Matula, The Base Conversion Theorem, Proceedings of the American Mathematical Society 19(3) 716--723 (June 1968).  See also other papers here by this author, and by I. Bennett Goldberg.

  * sort-key: Matula-1970 David W. Matula, A Formalization of Floating-Point Numeric Base Conversion, IEEE Transactions on Computers C-19(8) 681--692 (August 1970),

  * sort-key: Muller-2010 Jean-Michel Muller and eight others, Handbook of Floating-Point Arithmetic, Birkhäuser-Boston (2010), ISBN 0-8176-4704-X (xxiii + 572 pages).  This is a comprehensive treatise from a French team who are among the world's greatest experts in floating-point arithmetic, and among the most prolific writers of research papers in that field.  They have much to teach, and their book deserves a place on the shelves of every serious numerical programmer.

  * sort-key: Muller-2018 Jean-Michel Muller and eight others, Handbook of Floating-Point Arithmetic, Second edition, Birkhäuser-Boston (2018), ISBN 3-319-76525-6 (xxv + 627 pages).  This is a new edition of the preceding entry.

  * sort-key: Overton Michael Overton, Numerical Computing with IEEE Floating Point Arithmetic, Including One Theorem, One Rule of Thumb, and One Hundred and One Exercises, SIAM (2001), ISBN 0-89871-482-6 (xiv + 104 pages), This is a small volume that can be covered in a few hours.

  * sort-key: Steele-1990 Guy L. Steele Jr. and Jon L. White, How to Print Floating-Point Numbers Accurately, ACM SIGPLAN Notices 25(6) 112--126 (June 1990), https://doi.org/10.1145/93548.93559. See also the papers by Clinger.

  * sort-key: Steele-2004 Guy L. Steele Jr. and Jon L. White, Retrospective: How to Print Floating-Point Numbers Accurately, ACM SIGPLAN Notices 39(4) 372--389 (April 2004),  Reprint of 1990 paper, with additional commentary.

  * sort-key: Sterbenz Pat H. Sterbenz, Floating Point Computation, Prentice-Hall (1974), ISBN 0-13-322495-3 (xiv + 316 pages).  This often-cited book provides solid coverage of what floating-point arithmetic was like before the introduction of IEEE 754 arithmetic.