Copyright @copyright{} 2017,2019,2020 Free Software Foundation, Inc. @c Copyright @copyright{} 2020 Richard Stallman and Free Software Foundation, Inc. This is part of the GNU C Intro and Reference Manual and covered by its license. @node Floating Point in Depth @chapter Floating Point in Depth @menu * Floating Representations:: * Floating Type Specs:: * Special Float Values:: * Invalid Optimizations:: * Exception Flags:: * Exact Floating-Point:: * Rounding:: * Rounding Issues:: * Significance Loss:: * Fused Multiply-Add:: * Error Recovery:: @c * Double-Rounding Problems:: * Exact Floating Constants:: * Handling Infinity:: * Handling NaN:: * Signed Zeros:: * Scaling by the Base:: * Rounding Control:: * Machine Epsilon:: * Complex Arithmetic:: * Round-Trip Base Conversion:: * Further Reading:: @end menu @node Floating Representations @section Floating-Point Representations @cindex floating-point representations @cindex representation of floating-point numbers @cindex IEEE 754-2008 Standard Storing numbers as @dfn{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 @dfn{significand} (also called the @dfn{mantissa}), and a power of a fixed base. GNU C uses the floating-point representations specified by the @cite{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 @code{float}, @code{double}, and @code{long double}. GNU C supports the 16-bit floating point type @code{_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 @dfn{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 @minus{}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. @node Floating Type Specs @section Floating-Point Type Specifications The standard library header file @file{float.h} defines a number of constants that describe the platform's implementation of floating-point types @code{float}, @code{double} and @code{long double}. They include: @findex FLT_MIN @findex DBL_MIN @findex LDBL_MIN @findex FLT_HAS_SUBNORM @findex DBL_HAS_SUBNORM @findex LDBL_HAS_SUBNORM @findex FLT_TRUE_MIN @findex DBL_TRUE_MIN @findex LDBL_TRUE_MIN @findex FLT_MAX @findex DBL_MAX @findex LDBL_MAX @findex FLT_DECIMAL_DIG @findex DBL_DECIMAL_DIG @findex LDBL_DECIMAL_DIG @table @code @item FLT_MIN @itemx DBL_MIN @itemx LDBL_MIN Defines the minimum normalized positive floating-point values that can be represented with the type. @item FLT_HAS_SUBNORM @itemx DBL_HAS_SUBNORM @itemx LDBL_HAS_SUBNORM Defines if the floating-point type supports subnormal (or ``denormalized'') numbers or not (@pxref{subnormal numbers}). @item FLT_TRUE_MIN @itemx DBL_TRUE_MIN @itemx LDBL_TRUE_MIN Defines the minimum positive values (including subnormal values) that can be represented with the type. @item FLT_MAX @itemx DBL_MAX @itemx LDBL_MAX Defines the largest values that can be represented with the type. @item FLT_DECIMAL_DIG @itemx DBL_DECIMAL_DIG @itemx LDBL_DECIMAL_DIG Defines the number of decimal digits @code{n} such that any floating-point number that can be represented in the type can be rounded to a floating-point number with @code{n} decimal digits, and back again, without losing any precision of the value. @end table @node Special Float Values @section Special Floating-Point Values @cindex special floating-point values @cindex floating-point values, special IEEE floating point provides for special values that are not ordinary numbers. @table @asis @item infinities @code{+Infinity} and @code{-Infinity} are two different infinite values, one positive and one negative. These result from operations such as @code{1 / 0}, @code{Infinity + Infinity}, @code{Infinity * Infinity}, and @code{Infinity + @var{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. @xref{Handling Infinity}, for more about working with infinities. @item NaNs (not a number) @cindex QNaN @cindex SNaN 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 @code{0 / 0}, @code{sqrt (-1)}, @code{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 @emph{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 @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. @xref{Handling NaN}, for more about working with NaNs. @item subnormal numbers @cindex subnormal numbers @cindex underflow, floating @cindex floating underflow @anchor{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 @dfn{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 @dfn{normalized} representable value (i.e., one in which the leading significand bit is @code{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 @dfn{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 @dfn{subnormal}. You can use the library functions @code{fesetround} and @code{fegetround} to set and get the rounding mode. Rounding modes are defined (if supported by the platform) in @code{fenv.h} as: @code{FE_UPWARD} to round toward positive infinity; @code{FE_DOWNWARD} to round toward negative infinity; @code{FE_TOWARDZERO} to round toward zero; and @code{FE_TONEAREST} to round to the nearest representable value, the default mode. It is best to use @code{FE_TONEAREST} except when there is a special need for some other mode. @end table @node Invalid Optimizations @section Invalid Optimizations @cindex invalid optimizations in floating-point arithmetic @cindex floating-point arithmetic invalid optimizations Signed zeros, Infinity, and NaN invalidate some optimizations by programmers and compilers that might otherwise have seemed obvious: @itemize @bullet @item @code{x + 0} and @code{x - 0} are not the same as @code{x} when @code{x} is zero, because the result depends on the rounding rule. @xref{Rounding}, for more about rounding rules. @item @code{x * 0.0} is not the same as @code{0.0} when @code{x} is Infinity, a NaN, or negative zero. @item @code{x / x} is not the same as @code{1.0} when @code{x} is Infinity, a NaN, or zero. @item @code{(x - y)} is not the same as @code{-(y - x)} because when the operands are finite and equal, one evaluates to @code{+0} and the other to @code{-0}. @item @code{x - x} is not the same as @code{0.0} when @var{x} is Infinity or a NaN. @item @code{x == x} and @code{x != x} are not equivalent to @code{1} and @code{0} when @var{x} is a NaN. @item @code{x < y} and @code{isless (x, y)} are not equivalent, because the first sets a sticky exception flag (@pxref{Exception Flags}) when an operand is a NaN, whereas the second does not affect that flag. The same holds for the other @code{isxxx} functions that are companions to relational operators. @xref{FP Comparison Functions, , , libc, The GNU C Library Reference Manual}. @end itemize The @option{-funsafe-math-optimizations} option enables these optimizations. @node Exception Flags @section Floating Arithmetic Exception Flags @cindex floating arithmetic exception flags @cindex exception flags (floating point) @cindex sticky exception flags (floating point) @cindex floating overflow @cindex overflow, floating @cindex floating underflow @cindex underflow, floating @dfn{Sticky exception flags} record the occurrence of particular conditions: once set, they remain set until the program explicitly clears them. The conditions include @emph{invalid operand}, @emph{division-by_zero}, @emph{inexact result} (i.e., one that required rounding), @emph{underflow}, and @emph{overflow}. Some extended floating-point designs offer several additional exception flags. The functions @code{feclearexcept}, @code{feraiseexcept}, @code{fetestexcept}, @code{fegetexceptflags}, and @code{fesetexceptflags} provide a standardized interface to those flags. @xref{Status bit operations, , , libc, The GNU C Library Reference Manual}. One important use of those @anchor{fetestexcept}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 @emph{inexact result} flag with a call to @code{feclearexcept (FE_INEXACT)}, do the computation, and then test the flag with @code{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). @c ===================================================================== @ignore @node IEEE 754 Decimal Arithmetic @section IEEE 754 Decimal Arithmetic @cindex IEEE 754 decimal arithmetic One of the difficulties that users of computers for numerical work face, whether they realize it or not, is that the computer does not operate in the number base that most people are familiar with. As a result, input decimal fractions must be converted to binary floating-point values for use in computations, and then the final results converted back to decimal for humans. Because the precision is finite and limited, and because algorithms for correct round-trip conversion between number bases were not known until the 1990s, and are still not implemented on most systems and most programming languages, the result is frequent confusion for users, when a simple expression like @code{3.0*(1.0/3.0)} evaluates to 0.999999 instead of the expected 1.0. Here is an example from a floating-point calculator that allows rounding-mode control, with the mode set to @emph{round-towards-zero}: @example for (k = 1; k <= 10; ++k) (void)printf ("%2d\t%10.6f\n", k, k*(1.0/k)) 1 1.000000 2 1.000000 3 0.999999 4 1.000000 5 0.999999 6 0.999999 7 0.999999 8 1.000000 9 0.999999 10 0.999999 @end example Increasing working precision can sometimes help by reducing intermediate rounding errors, but the reality is that when fractional values are involved, @emph{no amount} of extra precision can suffice for some computations. For example, the nice decimal value @code{1/10} in C99-style binary representation is @code{+0x1.999999999999ap-4}; that final digit @code{a} is the rounding of an infinite string of @code{9}'s. Financial computations in particular depend critically on correct arithmetic, and the losses due to rounding errors can be large, especially for businesses with large numbers of small transactions, such as grocery stores and telephone companies. Tax authorities are particularly picky, and demand specific rounding rules, including one that instead of rounding ties to the nearest number, rounds instead in the direction that favors the taxman. Programming languages used for business applications, notably the venerable Cobol language, have therefore always implemented financial computations in @emph{fixed-point decimal arithmetic} in software, and because of the large monetary amounts that must be processed, successive Cobol standards have increased the minimum number size from 18 to 32 decimal digits, and the most recent one requires a decimal exponent range of at least @code{[-999, +999]}. The revised IEEE 754-2008 standard therefore requires decimal floating-point arithmetic, as well as the now-widely used binary formats from 1985. Like the binary formats, the decimal formats also support Infinity, NaN, and signed zero, and the five basic operations are also required to produce correctly rounded representations of infinitely precise exact results. However, the financial applications of decimal arithmetic introduce some new features: @itemize @bullet @item There are three decimal formats occupying 32, 64, and 128 bits of storage, and offering exactly 7, 16, and 34 decimal digits of precision. If one size has @code{n} digits, the next larger size has @code{2 n + 2} digits. Thus, a future 256-bit format would supply 70 decimal digits, and at least one library already supports the 256-bit binary and decimal formats. @item Decimal arithmetic has an additional rounding mode, called @emph{round-ties-to-away-from-zero}, meaning that a four-digit rounding of @code{1.2345} is @code{1.235}, and @code{-1.2345} becomes @code{-1.235}. That rounding mode is mandated by financial laws in several countries. @item The decimal significand is an @anchor{decimal-significand}@emph{integer}, instead of a fractional value, and trailing zeros are only removed at user request. That feature allows floating-point arithmetic to emulate the @emph{fixed-point arithmetic} traditionally used in financial computations. @end itemize @noindent We can easily estimate how many digits are likely to be needed for financial work: seven billion people on Earth, with an average annual income of less than US$10,000, means a world financial base that can be represented in just 15 decimal digits. Even allowing for alternate currencies, future growth, multiyear accounting, and intermediate computations, the 34 digits supplied by the 128-bit format are more than enough for financial purposes. We return to decimal arithmetic later in this chapter (@pxref{More on decimal floating-point arithmetic}), after we have covered more about floating-point arithmetic in general. @c ===================================================================== @end ignore @node Exact Floating-Point @section Exact Floating-Point Arithmetic @cindex exact floating-point arithmetic @cindex floating-point arithmetic, exact 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 @emph{exact}. It is easy to predict how many digits are needed for the results of arithmetic operations: @itemize @bullet @item addition and subtraction of two @var{n}-digit values with the @emph{same} exponent require at most @code{@var{n} + 1} digits, but when the exponents differ, many more digits may be needed; @item multiplication of two @var{n}-digit values requires exactly 2 @var{n} digits; @item although integer division produces a quotient and a remainder of no more than @var{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. @end itemize Whenever a result requires more than @var{n} digits, rounding is needed. @c ===================================================================== @node Rounding @section Rounding @cindex 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 @dfn{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 @code{FE_INEXACT} exception flag (@pxref{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 @dfn{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: @itemize * @item @emph{round-to-nearest}, with ties rounded to an even integer; @item @emph{round-up}, towards @code{+Infinity}; @item @emph{round-down}, towards @code{-Infinity}; @item @emph{round-towards-zero}. @end itemize Under those four rounding rules, a decimal value @code{-1.2345} that is to be rounded to a four-digit result would become @code{-1.234}, @code{-1.234}, @code{-1.235}, and @code{-1.234}, respectively. The default rounding rule is @emph{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 @emph{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, @emph{round-up} and @emph{round-down}, are essential for implementing @dfn{interval arithmetic}, whereby each arithmetic operation produces lower and upper bounds that are guaranteed to enclose the exact result. @xref{Rounding Control}, for details on getting and setting the current rounding mode. @node Rounding Issues @section Rounding Issues @cindex rounding issues (floating point) @cindex floating-point 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. @example #include #include int main (void) @{ double x = 100.0; double y; 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; @} @end example @noindent Here is the output: @example 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 @end example After 50 iterations, @code{y} has barely one correct digit, and soon after, there are no correct digits. @c ===================================================================== @node Significance Loss @section Significance Loss @cindex significance loss (floating point) @cindex floating-point significance loss A much more serious source of error in floating-point computation is @dfn{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 @emph{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 @code{@var{x}**(2*@var{n} + 1) / (2*@var{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 @var{k} terms of the series for @code{sin (@var{x})}, and compare the computed sums with known-to-be-accurate library functions: @example 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 @end example @noindent We have added underscores to the numbers to enhance readability. The @code{sinh (@var{x})} series with positive terms can be summed to high accuracy. By contrast, the series for @code{sin (@var{x})} suffers increasing significance loss, so that when @var{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 @code{long double} are important. @c Formerly mentioned @code{__float128} @c ===================================================================== @node Fused Multiply-Add @section Fused Multiply-Add @cindex fused multiply-add in floating-point computations @cindex floating-point fused multiply-add In 1990, when IBM introduced the POWER architecture, the CPU provided a previously unknown instruction, the @dfn{fused multiply-add} (FMA). It computes the value @code{x * y + z} with an @strong{exact} double-length product, followed by an addition with a @emph{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 @code{0.0} and @code{1.0}, and the normal @emph{multiply} and @emph{add} instructions are just wrappers around the FMA that compute @code{x * y + 0.0} and @code{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 @code{x * y - z} (FMS), @code{-x * y + z} (FNMA), and @code{-x * y - z} (FNMS). @c The IEEE 754-2008 revision requires implementations to provide @c the FMA, as a sixth basic operation. The functions @code{fmaf}, @code{fma}, and @code{fmal} implement fused multiply-add for the @code{float}, @code{double}, and @code{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 @option{-ffp-contract=fast} option to allow generation of FMA instructions, or @option{-ffp-contract=off} to disallow it. @c ===================================================================== @node Error Recovery @section Error Recovery @cindex error recovery (floating point) @cindex floating-point 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 @emph{round-to-nearest} rounding mode, we can determine the error in a sum like this: @example 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; @} @end example @noindent @cindex twosum Now, @code{x + y} is @emph{exactly} represented by @code{sum + err}. This basic operation, which has come to be called @dfn{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 @code{+} and @code{-} operators. We used the @code{volatile} qualifier (@pxref{volatile}) in the declaration of the variables, which forces the compiler to store and retrieve them from memory, and prevents the compiler from optimizing @code{err = y - ((x + y) - x)} into @code{err = 0}. For multiplication, we can compute the rounding error without magnitude tests with the FMA operation (@pxref{Fused Multiply-Add}), like this: @example volatile double err, prod, x, y; prod = x * y; /* @r{rounded product} */ err = fma (x, y, -prod); /* @r{exact product @code{= @var{prod} + @var{err}}} */ @end example 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 @code{x / y} with @code{quo + err} like this: @example volatile double err, quo, x, y; quo = x / y; err = fma (-quo, y, x) / y; @end example For square root, we can approximate @code{sqrt (x)} with @code{root + err} like this: @example volatile double err, root, x; root = sqrt (x); err = fma (-root, root, x) / (root + root); @end example 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. @c ===================================================================== @ignore @node Double-Rounding Problems @section Double-Rounding Problems @cindex double-rounding problems (floating point) @cindex floating-point double-rounding problems Most developers today use 64-bit x86 processors with a 64-bit operating system, with a Streaming SIMD Extensions (SSE) instruction set. In the past, using a 32-bit x87 instruction set was common, leading to issues described in this section. To offer a few more digits of precision and a wider exponent range, the IEEE 754 Standard included an optional @emph{temporary real} format, with 11 more bits in the significand, and 4 more bits in the biased exponent. Compilers are free to exploit the longer format, and most do so. That is usually a @emph{good thing}, such as in computing a lengthy sum or product, or in implementing the computation of the hypotenuse of a right-triangle as @code{sqrt (x*x + y*y)}: the wider exponent range is critical there for avoiding disastrous overflow or underflow. @findex fesetprec @findex fegetprec However, sometimes it is critical to know what the intermediate precision and rounding mode are, such as in tracking errors with the techniques of the preceding section. Some compilers provide options to prevent the use of the 80-bit format in computations with 64-bit @code{double}, but ensuring that code is always compiled that way may be difficult. The x86 architecture has the ability to force rounding of all operations in the 80-bit registers to the 64-bit storage format, and some systems provide a software interface with the functions @code{fesetprec} and @code{fegetprec}. Unfortunately, neither of those functions is defined by the ISO Standards for C and C++, and consequently, they are not standardly available on all platforms that use the x86 floating-point design. When @code{double} computations are done in the 80-bit format, results necessarily involve a @dfn{double rounding}: first to the 64-bit significand in intermediate operations in registers, and then to the 53-bit significand when the register contents are stored to memory. Here is an example in decimal arithmetic where such a double rounding results in the wrong answer: round @code{1_234_999} from seven to five to four digits. The result is @code{1_235_000}, whereas the correct representation to four significant digits is @code{1_234_000}. @cindex -ffloat-store One way to reduce the use of the 80-bit format is to declare variables as @code{volatile double}: that way, the compiler is required to store and load intermediates from memory, rather than keeping them in 80-bit registers over long sequences of floating-point instructions. Doing so does not, however, eliminate double rounding. The now-common x86-64 architecture has separate sets of 32-bit and 64-bit floating-point registers. The option @option{-float-store} says that floating-point computation should use only those registers, thus eliminating the possibility of double rounding. @end ignore @c ===================================================================== @node Exact Floating Constants @section Exact Floating-Point Constants @cindex exact specification of floating-point constants @cindex floating-point constants, exact specification of 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 @code{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. @cindex hexademical floating-point constants 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, @code{+0x1.fffffcp-1} is the number that is the IEEE 754 32-bit value closest to, but below, @code{1.0}. The significand is represented as a hexadecimal fraction, and the @emph{power of two} is written in decimal following the exponent letter @code{p} (the traditional exponent letter @code{e} is not possible, because it is a hexadecimal digit). In @code{printf} and @code{scanf} and related functions, you can use the @samp{%a} and @samp{%A} format specifiers for writing and reading hexadecimal floating-point values. @samp{%a} writes them with lower case letters and @samp{%A} writes them with upper case letters. For instance, this code reproduces our sample number: @example printf ("%a\n", 1.0 - pow (2.0, -23)); @print{} 0x1.fffffcp-1 @end example @noindent The @code{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. @c ===================================================================== @node Handling Infinity @section Handling Infinity @cindex infinity in floating-point arithmetic @cindex floating-point 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 @dfn{exception flags}, or by producing results with the special values Infinity and QNaN. In this section, we discuss Infinity; @pxref{Handling NaN} for the other. In GNU C, you can create a value of negative Infinity in software like this: @verbatim double x; x = -1.0 / 0.0; @end verbatim GNU C supplies the @code{__builtin_inf}, @code{__builtin_inff}, and @code{__builtin_infl} macros, and the GNU C Library provides the @code{INFINITY} macro, all of which are compile-time constants for positive infinity. GNU C also provides a standard function to test for an Infinity: @code{isinf (x)} returns @code{1} if the argument is a signed infinity, and @code{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 finite value. Infinities propagate in addition, subtraction, multiplication, and square root, but in division, they disappear, because of the rule that @code{finite / Infinity} is @code{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. @c ===================================================================== @node Handling NaN @section Handling NaN @cindex NaN in floating-point arithmetic @cindex not a number @cindex floating-point 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 @emph{unequal to everything, including themselves}! Thus, you can write a test for a NaN like this: @example if (x != x) printf ("x is a NaN\n"); @end example @noindent 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 @code{isnan} function declared in @code{math.h}: @example if (isnan (x)) printf ("x is a NaN\n"); @end example @noindent @xref{Floating Point Classes, , , libc, 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 @code{0.0 / 0.0} is the commonest way, but @code{Infinity - Infinity}, @code{Infinity / Infinity}, and @code{sqrt (-1.0)} also work. Functions that receive out-of-bounds arguments can choose to return a stored NaN value, such as with the @code{NAN} macro defined in @code{math.h}, but that does not set the @emph{invalid operand} exception flag, and that can fool some programs. @cindex NaNs-always-propagate rule 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 @emph{NaNs-always-propagate} rule: @itemize @bullet @item @code{pow (x, 0.0)} always returns @code{1.0}, even if @code{x} is 0.0, Infinity, or a NaN. @item @code{pow (1, y)} always returns @code{1}, even if @code{y} is a NaN. @item @code{hypot (INFINITY, y)} and @code{hypot (-INFINITY, y)} both always return @code{INFINITY}, even if @code{y} is a Nan. @item If just one of the arguments to @code{fmax (x, y)} or @code{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 @code{x}, or @code{y}, or some other quiet NaN. @end itemize 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 @code{sqrt (-1.0)} and @code{sin (Infinity)}. However, note that a result that is merely too big to represent should always produce an Infinity, such as with @code{exp (1000.0)} (too big) and @code{exp (Infinity)} (truly infinite). @c ===================================================================== @node Signed Zeros @section Signed Zeros @cindex signed zeros in floating-point arithmetic @cindex floating-point 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 @dfn{branch cuts} in complex arithmetic (@pxref{Complex Arithmetic}). The key point about signed zeros is that in comparisons, their sign does not matter: @code{0.0 == -0.0} must @emph{always} evaluate to @code{1} (true). However, they are not @emph{the same number}, and @code{-0.0} in C code stands for a negative zero. @c ===================================================================== @node Scaling by the Base @section Scaling by Powers of the Base @cindex scaling floating point by powers of the base @cindex floating-point 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 @emph{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: @example y = x * pow (2.0, (double)k); /* @r{Undesirable scaling: avoid!} */ @end example @noindent 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 @code{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: @example y = ldexp (x, k); /* @r{Traditional pre-C99 style.} */ y = scalbn (x, k); /* @r{C99 style.} */ @end example @noindent Both functions return @code{x * 2**k}. @xref{Normalization Functions, , , libc, The GNU C Library Reference Manual}. @c ===================================================================== @node Rounding Control @section Rounding Control @cindex rounding control (floating point) @cindex floating-point rounding control Here we describe how to specify the rounding mode at run time. System header file @file{fenv.h} provides the prototypes for these functions. @xref{Rounding, , , libc, The GNU C Library Reference Manual}. @noindent That header file also provides constant names for the four rounding modes: @code{FE_DOWNWARD}, @code{FE_TONEAREST}, @code{FE_TOWARDZERO}, and @code{FE_UPWARD}. The function @code{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 @code{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: @example @{ struct interval_double @{ double hi, lo; @} v; 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"); @} @end example @noindent The @code{volatile} qualifier (@pxref{volatile}) is essential on x86 platforms to prevent an optimizing compiler from producing the same value for both bounds. @ignore We no longer discuss the double rounding issue. The code also needs to be compiled with the option @option{-ffloat-store} that prevents use of higher precision for the basic operations, because that would introduce double rounding that could spoil the bounds guarantee of interval arithmetic. @end ignore @c ===================================================================== @node Machine Epsilon @section Machine Epsilon @cindex machine epsilon (floating point) @cindex floating-point machine epsilon In any floating-point system, three attributes are particularly important to know: @dfn{base} (the number that the exponent specifies a power of), @dfn{precision} (number of digits in the significand), and @dfn{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 @code{1.0} such that the sum differs from @code{1.0}? That number is called the @dfn{machine epsilon}. We could define the needed machine-epsilon constants for @code{float}, @code{double}, and @code{long double} like this: @example static const float epsf = 0x1p-23; /* @r{about 1.192e-07} */ static const double eps = 0x1p-52; /* @r{about 2.220e-16} */ static const long double epsl = 0x1p-63; /* @r{about 1.084e-19} */ @end example @noindent Instead of the hexadecimal constants, we could also have used the Standard C macros, @code{FLT_EPSILON}, @code{DBL_EPSILON}, and @code{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 @code{1.0} with a user-supplied value: @example double macheps (double x) @{ /* @r{Return machine epsilon for @var{x},} */ @r{such that @var{x} + macheps (@var{x}) > @var{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; /* @r{Always exact!} */ @} return (eps); @} @end example @noindent If we call that function with arguments from @code{0} to @code{10}, as well as Infinity and NaN, and print the returned values in hexadecimal, we get output like this: @example 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 @end example @noindent Notice that @code{macheps} has a special test for a NaN to prevent an infinite loop. @ignore We no longer discuss double rounding. To ensure that no expressions are evaluated with an intermediate higher precision, we can compile with the @option{-fexcess-precision=standard} option, which tells the compiler that all calculation results, including intermediate results, are to be put on the stack, forcing rounding. @end ignore 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 @code{2**(-1074)}, which is about @code{4.941e-324}. No special test is needed for an Infinity, because the @code{eps}-reduction loop then terminates at the first iteration. Our @code{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: @example #include /* @r{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); @end example @noindent These return the machine number nearest @var{x} in the direction of @var{y}. For example, @code{nextafter (1.0, 2.0)} produces the same result as @code{1.0 + macheps (1.0)} and @code{1.0 + DBL_EPSILON}. @xref{FP Bit Twiddling, , , libc, 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 @var{x} is smaller than that just above @var{x} by a factor @code{1 / base}. For example, @code{macheps (1.0)} returns @code{+0x1p-52}, whereas @code{macheps (-1.0)} returns @code{+0x1p-53}. Some authors distinguish those cases by calling them the @emph{positive} and @emph{negative}, or @emph{big} and @emph{small}, machine epsilons. You can produce their values like this: @example eps_neg = 1.0 - nextafter (1.0, -1.0); eps_pos = nextafter (1.0, +2.0) - 1.0; @end example If @var{x} is a variable, such that you do not know its value at compile time, then you can substitute literal @var{y} values with either @code{-inf()} or @code{+inf()}, like this: @example eps_neg = x - nextafter (x, -inf ()); eps_pos = nextafter (x, +inf() - x); @end example @noindent In such cases, if @var{x} is Infinity, then @emph{the @code{nextafter} functions return @var{y} if @var{x} equals @var{y}}. Our two assignments then produce @code{+0x1.fffffffffffffp+1023} (about 1.798e+308) for @var{eps_neg} and Infinity for @var{eps_pos}. Thus, the call @code{nextafter (INFINITY, -INFINITY)} can be used to find the largest representable finite number, and with the call @code{nextafter (0.0, 1.0)}, the smallest representable number (here, @code{0x1p-1074} (about 4.491e-324), a number that we saw before as the output from @code{macheps (0.0)}). @c ===================================================================== @node Complex Arithmetic @section Complex Arithmetic @cindex complex arithmetic in floating-point calculations @cindex floating-point arithmetic with complex numbers We've already looked at defining and referring to complex numbers (@pxref{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 @emph{pervasive}, and affects @emph{every one} of the basic operations, and @emph{almost all} of the math-library functions. To understand why, recall the rules for complex multiplication and division: @example a = u + I*v /* @r{First operand.} */ b = x + I*y /* @r{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) @end example @noindent There are four critical observations about those formulas: @itemize @bullet @item the multiplications on the right-hand side introduce the possibility of premature underflow or overflow; @item the products must be accurate to twice working precision; @item there is @emph{always} one subtraction on the right-hand sides that is subject to catastrophic significance loss; and @item complex multiplication has up to @emph{six} rounding errors, and complex division has @emph{ten} rounding errors. @end itemize @cindex branch cuts Another point that needs careful study is the fact that many functions in complex arithmetic have @dfn{branch cuts}. You can view a function with a complex argument, @code{f (z)}, as @code{f (x + I*y)}, and thus, it defines a relation between a point @code{(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 @emph{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 @emph{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 @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 @code{Nan + I*NaN}. By default these checks are enabled. You can explicitly enable them with the @option{-fno-cx-limited-range} option. @ignore @c ===================================================================== @node More on Decimal Floating-Point Arithmetic @section More on Decimal Floating-Point Arithmetic @cindex decimal floating-point arithmetic @cindex floating-point arithmetic, decimal Proposed extensions to the C programming language call for the inclusion of decimal floating-point arithmetic, which handles floating-point numbers with a specified radix of 10, instead of the unspecified traditional radix of 2. The proposed new types are @code{_Decimal32}, @code{_Decimal64}, and @code{_Decimal128}, with corresponding literal constant suffixes of @code{df}, @code{dd}, and @code{dl}, respectively. For example, a 32-bit decimal floating-point variable could be defined as: @example _Decimal32 foo = 42.123df; @end example We stated earlier (@pxref{decimal-significand}) that the significand in decimal floating-point arithmetic is an integer, rather than fractional, value. Decimal instructions do not normally alter the exponent by normalizing nonzero significands to remove trailing zeros. That design feature is intentional: it allows emulation of the fixed-point arithmetic that has commonly been used for financial computations. One consequence of the lack of normalization is that there are multiple representations of any number that does not use all of the significand digits. Thus, in the 32-bit format, the values @code{1.DF}, @code{1.0DF}, @code{1.00DF}, @dots{}, @code{1.000_000DF}, all have different bit patterns in storage, even though they compare equal. Thus, programmers need to be careful about trailing zero digits, because they appear in the results, and affect scaling. For example, @code{1.0DF * 1.0DF} evaluates to @code{1.00DF}, which is stored as @code{100 * 10**(-2)}. In general, when you look at a decimal expression with fractional digits, you should mentally rewrite it in integer form with suitable powers of ten. Thus, a multiplication like @code{1.23 * 4.56} really means @code{123 * 10**(-2) * 456 * 10**(-2)}, which evaluates to @code{56088 * 10**(-4)}, and would be output as @code{5.6088}. Another consequence of the decimal significand choice is that initializing decimal floating-point values to a pattern of all-bits-zero does not produce the expected value @code{0.}: instead, the result is the subnormal values @code{0.e-101}, @code{0.e-398}, and @code{0.e-6176} in the three storage formats. GNU C currently supports basic storage and manipulation of decimal floating-point values on some platforms, and support is expected to grow in future releases. @c ??? Suggest chopping the rest of this section, at least for the time @c ??? being. Decimal floating-point support in GNU C is not yet complete, @c ??? and functionality discussed appears to not be available on all @c ??? platforms, and is not obviously documented for end users of GNU C. --TJR The exponent in decimal arithmetic is called the @emph{quantum}, and financial computations require that the quantum always be preserved. If it is not, then rounding may have happened, and additional scaling is required. The function @code{samequantumd (x,y)} for 64-bit decimal arithmetic returns @code{1} if the arguments have the same exponent, and @code{0} otherwise. The function @code{quantized (x,y)} returns a value of @var{x} that has been adjusted to have the same quantum as @var{y}; that adjustment could require rounding of the significand. For example, @code{quantized (5.dd, 1.00dd)} returns the value @code{5.00dd}, which is stored as @code{500 * 10**(-2)}. As another example, a sales-tax computation might be carried out like this: @example decimal_long_double amount, rate, total; amount = 0.70DL; rate = 1.05DL; total = quantizedl (amount * rate, 1.00DL); @end example @noindent Without the call to @code{quantizedl}, the result would have been @code{0.7350}, instead of the correctly rounded @code{0.74}. That particular example was chosen because it illustrates yet another difference between decimal and binary arithmetic: in the latter, the factors each require an infinite number of bits, and their product, when converted to decimal, looks like @code{0.734_999_999@dots{}}. Thus, rounding the product in binary format to two decimal places always gets @code{0.73}, which is the @emph{wrong} answer for tax laws. In financial computations in decimal floating-point arithmetic, the @code{quantized} function family is expected to receive wide use whenever multiplication or division change the desired quantum of the result. The function call @code{normalized (x)} returns a value that is numerically equal to @var{x}, but with trailing zeros trimmed. Here are some examples of its operation: @multitable @columnfractions .5 .5 @headitem Function Call @tab Result @item normalized (+0.00100DD) @tab +0.001DD @item normalized (+1.00DD) @tab +1.DD @item normalized (+1.E2DD) @tab +1E+2DD @item normalized (+100.DD) @tab +1E+2DD @item normalized (+100.00DD) @tab +1E+2DD @item normalized (+NaN(0x1234)) @tab +NaN(0x1234) @item normalized (-NaN(0x1234)) @tab -NaN(0x1234) @item normalized (+Infinity) @tab +Infinity @item normalized (-Infinity) @tab -Infinity @end multitable @noindent The NaN examples show that payloads are preserved. Because the @code{printf} and @code{scanf} families were designed long before IEEE 754 decimal arithmetic, their format items do not support distinguishing between numbers with identical values, but different quanta, and yet, that distinction is likely to be needed in output. The solution adopted by one early library for decimal arithmetic is to provide a family of number-to-string conversion functions that preserve quantization. Here is a code fragment and its output that shows how they work. @example decimal_float x; x = 123.000DF; printf ("%%He: x = %He\n", x); printf ("%%Hf: x = %Hf\n", x); printf ("%%Hg: x = %Hg\n", x); printf ("ntosdf (x) = %s\n", ntosdf (x)); %He: x = 1.230000e+02 %Hf: x = 123.000000 %Hg: x = 123 ntosdf (x) = +123.000 @end example @noindent The format modifier letter @code{H} indicates a 32-bit decimal value, and the modifiers @code{DD} and @code{DL} correspond to the two other formats. @c ===================================================================== @end ignore @node Round-Trip Base Conversion @section Round-Trip Base Conversion @cindex round-trip base conversion @cindex base conversion (floating point) @cindex floating-point 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: @multitable @columnfractions .5 .1 .1 .1 .1 .1 @item binary in @tab 24 @tab 53 @tab 64 @tab 113 @tab 237 @item decimal out @tab 9 @tab 17 @tab 21 @tab 36 @tab 73 @end multitable @multitable @columnfractions .5 .1 .1 .1 .1 @item decimal in @tab 7 @tab 16 @tab 34 @tab 70 @item binary out @tab 25 @tab 55 @tab 114 @tab 234 @end multitable We can compute the table numbers with these two functions: @example int matula(int nbits) @{ /* @r{Return output decimal digits needed for nbits-bits input.} */ return ((int)ceil((double)nbits / log2(10.0) + 1.0)); @} int goldberg(int ndec) @{ /* @r{Return output bits needed for ndec-digits input.} */ return ((int)ceil((double)ndec / log10(2.0) + 1.0)); @} @end example 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. @c ===================================================================== @node Further Reading @section 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: @c ===================================================================== @c Each bibliography item has a sort key, so the bibliography can be @c sorted in emacs with M-x sort-paragraphs on the region with the items. @c ===================================================================== @itemize @bullet @item @c sort-key: Abbott Paul H. Abbott and 15 others, @cite{Architecture and software support in IBM S/390 Parallel Enterprise Servers for IEEE Floating-Point arithmetic}, IBM Journal of Research and Development @b{43}(5/6) 723--760 (1999), @uref{https://doi.org/10.1147/rd.435.0723}. This article gives a good description of IBM's algorithm for exact decimal-to-binary conversion, complementing earlier ones by Clinger and others. @item @c sort-key: Beebe Nelson H. F. Beebe, @cite{The Mathematical-Function Computation Handbook: Programming Using the MathCW Portable Software Library}, Springer (2017), ISBN 3-319-64109-3 (hardcover), 3-319-64110-7 (e-book) (xxxvi + 1114 pages), @uref{https://doi.org/10.1007/978-3-319-64110-2}. 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 @cite{NIST Handbook of Mathematical Functions}, Cambridge (2018), ISBN 0-521-19225-0. See @uref{http://www.math.utah.edu/pub/mathcw} for compilers and libraries. @item @c sort-key: Clinger-1990 William D. Clinger, @cite{How to Read Floating Point Numbers Accurately}, ACM SIGPLAN Notices @b{25}(6) 92--101 (June 1990), @uref{https://doi.org/10.1145/93548.93557}. See also the papers by Steele & White. @item @c sort-key: Clinger-2004 William D. Clinger, @cite{Retrospective: How to read floating point numbers accurately}, ACM SIGPLAN Notices @b{39}(4) 360--371 (April 2004), @uref{http://doi.acm.org/10.1145/989393.989430}. Reprint of 1990 paper, with additional commentary. @item @c sort-key: Goldberg-1967 I. Bennett Goldberg, @cite{27 Bits Are Not Enough For 8-Digit Accuracy}, Communications of the ACM @b{10}(2) 105--106 (February 1967), @uref{http://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. @item @c sort-key: Goldberg-1991 David Goldberg, @cite{What Every Computer Scientist Should Know About Floating-Point Arithmetic}, ACM Computing Surveys @b{23}(1) 5--58 (March 1991), corrigendum @b{23}(3) 413 (September 1991), @uref{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. @item @c sort-key: Juffa Norbert Juffa and Nelson H. F. Beebe, @cite{A Bibliography of Publications on Floating-Point Arithmetic}, @uref{http://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. @item @c sort-key: Kahan William Kahan, @cite{Branch Cuts for Complex Elementary Functions, or Much Ado About Nothing's Sign Bit}, (1987), @uref{http://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, @cite{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. @item @c sort-key: Knuth Donald E. Knuth, @cite{A Simple Program Whose Proof Isn't}, in @cite{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, @uref{https://doi.org/10.1007/978-1-4612-4476-9}. 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. @item @c sort-key: Matula-1968a David W. Matula, @cite{In-and-out conversions}, Communications of the ACM @b{11}(1) 57--50 (January 1968), @uref{https://doi.org/10.1145/362851.362887}. @item @c sort-key: Matula-1968b David W. Matula, @cite{The Base Conversion Theorem}, Proceedings of the American Mathematical Society @b{19}(3) 716--723 (June 1968). See also other papers here by this author, and by I. Bennett Goldberg. @item @c sort-key: Matula-1970 David W. Matula, @cite{A Formalization of Floating-Point Numeric Base Conversion}, IEEE Transactions on Computers @b{C-19}(8) 681--692 (August 1970), @uref{https://doi.org/10.1109/T-C.1970.223017}. @item @c sort-key: Muller-2010 Jean-Michel Muller and eight others, @cite{Handbook of Floating-Point Arithmetic}, Birkh@"auser-Boston (2010), ISBN 0-8176-4704-X (xxiii + 572 pages), @uref{https://doi.org/10.1007/978-0-8176-4704-9}. 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. @item @c sort-key: Muller-2018 Jean-Michel Muller and eight others, @cite{Handbook of Floating-Point Arithmetic}, Second edition, Birkh@"auser-Boston (2018), ISBN 3-319-76525-6 (xxv + 627 pages), @uref{https://doi.org/10.1007/978-3-319-76526-6}. This is a new edition of the preceding entry. @item @c sort-key: Overton Michael Overton, @cite{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), @uref{http://www.ec-securehost.com/SIAM/ot76.html}. This is a small volume that can be covered in a few hours. @item @c sort-key: Steele-1990 Guy L. Steele Jr. and Jon L. White, @cite{How to Print Floating-Point Numbers Accurately}, ACM SIGPLAN Notices @b{25}(6) 112--126 (June 1990), @uref{https://doi.org/10.1145/93548.93559}. See also the papers by Clinger. @item @c sort-key: Steele-2004 Guy L. Steele Jr. and Jon L. White, @cite{Retrospective: How to Print Floating-Point Numbers Accurately}, ACM SIGPLAN Notices @b{39}(4) 372--389 (April 2004), @uref{http://doi.acm.org/10.1145/989393.989431}. Reprint of 1990 paper, with additional commentary. @item @c sort-key: Sterbenz Pat H. Sterbenz, @cite{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 @emph{before} the introduction of IEEE 754 arithmetic. @end itemize