Floating COMP3320/COMP6464: Point Numbers Alistair Rendell
COMP3320 Lecture 3-0 Copyright c 2008 The Australian National University
3.1
Rational
It is extremely important that all scientic programmers have a basic understanding of numerical data representation and error propagation catastropic failure of the Ariane 5 rocket in 1996 and the inability of Patriot missile systems to reach their targets during the 1991 Gulf war were both attributed to numerical computing errors Ariane 5 Explosion (1996) - data conversion of too large number (trying to stu a 64-bit number into a 16-bit eld) - problem came to light as Ariane 5 was faster than Ariane 4.
Patriot-Scud fails again (1991) Mars Orbiter loss (1999)
- due to inaccurate calculation of the time since computer boot time
- mixture of miles and meters
See http://wwwzenger.informatik.tu-muenchen.de/persons/huckle/bugse.html
COMP3320 Lecture 3-1 Copyright c 2008 The Australian National University
3.2
Floating Point Numbers
Most scientic problems are concerned with continuous phenomena, e.g time distance, velocity, temperature - i.e. real numbers not integers Computers operate in an exact discrete world
An understanding of how oating point arithmetic is handled on computers is fundamental to scientic computing - all computers (CPUs) of interest have separate oating point and integer arithmetic units Numerical analysis deals with the design and analysis of the behaviour of algorithms for solving scientic problems
COMP3320 Lecture 3-2 Copyright c 2008 The Australian National University
3.3
Sources of Error
Before computation begins
- modeling - empirical measurements - previous computations - mistakes
During computation
- approximation error (truncation or discretization error) - rounding error
Accuracy of nal result inevitably reects combination of approximations, and perturbations may be amplied by mature of problem or algorithm
COMP3320 Lecture 3-3 Copyright c 2008 The Australian National University
3.4
Approximation Error
Dierence between true result (for actual input) and result that would be produced by given algorithm using exact arithmetic, e.g. - approximation of a series x2 x3 xn = 1 + x + + + e = n=0 n! 2 6
x
ex - numerical derivative
N is large n=0
xn n!
f (x + h) f (x) f (x) = lim h0 h f (x + h) f (x) f (x) h (h is small)
COMP3320 Lecture 3-4 Copyright c 2008 The Australian National University
3.5
Rounding Error
Dierence between result produced by given algorithm using exact aritmetic and result produced by same algorithm using rounded arithmetic - due to inexactness in representation of real numbers and arithmetic operations upon them - eects like 1/3 + 1/3 + 1/3 = 0.9999
COMP3320 Lecture 3-5 Copyright c 2008 The Australian National University
3.6
Total Computation Error
Total Computation Error = Approximation Error + Rounding Error - but usually one of these is dominant (see labs)
COMP3320 Lecture 3-6 Copyright c 2008 The Australian National University
3.7
Absolute and Relative Error Error Absolute = Value Approximate Value True Value Approximate = (Value True)(1 + Relative Error) Va = Vt(1 + )
Absolute Error
Relative Error Error Relative = (Error Absolute)/(Value True) (Absolute Error)/(Value Approximate) True value is usually unknown
- so we must estimate or bound error rather than compute it exactly
Likewise relative error is often taken to be relative to approximate value, rather than true value
COMP3320 Lecture 3-7 Copyright c 2008 The Australian National University
3.8
Floating Point Numbers
Characterised by
- : base or radix - t: the precision (signicant gures) - [L, U ]: exponent range
Number (e.g. 1.23456 1013) represented as d1 d2 d3 dt1 e x = (d0 + 1 + 2 + 3 + + t1 ) 0 di 1 i = 0, 1, 2, , t 1 LeU d0d1d2 dt1 mantissa (d1d2 dt1 behind decimal place) e is the exponent
COMP3320 Lecture 3-8 Copyright c 2008 The Australian National University
3.9
Typical Floating Point Systems
Most computers use binary ( = 2) arithmetic System t L U IEEE SP 2 24 -126 127 IEEE DP 2 53 -1022 1023 Cray XMP 2 48 -16383 16384 HP Calculator 10 12 -499 499 IBM Mainframe 16 6 -64 63 IEEE by far most common
Note trade-o between t and [L, U ]
COMP3320 Lecture 3-9 Copyright c 2008 The Australian National University
3.10
Normalization
Floating point system is normalized if leading digit d0 is always nonzero - mantissa m of nonzero oating point number always 1m< Why - unique representation of each number - no digits wasted on leading zero - in binary system the leading bit need not be stored
COMP3320 Lecture 3-10 Copyright c 2008 The Australian National University
3.11
Properties of Floating Point Systems
FP number system is nite and discrete Number of normalized FP numbers 2( 1) t1(U L + 1) + 1
Smallest and largest positive normalized number - underow level = UFL = l - overow level = OFL = U +1(1 t)
FP numbers are equally spaced only between successive powers of Not all real numbers are exactly representable - those that are, are called machine numbers
COMP3320 Lecture 3-11 Copyright c 2008 The Australian National University
3.12
IEEE Single Precision (32 Bits)
= 2, t = 24 (normalized), L = 126, U = +127 2(2)23(127+126+1)+1 numbers 1 2 3 4 Sign(1) Special Cases Exponent(8) 9 10 11 Mantissa(23) 32
- zero - innity (plus and minus) - not a number (nan) - sub-normal numbers
COMP3320 Lecture 3-12 Copyright c 2008 The Australian National University
3.13
Special Cases: Single Precision 0 0000 1110 1010 0011 1111 1011 0100 111 Valid IEEE Single Precision Number 0 0000 0000 0000 0000 0000 0000 0000 000 IEEE +0 1 0000 0000 0000 0000 0000 0000 0000 000 IEEE -0 0 1111 1111 0000 0000 0000 0000 0000 000 IEEE + 1 1111 1111 0000 0000 0000 0000 0000 000 IEEE 1111 1111 NON ZERO IEEE Not A Number (NAN)
COMP3320 Lecture 3-13 Copyright c 2008 The Australian National University
3.14
Example Floating Point System
4.0
= 2, t = 3(normalized), L = 1, U = +1 (25 numbers)
Sign d0 d1 d2 0 1 1 1 0 1 1 0 0 1 0 1 0 1 0 0 0 1 1 1 0 1 1 0 0 1 0 1 0 1 0 0 0 1 1 1 0 1 1 0 0 1 0 1 0 1 0 0 0 0 0 0 e decimal 1 3.5 1 3.0 1 2.5 1 2.0 0 1.75 0 1.50 0 1.25 0 1.0 -1 0.875 -1 0.75 -1 0.625 -1 0.5 0 0.0
Overflow limit 3.0
2.0
1.0 Underflow Limit 0.0 Underflow Limit 1.0
2.0
3.0 Overflow limit 4.0
OFL = (1.11)2 21 = 3.510, UFL = (1.00)2 21 = 0.510 Note grainy and unequally spaced
COMP3320 Lecture 3-14 Copyright c 2008 The Australian National University
3.15
Rounding Rules
Non-machine numbers must be approximated to nearby FP number - termed rounding and introduces rounding error - denote as f l(x)
Rounding rules
- chop: truncate x after (t 1) digit (round toward zero or round down) - round to nearest
IEEE uses rounding to nearest
COMP3320 Lecture 3-15 Copyright c 2008 The Australian National University
3.16
Machine Precision
Accuracy of FP system characterized by quantity called unit round-o error or machine precision, or machine epsilon, denoted by mach - chopping: = mach 1t - nearest: mach = 1/2 1t Maximum relative error in representing real number x is given by: f l(x) x x
mach
COMP3320 Lecture 3-16 Copyright c 2008 The Australian National University
3.17
Machine Precision = 224 107 53 1016 mach = 2
For IEEE FP systems - Single Precision: - Double Precision:
mach
IEEE SP and DP have 7 and 16 decimal digits of precision respectively Machine precision and UFL are VERY dierent 0 < UFL <
mach
< OFL
COMP3320 Lecture 3-17 Copyright c 2008 The Australian National University
3.18
Machine Precision: Toy System
mach
Round to nearest:
= 1/2 1t = 0.125 f l(x) x x
Exponent -1 0.500 -1 0.500 -1 0.500 -1 0.500 0 1.000 0 1.000 0 1.000 0 1.000 1 2.000 1 2.000 1 2.000 1 2.000
Err1 =
Err2 =
x 0.500 0.625 0.750 0.875 1.000 1.250 1.500 1.750 2.000 2.500 3.000 3.500 f l(x) 0.562 0.688 0.812 0.938 1.125 1.375 1.625 1.875 2.250 2.750 3.250
f l(x) x f l(x)
Err1 0.111 0.091 0.077 0.067 0.111 0.091 0.077 0.067 0.111 0.091 0.077 Err2 0.125 0.100 0.083 0.071 0.125 0.100 0.083 0.071 0.125 0.100 0.083
1 1 1 1 1 1 1 1 1 1 1 1
Mantissa 0 0 1.000 0 1 1.250 1 0 1.500 1 1 1.750 0 0 1.000 0 1 1.250 1 0 1.500 1 1 1.750 0 0 1.000 0 1 1.250 1 0 1.500 1 1 1.750
COMP3320 Lecture 3-18 Copyright c 2008 The Australian National University
3.19
Subnormal Numbers and Gradual Underow
Normalization causes gap around zero
4.0 3.0 2.0 1.0 0.0 1.0 2.0 3.0 4.0
Subnormal or denormalized numbers permit leading digit to be zero only for minimum exponent.
4.0 3.0 2.0 1.0 0.0 1.0 2.0 3.0 4.0
Subnormal Numbers
Subnormal numbers extend range of magnitude representable, but they have less precision and unit roundo is the same Allows system to exhibit gradual underow
COMP3320 Lecture 3-19 Copyright c 2008 The Australian National University
3.20
Errors and Floating Point Arithmetic
Require, x op y = f l(x op y)
- IEEE systems achieve this as long as (x op y) is within the range of the FP system
Some familiar laws of real arithmetic are not valid!
- FP addition and multipication are commutative (a b = b a), but NOT associative ((a + b) + c = a + (b + c)) - e.g if is a small positive number slightly less than mach (1 + ) + = 1 + = 1 BUT 1 + ( + ) > 1
COMP3320 Lecture 3-20 Copyright c 2008 The Australian National University
3.21
Cancellation
Subtraction of two t-digit numbers with same sign and similar magnitude yields fewer than t digits - leading digits of two numbers cancel, e.g. 192.403 192.275 000.128 = 1.28000 101 Exact result but cancellation often implies serious loss of information - e.g. <
mach
correct result of 2 entirely lost
(1 + ) (1 ) = 1 1 = 0
COMP3320 Lecture 3-21 Copyright c 2008 The Australian National University
3.22
Error Propagation: Addition/Subtraction
Addition/subtraction
- add or subtract absolute errors b = 2.0 0.1 c = 3.0 0.3 a = b + c = 5.0 0.4 = 0.1/2.0 = 0.05 c = 0.3/3.0 = 0.10 a = 0.4/5.0 = 0.08
b c
relative error averaged
a
=
b b+c b+c
COMP3320 Lecture 3-22 Copyright c 2008 The Australian National University
3.23
Error Propagation: Multiplication/Division
Multiplication/Division b = 2.0 0.1 c = 3.0 0.3 a = = = =
- add or subtract relative errors = 0.1/2.0 = 0.05 c = 0.3/3.0 = 0.10
b
bc (2.0 0.1) (3.0 0.3) 6.0 (0.6 + 0.3) 6.0 0.9
a
= 0.9/6.0 = 0.15
relative error summed
COMP3320 Lecture 3-23 Copyright c 2008 The Australian National University
3.24
Real Life Error Progagation errors = N N = Number of operations = Error per operation
Good algorithms usually give
1GHz UltraSparc performs 2 109 oating point operations per second - Compute for 1 hour N = 7.2 1012 operations - N 3 106 - compare to single precision mach 107
Virtually all accuracy is lost Need to be very carefull about error propagation
COMP3320 Lecture 3-24 Copyright c 2008 The Australian National University
3.25
Assessing Rounding Errors
Examine eect of using higher precision
Consider famous example from Rump in 1988 (IBM S/370) a = 77617.0, b = 33096.0 single precision f = 1.172603 . . . double precision f = 1.1726039400531 . . . extended precision f = 1.172603940053178 . . . However true result is Even the sign is wrong!
f = 333.75b6 + a2(11a2b2 b6 121b4 2) + 5.5b8 + a/(2b)
f = .82739605994682135 5 1017
While this is an extreme example of catastrophic cancellation in an inherently unstable function there must be a better way!
COMP3320 Lecture 3-25 Copyright c 2008 The Australian National University
3.26
Interval Arithmetic: A New Floating Point Paradigm
We are all familiar with representation of uncertainity using the (x ) notation, but this is cumbersome for computation, instead we dene an interval as - a range of numbers bounded by the intervals endpoints - e.g. the range 378 to 432 would be written as [378,432]. Evaluating a function in an interval, gives another interval
- how long to double your money if invested at a compound interest rate of between [3,6]% per annum - answer is the interval [11yr. 7mo., 23yr.2mo.]
Interval computations do not always increase interval width - evaluate f (x) = ln(x)/x over interval X=[2.716,2.718] - answer is the interval [0.3678793,0.3678795]
COMP3320 Lecture 3-26 Copyright c 2008 The Australian National University
3.27
Improved Numerical Accuracy
Two common ways to improve numerical accuracy and reduce interval width
- reduce interval width of input values - increase the number of interval data values used to compute interval bounds
First is obvious the second is more subtle
- given n interval observations of the same true value, the best way to compute a narrow interval bound is not to compute the average, but rather the intersection - every observation must contain the correct value, so must their intersection
COMP3320 Lecture 3-27 Copyright c 2008 The Australian National University
3.28
Support for Interval Arithmetic
C++ and Fortran95 contain language constructs to enable abstract mathematical data types - use of classes/modules, user dened types and operator overloading This is okay for matrices, but less good for intervals
- Suns Fortran95 has been augmented to include an interval data type - there is a proposal to include intervals in future Fortran standards - UltraSPARC III includes some hardware support for interval instructions
COMP3320 Lecture 3-28 Copyright c 2008 The Australian National University
3.29
Use of Intervals
The use of intervals in numerical computation is still very much an area of research, although there is a growing argument that projects like ASCI should require intervals (see Gustafson). Use to track rounding errors - assign to each variable an interval related to machine precision - propagate computation
Use to accurately reect uncertainty Fundamentally new approaches
- used in rening the value of the gravitational constant G
- solution of non-linear equations - global optimisation
We may briey explore the use of intervals in the labs
COMP3320 Lecture 3-29 Copyright c 2008 The Australian National University