The Perils of Floating PointbyBruceM. BushCopyright (c) 1996 Lahey Computer Systems, Inc.Permission to copy is granted with acknowledgement of the source.
Many great engineering and scientific advances of recent decades would not have been possible without the floating-point capabilities of digitalcomputers. Still, some results of floating-point calculations lookpretty strange, even to people with years of mathematical experience.I will attempt to explain the causes of some of these strange resultsand give some suggestions where appropriate.Floating-point representations and arithmetic are inexact, but I don'tbelieve that is particularly troublesome to most programmers. Many inputvalues are measurements, which are inherently inexact, so the questionabout the output values isn't whether there is error, but how much errorshould be expected. However, when you can compute a more accurate resultin your head than your computer can with its floating-point, you start toget suspicious. I have programmed my examples in FORTRAN for a couple of reasons:
1. More floating-point calculations are performed in FORTRAN than anyother computer language.
2. I work for a company, Lahey Computer Systems, that develops and sellsFORTRAN language systems.
Binary Floating-PointAt the heart of many strange results is one fundamental: floating-pointon computers is usually base 2, whereas the external representation isbase 10. We expect that 1/3 will not be exactly representable, but itseems intuitive that .01 would be. Not so! .01 in IEEE single-precisionformat is exactly 10737418/1073741824 or approximately 0.009999999776482582. You might not even notice this difference until you see a bit of code like the following:REAL XDATA X /.01/IF ( X * 100.d0 .NE. 1.0 ) THEN PRINT *, 'Many systems print this surprising result. 'ELSE PRINT *, 'And some may print this.'ENDIF Base-10 floating-point implementations don't have this anomaly. However,base-10 floating-point implementations are rare because base-2 (binary)arithmetic is so much faster on digital computers.
InexactnessFloating-point arithmetic on digital computers is inherently inexact.The 24 bits (including the hidden bit) of mantissa in a 32-bitfloating-point number represent approximately 7 significant decimaldigits. Unlike the real number system, which is continuous, afloating-point system has gaps between each number. If a numberis not exactly representable, then it must be approximated by oneof the nearest representable values.Because the same number of bits are used to represent all normalizednumbers, the smaller the exponent, the greater the density ofrepresentable numbers. For example, there are approximately 8,388,607single-precision numbers between 1.0 and 2.0, while there are onlyabout 8191 between 1023.0 and 1024.0. On any computer, mathematically equivalent expressions can producedifferent values using floating-point arithmetic. In the followingexample, Z and Z1 will typically have different values because (1/Y)or 1/7 is not exactly representable in binary floating-point: REAL X, Y, Y1, Z, Z1DATA X/77777/, Y/7/Y1 = 1 / YZ = X / YZ1 = X * Y1IF (Z .NE. Z1) PRINT *, 'Not equal!'END
Insignificant DigitsThe following code example illustrates the phenomenon of meaninglessdigits that could seem to be significant:REAL A, YDATA Y /1000.2/ ! About 7 digits of precision in YA = Y - 1000.0 ! About 3 significant digits in resultPRINT *, A ! Prints 0.200012END A single-precision (REAL) entity can represent a maximum of about 7decimal digits of precision, so the subtraction above represents(1000.200 - 1000.000). The result, therefore, can only representabout 3 decimal digits. The program, however, will happily printout "0.200012". Because 1000.2 is not exactly representable inbinary floating-point and 1000.0 is, the result A is a little largerthan 0.2. The computer doesn't know that the digits beyond ".200"have no meaning.Perhaps someday the computer will keep track of the number of bitsin a result that are truly significant. For now, it is still theresponsibility of the programmer. If you stay aware of the numberof decimal digits represented by a data type, approximating thenumber of significant digits is a straight-forward, but perhapstime-consuming, task. Give the most attention to: 1. subtractions of numbers that are nearly equal,
2. additions of numbers whose magnitudes are nearly equal, butwhose signs are opposite, and
3. additions and subtractions of numbers that differ greatly in magnitude.
Crazy ConversionsConversions to integer can unmask inaccuracies in a floating-pointnumber, as is demonstrated by the next example. The closestsingle-precision floating-point number to 21.33 is slightly less than21.33, so when it is multiplied by 100., the result Y is slightly lessthan 2133.0. If you print Y in a typical floating-point format,rounding causes it to be displayed as 2133.00. However, if you assignY to an integer I, no rounding is done, and the number istruncated to 2132.REAL X, YINTEGER IX = 21.33 ! Slightly less than 21.33Y = X * 100. ! Slightly less than 2133.0I = Y ! Truncates to 2132PRINT *, Y, I ! Prints "2133.00 2132"END The following program prints "1.66661000251770" when compiled withLahey's LF90:DOUBLE PRECISION DREAL XX = 1.66661 ! Assign to single precisionD = X ! Convert to double precisionPRINT *, DEND You ask, "Why do you extend the single-precision number with theseemingly random '000251770'?" Well, the number isn't extendedwith random values; the computer's floating-point does theconversion by padding with zeros in the binary representation.So D is exactly equal to X, but when it is printed out to 15decimal digits, the inexactness shows up. This is also anotherexample of insignificant digits. Remember that assigning asingle-precision number to a double-precision number doesn'tincrease the number of significant digits.
Too Many DigitsYou may decide to check the previous program example by printing outboth D and X in the same format, like this:30 FORMAT (X, 2G25.15)PRINT 30, X, D In some FORTRAN implementations, both numbers print out the same.You may walk away satisfied, but you are actually being misled bylow-order digits in the display of the single-precision number.In Lahey FORTRAN the numbers are printed out as: 1.66661000000000 1.66661000251770 The reason for this is fairly simple: the formatted-I/O conversionroutines know that the absolute maximum decimal digits that have anysignificance when printing a single-precision entity is 9. The restof the field is filled with the current "precision-fill" character,which is "0" by default. The precision-fill character can be changedto any ASCII character, e.g., asterisk or blank. Changing theprecision-fill character to "*" emphasizes the insignificance of thelow-order digits: 1.66661000****** 1.66661000251770
Too Much PrecisionThe IEEE single-precision format has 24 bits of mantissa, 8 bits ofexponent, and a sign bit. The internal floating-point registers inIntel microprocessors such as the Pentium have 64 bits of mantissa,15 bits of exponent and a sign bit. This allows intermediate calculationsto be performed with much less loss of precision than many otherimplementations. The down side of this is that, depending upon howintermediate values are kept in registers, calculations that look thesame can give different results.In the following example, a compiler could generate code thatcalculates A/B, stores the intermediate result into a single-precisiontemporary, calculates X/Y, performs a reversed subtract of the temporary,then stores the result. Z will not be zero, because precision will belost in storing into a single-precision temporary. If the generated codekeeps the intermediate result in registers, no precision will be lost,and Z will be zero. REAL A, B, X, Y, ZDATA A/10./, B/3./, X/10./, Y/3./Z = (A/B) - (X/Y)PRINT *, Z ! Could be zero or not.END The next example illustrates a variation on the previous example.A compiler can still generate code that keeps the intermediate resultC in a register, which means that Z will be zero. If precision islost by storing A/B into C, then Z will be nonzero.REAL A, B, C, X, Y, ZDATA A/10./, B/3./, X/10./, Y/3./C = A/BZ = C - (X/Y)PRINT *, Z ! Could be zero or not.END The slight variation of adding the statement label 100 foils theoptimization of keeping C in a register, so Z will probably be nonzerowith almost any compiler.REAL A, B, C, X, Y, ZDATA A/10./, B/3./, X/10./, Y/3./C = A/B100 Z = C - (X/Y)PRINT *, ZIF ( ... ) GO TO 100END
Safe ComparisonsDifferent computers use different numbers of bits to storefloating-point numbers. Even when the same IEEE formats are usedfor storing numbers, differences in calculations can occur becauseof the size of intermediate registers. To increase portability andto ensure consistent results, I recommend against comparing for exactequality of real numbers in FORTRAN. A better technique is to comparethe absolute value of the difference of two numbers with an appropriateepsilon to get relations like approximately equal, definitely greaterthan, etc.
Example:REAL EPSILON, X, A, BPARAMETER (EPSILON = .000001)DATA A/13.9/, B/.000005/X = (A * B) / BIF (ABS(X - A) .LE. (ABS(X)*EPSILON)) THEN PRINT *, 'X is approximately equal to A'ENDIFIF ((X - A) .GT. (ABS(X)*EPSILON)) THEN PRINT *, 'X is definitely greater than A'ENDIFIF ((A - X) .GT. (ABS(X)*EPSILON)) THEN PRINT *, 'X is definitely less than A'ENDIFEND Multiplying the epsilon by one of the comparands adjusts the comparisonto the range of the numbers, allowing a single epsilon to be used formany, or perhaps all compares. For the most predictable results, usean epsilon half as large and multiply it by the sum of the comparands,as in the following example:REAL EPSILON, X, A, BPARAMETER (EPSILON = .0000005) ! Smaller epsilonDATA A/13.9/, B/.000005/X = (A * B) / BIF (ABS(X - A) .LE. (ABS(X+A)*EPSILON)) THEN PRINT *, 'X is approximately equal to A'ENDIF Even comparisons of greater-than, less-than-or-equal-to, etc., canproduce unexpected results, because a floating-point computation canproduce a value that is not mathematically possible. In the followingexample X is always mathematically greater than J, so X/J should alwaysbe greater than 1.0. For large values of J, however, the addition ofdelta is not representable by X, because of limited mantissa size.REAL X, DELTADATA DELTA/.001/DO 10 J = 1, 100 000 X = J + DELTA ! Make X bigger than J CALL SUB (X) ! Force X out of register IF ( X/J .LE. 1.0 ) THEN ! X/J always > 1 ? PRINT *, 'Error !' STOP END IF10 CONTINUEENDSUBROUTINE SUB (X)END
Programming with the PerilsThere are no easy answers. It is the nature of binary floating-pointto behave the way I have described. In order to take advantage of thepower of computer floating-point, you need to know its limitations andwork within them. Keeping the following things in mind when programmingfloating-point arithmetic should help a lot:1. Only about 7 decimal digits are representable in single-precisionIEEE format, and about 16 in double-precision IEEE format.
2. Every time numbers are transferred from external decimal to internalbinary or vice-versa, precision can be lost.
3. Always use safe comparisons.
4. Beware of additions and subtractions that can quickly erode thetrue significance in a result. The computer doesn't know what bitsare truely significant.5. Conversions between data types can be tricky. Conversions to double-precision don't increase the number of truely significant bits. Conversions to integer always truncate toward zero, even if the floating-point number is printed as a larger integer.
6. Don't expect identical results from two different floating-pointimplementations. I hope that I have given you a little more awareness of what ishappening in the internals of floating-point arithmetic, and thatsome of the strange results you have seen make a little more sense. While some of the "perils" can be avoided, many just need to be understood and accepted.
IEEE Standard Floating-Point FormatsThe Institute of Electrical and Electronics Engineers, Inc. (IEEE) hasdefined standards for floating-point representations and computationalresults (IEEE Std 754-1985). This section is an overview of the IEEEstandard for representing floating-point numbers. The data containedherein helps explain some of the details in the rest of the article,but is not required for understanding the basic concepts.Most binary floating-point numbers can be represented as 1.ffffff x 2^n,where the 1 is the integer bit, the f's are the fractional bits, and then is the exponent. The combination of the integer bit and the fractionalbits is called the mantissa (or significand). Because most numbers can havetheir exponent adjusted so that there is a 1 in the integer bit (a processcalled normalizing), the 1 does not need to be stored, effectively allowingfor an extra bit of precision. This bit is called a hidden bit. Numbersare represented as sign-magnitude, so that a negative number has the samemantissa as a positive number of the same magnitude, but with a sign bitof 1. A constant, called a bias, is added to the exponent so that allexponents are positive. The value 0.0, represented by a zero exponent and zero mantissa, can havea negative sign. Negative zeros have some subtle properties that willnot be evident in most programs. A zero exponent with a nonzero mantissais a "denormal." A denormal is a number whose magnitude is too small tobe represented with an integer bit of 1 and can have as few as onesignificant bit. Exponent fields of all ones (largest exponent) represent special numericresults. A mantissa of zero represents infinity (positive or negative);a nonzero mantissa represents a NAN (not-a-number). NANs, which occur asa result of invalid numeric operations, are not discussed further in thisarticle. The IEEE Standard defines 32-bit and 64-bit floating-point representations.The 32-bit (single-precision) format is, from high-order to low-order, asign bit, an 8-bit exponent with a bias of 127, and 23 bits of mantissa.The 64-bit (double-precision) format is, a sign bit, an 11-bit exponentwith a bias of 1023, and 52 bits of mantissa. With the hidden bit,normalized numbers have an effective precision of 24 and 53 bits,respectively. Single-precision format
31, 30-23, 22-0
S, Exponent, Significand Double-precision format
63, 62-52, 51-0
S, Exponent, Significand
BibliographyAmerican National Standards Institute (1978), "American NationalStandard, Programming Language FORTRAN", ANSI X3.9-1978, ISO1539-1980 (E).IEEE Computer Society (1985), "IEEE Standard for Binary Floating-PointArithmetic", IEEE Std 754-1985.
|
Formed in 2009, the Archive Team (not to be confused with the archive.org Archive-It Team) is a rogue archivist collective dedicated to saving copies of rapidly dying or deleted websites for the sake of history and digital heritage. The group is 100% composed of volunteers and interested parties, and has expanded into a large amount of related projects for saving online and digital history.
