11CS 594 Spring 2003Lecture 5:Computer Arithmetic and CacheJack Dongarra Computer Science DepartmentUniversity of Tennessee 2Plan for Todayu Talk about computer arithmetic u Begin Cache discussion3A Little Historyu Von Neumann and Goldstine - 1947 Ø “Can’t expect to solve most big [n>15] linear systems without ca rrying many decimal digits [d>8], otherwise the computed answer would be com pletely inaccurate.” - WRONG!u Turing - 1949Ø “Carrying d digits is equivalent to changing the input data in the d-th place and then solving Ax=b. So if A is only known to d digits, the answer is as accurate as the data deserves.”Ø Backward Error Analysis u Rediscovered in 1961 by Wilkinson and publicizedu Starting in the 1960s- many papers doing backward error analysis of various algorithmsu Many years where each machine did FP arithmetic slightly differe ntlyØ Both rounding and exception handling differedØ Hard to write portable and reliable softwareØ Motivated search for industry-wide standard, beginning late 1970sØ First implementation: Intel 8087u ACM Turing Award 1989 to W. Kahan for design of the IEEE Floating Point Standards 754 (binary) and 854 (decimal)Ø Nearly universally implemented in general purpose machines4Defining Floating Point Arithmeticu Representable numbersØ Scientific notation: +/ - d.d…d x rexpØ sign bit +/-Ø radix r (usually 2 or 10, sometimes 16)Ø significand d.d…d (how many base-r digits d?)Ø exponent exp (range?)Ø others?u Operations:Ø arithmetic: +,-,x,/,... » how to round result to fit in formatØ comparison (<, =, >)Ø conversion between different formats » short to long FP numbers, FP to integerØ exception handling » what to do for 0/0, 2*largest_number, etc.Ø binary/decimal conversion» for I/O, when radix not 10u Language/library support for these operations5IEEE Floating Point Arithmetic Standard 754 - Normalized Numbersu Normalized Nonzero Representable Numbers: +- 1.d…d x 2expØ Macheps = Machine epsilon = 2-#significand bits= relative error in each operationØ OV = overflow threshold = largest numberØ UN = underflow threshold = smallest numberFormat # bits #significand bits macheps #exponent bits exponent range---------- -------- ----------------------- ------------ -------------------- ----------------------Single 32 23+1 2-24(~10-7) 8 2-126- 2127(~10+-38)Double 64 52+1 2-53(~10-16) 11 2-1022- 21023(~10+-308)Double >=80 >=64 <=2-64(~10-19) >=15 2-16382- 216383(~10+-4932)Extended (80 bits on all Intel machines)6IEEE Floating Point Arithmetic Standard 754 -“Denorms ”u Denormalized Numbers: + -0.d…d x 2min_expØ sign bit, nonzero significand , minimum exponentØ Fills in gap between UN and 0u Underflow ExceptionØ occurs when exact nonzero result is less than underflow threshol d UNØ Ex: UN/3Ø return a denorm, or zerou Why bother? Ø Necessary so that following code never divides by zero if (a != b) then x = a/(a-b)27IEEE Floating Point Arithmetic Standard 754 - +- Infinityu+- Infinity: Sign bit, zero significand, maximum exponentuOverflow ExceptionØ occurs when exact finite result too large to represent accuratelyØ Ex: 2*OVØ return +- infinityuDivide by zero ExceptionØ return +- infinity = 1/+-0 Ø sign of zero important!uAlso return +- infinity forØ 3+infinity, 2*infinity, infinity*infinityØ Result is exact, not an exception!8IEEE Floating Point Arithmetic Standard 754 - NAN (Not A Number)uNAN: Sign bit, nonzero significand, maximum exponentuInvalid ExceptionØ occurs when exact result not a well-defined real numberØ 0/0Ø sqrt(-1)Ø infinity-infinity, infinity/infinity, 0*infinityØ NAN + 3Ø NAN > 3?Ø Return a NAN in all these casesuTwo kinds of NANsØ Quiet - propagates without raising an exceptionØ Signaling - generate an exception when touched» good for detecting uninitializeddata9Error AnalysisuBasic error formulaØ fl(a op b) = (a op b)*(1 + d) where» op one of +,-,*,/» | d| <= macheps» assuming no overflow, underflow, or divide by zerouExample: adding 4 numbersØ fl(x1+x2+x3+x4) = {[(x1+x2)*(1+d1) + x3]*(1+d2) + x4}*(1+d3)= x1*(1+d1)*(1+d2)*(1+d3) + x2*(1+d1)*(1+d2)*(1+d3) +x3*(1+d2)*(1+d3) + x4*(1+d3)= x1*(1+e1) + x2*(1+e2) + x3*(1+e3) + x4*(1+e4)where each |ei| <~ 3*machepsØ get exact sum of slightly changed summands xi*(1+ei)Ø Backward Error Analysis - algorithm called numerically stable if it gives the exact result for slightly changed inputsØ Numerical Stability is an algorithm design goal10Backward erroru Approximate solution is exact solution to modified problem.u How large a modification to original problem is required to give result actually obtained?u How much data error in initial input would be required to explain all the error in computed results?u Approximate solution is good if it is exact solution to “nearby” problem.ff(x)f’(x)f’fx’xForward errorBackward error11Sensitivity and Conditioningu Problem is insensitive or well conditioned if relative change in input causes commensurate relative change in solution.u Problem is sensitive or ill-conditioned, if relative change in solution can be much larger than that in input data.Cond = |Relative change in solution| / |Relative change in input data|= |[f(x’) – f(x)]/f(x)| / |(x’ – x)/x|u Problem is sensitive, or ill-conditioned, if cond >> 1.u When function f is evaluated for approximate input x’ = x+h instead of true input value of x.u Absolute error = f(x + h) – f(x) » h f’(x)u Relative error =[ f(x + h) – f(x) ] / f(x) » h f’(x) / f(x)12Sensitivity: 2 Examplescos(p/2) and 2-d System of Equationsu Consider problem of computing cosine function for arguments near p/2.u Let x » p/2 and let h be small perturbation to x. Then absolute error = cos(x+h) – cos(x)» -h sin(x) » -h,relative error » -h tan(x) » 8u So small change in x near p /2 causes large relative change in cos(x) regardless of method used.u cos(1.57079) = 0.63267949 x 10-5 u cos(1.57078) = 1.64267949 x 10-5u Relative change in output is a quarter million times greater than relative change in input..313Sensitivity: 2 Examplescos(p/2) and 2-d System of Equationsu Consider problem of computing cosine function for arguments near p/2.u Let x » p/2 and let h be small perturbation to x. Then absolute error = cos(x+h) – cos(x)» -h sin(x) »
View Full Document
Unlocking...