Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 92
Текст из файла (страница 92)
An interface betweenBrent’s package and the Fortran precompiled Augment [254, 1979] is describedby Brent, Hooper, and Yohe [147, 1980]; it makes Brent’s package much easierto use.25.9 M ULTIPLE P RECISION A RITHMETIC505Recently, Bailey has written a multiple precision arithmetic package MPFUN [44, 1993], which consists of about 10,000 lines of Fortran 77 code in 87subprograms.
In MPFUN, a multiprecision (MP) number is a vector of singleprecision floating point numbers; it represents a number in base 224 (for IEEEarithmetic). Complex multiprecision numbers are also supported.MPFUN routines are available to perform the basic arithmetic operations,evaluate nth roots and transcendental functions, compare MP numbers, produce a random MP number, solve polynomial equations, and perform otheroperations. For many of these routines both simple and advanced algorithmsare provided; the advanced algorithms are intended for computations at precision levels above 1000 digits.
One advanced algorithm is a fast Fouriertransform technique for evaluating the convolution involved in a multiplication. Another is used for division: x/y is evaluated as x * (1/y ), where 1/y isevaluated by Newton’s method (see Problem 2.26), which involves only multiplications. An interesting aspect of Newton’s method in this context is thatit can be implemented with a precision level that doubles at each stage, sincethe iteration damps out any errors. Another interesting feature of a variableprecision environment is that numerical instability can be tackled simply byincreasing the precision. As an example, MPFUN does complex multiplicationusing the 3M methodwhich uses three real multiplications instead of the usual four.
As we saw in§22.2.4, the 3M method produces a computed imaginary part that can havelarge relative error, but this instability is not a serious obstacle to its use inMPFUN.Bailey provides a translator that takes Fortran 77 source containing directives in comments that specify the precision level and which variables are tobe treated as multiprecision, and produces a Fortran 77 program containingthe appropriate multiprecision subroutine calls. He also provides a Fortran 90version of the package that employs derived data types and operator extensions [45, 1994]. This Fortran 90 version is a powerful and easy to use toolfor doing high precision numerical computations.Bailey’s own use of the packages includes high precision computation ofconstants such as π [42, 1988] and Euler’s constant γ. A practical applicationto a fluid flow problem is described by Bailey, Krasny, and Pelz [46, 1993].They found that converting an existing code to use MPFUN routines increasedexecution times by a factor of about 400 at 56 decimal digit precision, andthey comment that this magnitude of increase is fairly typical for multipleprecision computation.Another recent Fortran 77 program for multiple precision arithmetic is theFM package of Smith [926, 1991].
It is functionally quite similar to Bailey’sMPFUN. No translator is supplied for converting from standard Fortran 77506S OFTWARE I SSUESINFLOATING P OINT A RITHMETICsource code to code that invokes FM, but Smith notes that the general purposeprecompiled Augment [254, 1979] can be used for this purpose.Fortran routines for high precision computation can also be found in Numerical Recipes [842, 1992, §20.6], and high precision numerical computationis supported by many symbolic manipulation packages, including Maple [199,1991] and Mathematical [1109, 1991].A Pascal-like programming language called Turing [579, 1988] developedat the University of Toronto in the 1980s is the basis of an extension calledNumerical Turing, developed by Hull and his co-workers [591, 1985]. Numerical Turing includes decimal floating point arithmetic whose precision can bedynamically varied throughout the course of a program, a feature argued forin [590, 1982] and implemented in hardware in [231, 1983].An extension to the level-2 BLAS (see Appendix D) is proposed in [313,1988, App.
B], comprising routines having the same specifications as those inthe standard BLAS but which calculate in extended precision.25.10. Patriot Missile Software ProblemA report from the United States General Accounting Office begins “On February 25, 1991, a Patriot missile defense system operating at Dhahran, SaudiArabia, during Operation Desert Storm failed to track and intercept an incoming Scud.
This Scud subsequently hit an Army barracks, killing 28 Americans”[1036, 1992]. The report finds that the failure to track the Scud missile wascaused by a precision problem in the software.The computer used to control the Patriot missile is baaed on a 1970s designand uses 24-bit arithmetic. The Patriot system tracks its targets by measuringthe time it takes for radar pulses to bounce back from them. Time is recordedby the system clock in tenths of a second, but is stored as an integer. To enabletracking calculations the time is converted to a 24 bit floating point number.Rounding errors in the time conversions cause shifts in the system’s “rangegate”, which is used to track the target.On February 11, 1991 the Patriot Project Office received field data identifying a 20% shift in the Patriot system’s range gate after the system had beenrunning continuously for 8 hours.
This data implied that after 20 consecutivehours of use the system would fail to track and intercept a Scud. Modifiedsoftware that compensated for the inaccurate time calculations was releasedon February 16 by army officials. On February 25, Alpha Battery, which wasprotecting the Dhahran Air Base, had been in continuous operation for over100 hours. The inaccurate time calculations caused the range gate to shift somuch that the system could not track the incoming Scud. On February 26,the next day, the modified software arrived in Dhahran.
Table 25.2, takenfrom [1036, 1992], shows clearly how, with increasing time of operation, the25.11 N OTESANDR EFERENCES507Table 25.2. Effect of extended run time on Patriot missile operation.Hours Seconds Calculated time(seconds)000136003599.99662880028799.9725820a7200071999.931348172800172799.8352259200259199.752872100 b360000359999.6567 cabcInaccuracy Approximate shift inrange gate (meters)(seconds)007.003455.0275137.0687.1648330494.2472687.3433For continuous operation exceeding 20 hours target is outside range gate.Alpha battery ran continuously for about 100 hours.Corrected value.Patriot lost track of its target.
Note that the numbers in Table 25.2 are consistent with a relative error of 2–20 in the computer’s representation of 0.1,this constant being used to convert from the system clock’s tenths of a secondto seconds (2–20 is the relative error introduced by chopping 0.1 to 23 bitsafter the binary point).25.11. Notes and ReferencesThe issues discussed in this chapter are rarely treated in textbooks. One bookthat is worth consulting is Miller’s Engineering of Numerical Software [759,1984], which is concerned with the designing of reliable numerical softwareand presents practical examples.Kahan’s paper [627, 1972] is not, as the title suggests, a survey of erroranalysis, but treats in detail some specific problems, including solution of aquadratic equation and summation, giving error analysis, programs illustrating foibles of Fortran compilers, and ingenious workarounds.The collection of papers improving Floating-Point Programming [1061,1990] contains, in addition to an introduction to floating point arithmetic andBrown’s model, descriptions of recent work on embedding interval arithmeticwith a super-accurate inner product in Pascal and ADA (see §24.4).There are several algorithms for evaluating a continued fraction; see, forexample, Blanch [123, 1964] and Press, Teukolsky, Vetterling, and Flannery[842, 1992, §5.2].
Rounding error analysis for the “bottom up” algorithm isgiven by Jones and Thron [617, 1974].A standard set of names and definitions of machine parameters was proposed by the IFIP Working Group on Numerical Software (WG 2.5) [383,S OFTWARE I SSUES508INFLOATING P OINT A RITHMETIC1978], though these do not appear to have been widely adopted.The use of high precision computations to investigate the Riemann hypothesis and other mathematical problems and conjectures is described byVarga [1052, 1990], who has made use of Brent’s package [145, 1978], [146,1978].Problems25.1.
(a) The following MATLAB code attempts to computeµ := min{floating point x : fl(1 + x) > 1},under the assumption that the base β = 2.function x = mu()x = 1;xp1 = x + l;while xp1 > 1x = x/2 ;xp1 = x + 1;endx = 2*x;On my workstation, running this code gives>> muans =2.2204e-016Try this code, or a translation of it into another language, on any machinesavailable to you. Are the answers what you expected?(b) Under what circumstances might the code in (a) fail to compute µ?(Hint: consider optimizing compilers.)(c) What is the value of µ in terms of the arithmetic parameters β andt ? Note: the answer depends subtly on how rounding is performed, and inparticular on whether double rounding is in effect; see Higham [550, 1991]and Moler [771, 1990].
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
