As with numerical analysis, the intricacies of how floating-point numbers are represented, and the quirks of their representation on different platforms, are a maelstrom into which you can fall, dazed, confused and unproductive. However (and again as with numerical analysis), if your codes become elaborate enough, then you are going to have to take the plunge.
Even if you have no plans to venture into the deep, there is a minimum awareness of the problems which can help you write more robust code.
What follows here is rather detailed, but it does, I believe, represent the majority of what you might need to know about floating point representation; there is a good alternative introduction in Numerical Recipes [nr], and an excellent and detailed introduction in chapter 2 of Sun's Numerical Computation Guide [sunncg] (appendix E of the same book claims to be `What every computer scientist should know about floating point arithmetic', and is an edited reprint of [goldberg91]. It makes interesting reading, but it's probably more than many natural scientists need to know). My account here is not intended to supplant either of these sources. Below, I'll talk exclusively of IEEE base-2 floating point numbers, as defined in IEEE standard IEEE-754; there are other standards, but they're now of historical interest only, as just about all modern machines other than VAXes and older Crays use IEEE. Most of what I say concerning accuracy will be about single precision; exactly the same issues arise with double precision, but you can sweep them under the carpet for longer.
The IEEE standard leaves it open to chip manufacturers to decide the order in memory of the four or eight bytes of a floating point number (the same is true of integers, though they're not part of the IEEE spec). This is denoted by the endian-ness (or `bytesex') of the machine. Alpha and Intel chips are little-endian; Sparcs, the Motorola 68k family used in Macs, and Motorola PowerPC chips, are big-endian (there's also a `PGP-endian' ordering, but that really isn't likely to trouble anyone any more). This generally does not matter if you stick to a single platform (other than in the context of programs such as those described in Appendix A.1 and Appendix A.2), but it makes a difference if you want to share raw data files (see Section 18.104.22.168) between architectures.
The central insight is that you cannot represent an infinite number of reals in a finite number of bits without approximation. It is a consequence of this that the result of any calculation will lose precision when it is represented as a floating-point bit pattern, and that the extent and importance of this loss of precision depends on the operands. Both of these seem obvious by themselves, but the consequences can be sufficiently non-obvious to trip you up.
A non-special IEEE floating-point single-precision
longword represents a number
where (the sign), (the mantissa or significand) and (the exponent) are base-2 numbers represented in 1, 8 and 23 bits respectively, and ; the offset 127 is called the bias. These are encoded into the 32 bits of the longword as shown in table Table 1.
|0||0 1 1 1 1 1 1 1||0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0|
|0||1 0 0 0 0 0 0 0||0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0|
|0||1 0 0 0 0 0 0 0||1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0|
|0||0 1 1 0 1 0 0 0||0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0|
|0||0 1 1 1 1 1 1 1||0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1|
|0||0 1 1 1 1 1 1 1||0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 1 1 0 0 1 1|
|0||0 1 1 1 0 0 1 0||1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0|
IEEE floating-point representations of numbers. This table is adapted from table 2--3 in [sunncg] and figure 1.3.1 in [nr], which have fuller discussions of the issues here. For discussion of and see the text.
Firstly, it is important to distinguish between the range and the accuracy of floating-point numbers. The smallest and largest positive numbers which can be represented by normal IEEE single-precision numbers are and , but this is very different from the accuracy to which these numbers are represented. These extreme values differ from the next representable ones up and down by and respectively, so that any number which differs from them by less than that amount is unrepresentable. This accuracy limitation is expressed by the machine epsilon, (for IEEE single-precision), which is the maximum relative error incurred when a real number is represented in floating-point. It is a consequence of this that is the smallest number such that , where the operator represents addition of floating-point numbers.[Note 1]
Another way of thinking about this is as follows. When two floating-point numbers are added, the exponent of the smaller one must be increased to match that of the larger, and the mantissa (or significand), , must be shifted rightwards to compensate. If the two numbers differ greatly in magnitude, then the smaller will lose significance in this operation, and if it is less than a factor of the larger one, the shift will push the number right off the end of the mantissa, and contribute nothing to the sum. You can easily confirm that .
The immediate practical consequence of this is that if, in the bowels of some calculation, you are adding up many small numbers (doing an inner product for very large vectors, for example), then the final total may be grossly inaccurate if you've fallen into this trap. Similarly, if you are subtracting numbers which are almost equal, perhaps in the course of debiasing or removing some DC offset, the subtraction may allow all or most of the leading accurate digits to cancel, leaving the result to be determined by digits possibly heavily contaminated by roundoff. [RW]
Consider the simple calculation , with , and . Here, the rounding in the two multiplications conspires to give a difference which has a huge relative error of . We can, however, address the problem by rewriting the calculation in ways which are equivalent for real arithmetic, but inequivalent for machine arithmetic. If we do the subtraction before the multiplication, and calculate , then we retain a little precision in the subtraction and end up with a relative error of . Finally, if we instead calculate , then we do the subtraction with a full 23 bits of precision (compare Table 1), and lose as little as possible of this in the final multiplication, and end up with a relative error of only . This is clearly a rather synthetic example, but it is not at all unrealistic, as it is equally clearly closely related to the common problem of calculating the discriminant of the quadratic formula, when the quadratic is nearly degenerate.
These problems are demonstrated in
the short program
fpdemo.c in Appendix A.3. Note that the numbers , and are (a)
specifically chosen to demonstrate the effect, and (b)
are calculated within the program, rather than being
initialised from the decimal representations 1.000244,
and so on. If the latter were not true, then the
improvement in relative error would disappear, since all
precision would have been lost in this initialisation.
If you wish to explore the representation of
floating-point numbers on your machine, you can do so
using the example program
fpp.c in Appendix A.2.
These problems are unfortunately both easy to run into, and hard to notice if you're not in the habit of looking for them whilst writing your code. A brute-force way around them is to do sensitive parts of your calculation in double precision. If you are doing the calculation in double precision and still running into the problem, then you will have to rearrange the calculation somehow, to contain the loss of precision.[Note 2] Note, however, that a compiler's optimizer can frustrate your efforts here, if you've given it a sufficiently high optimization level that it starts reorganising your code: in general , and , but a high optimization level may instruct the compiler to ignore this fact, which can be disastrous for the accuracy of your code. Moral: make sure you know what the various optimization levels actually do before you invoke them.
As alluded to above, it is generally possible to
sweep accuracy problems under the carpet by using
double-precision floats. These occupy eight bytes of
storage rather than four, and the analogue of Eqn.(1) is
where , and are respectively 1, 52 and 11 bits long, and the bias is 1023. Thus the smallest and largest normal positive double-precision IEEE numbers are and , and the machine epsilon for double precision is .
On some platforms such as Solaris and Compaq, with the manufacturer's compilers, it is possible to use quadruple precision. You should not use double or quadruple precision automatically, however, for a variety of practical and principled reasons.
Firstly, if your program already suffers from rather prodigal use of memory (ie, has very large arrays), then the doubling in the size of real arrays will only make the problem worse and, by probably increasing the amount of swapping, might make your program significantly slower.
Secondly, although there are some issues to do with the relative speed of single- and double-precision calculations, these are probably not significant enough to be worth worrying about in all but the most long-running codes, and in any case would generally be swept up by the compiler (I'd welcome any correction or amplification on this point).[Note 3]
Thirdly, if your results change if you switch to double-precision, then there is an accuracy problem in your code, which would probably bear some investigation. If there is some part of your code which genuinely requires the extra precision -- say because you have to add numbers of hugely different scale in accumulating a large inner product -- then there is not much you can do about it, and that part of your code, at least, must have the extra precision. Failing such an explanation, however, you can regard such behaviour as a miner's canary, indicating some potential problem with your code, which you might benefit from investigating further.
Fourthly, remember that some algorithms
are designed to work with single precision, and will no
longer work properly if you search-and-replace
real*4. As just
one example, one of the Numerical Recipes algorithms
svdcmp) includes a convergence test of,
if (val+err==val). Whatever
the merits or otherwise of this as a test, it is a cheap
way of establishing whether
err has reached
machine precision, which will not work in the expected
err are double
For further details on floating point representation, see, for example [hauser96], and informal but very useful discussions (both available on the web) by William Kahan [kahan96] and Jim Demmel [demmel]. The former is one of the authors of the IEEE-754 standard for floating-point numbers.
As well as defining the format in which normal floating-point numbers are stored, the IEEE-754 standard includes the definitions for several other (classes of) special values.
When the exponent (of a single-precision number) is zero,
the longword is taken to represent the number
(compare Eqn.(1)). Unlike the usual floating-point numbers, which have an implied leading 1 in the significand and 23 bits of precision, and which are referred to as `normalised', these have an implied leading 0 and less than 23 bits of precision. These are `denormal' or `subnormal' numbers, and are the result of an underflow, such as dividing the smallest normalised number by two. This behaviour is known as `gradual underflow', and allows the precision in a calculation to degrade gracefully when it underflows. It is distinct from `Store-0 underflow' common before the IEEE standard, in which any expression which underflowed was replaced by zero. This was, and remains, one of the more contentious parts of the IEEE standard. Be warned that older Crays, which use their own floating-point format, have a Store-0 underflow policy, and that the Alpha chip, although it generally implements IEEE floats, has a Store-0 underflow as its default mode, and will neither produce, nor accept as input, denormalised numbers.
If the significand as well as the exponent is zero, the longword represents the number . Note that the IEEE zero is signed, which allows and to be positive and negative infinity; this can be important when doing calculations near branch points.
If the exponent is 255 and the significand is zero, the longword represents the value . That is, infinity has a special value, distinct from the largest possible normalised number. Positive infinity is generated by, for example , or , where infinity is the mathematically correct, but otherwise unrepresentable, value of a calculation.
If the exponent is 255 and the significand is non-zero, the
longword is Not a Number, usually represented by the string
`NaN'. A NaN is the result of an operation on invalid operands, such
as or . NaNs have the properties that any
operation which has a NaN as an operand has a NaN as its result; and
any comparison on a NaN, such as
evaluates to False, including
NaN==NaN. The only exception
is that, when
x is a NaN, then
x != x is
true. Why would you
want to use such a peculiar number? Generally you don't, and its
appearance in a calculation is an error, which is why processors can
be set to dump core when a NaN or an Infinity is produced. However,
it can be useful to turn off this behaviour if it is not off by
default, and rather than elaborately avoid producing a NaN at each
stage, make a check once at the end of a calculation, possibly
invoking an alternative (possibly more expensive or special-cased)
algorithm if any NaNs are found.
Different compilers handle exceptional values in different ways; see Table 2. Of the three Starlink platforms, only the Alpha traps on exceptional values by default.
|Compiler||traps enabled by default||see||suppress with|
|Sun cc & f77||none|
|Alpha cc||overflow, divide-by-zero, invalid operation||-fptm, -ieee||-ieee|
|Alpha f77||overflow, divide-by-zero, invalid operation||-fpe||-fpe1|
Summary of compiler IEEE traps
We have only scratched the surface of a rather
intricate topic here. Both Suns and Alphas have
man-pages, which hint at the broad range of
floating-point options available when compiling and
linking your code. See Section 22.214.171.124 for
discussion of the efficiency tradeoffs when using IEEE