Nash - Compact Numerical Methods for Computers (523163), страница 3
Текст из файла (страница 3)
Computing the eigensolutions of a generalsquare matrix is a problem with many inherent difficulties, but we shall not dwell onthese at length.Constrained optimisation is still a topic where I would be happier to offer morematerial, but general methods of sufficient simplicity to include in a handbook of thissort have eluded my research efforts. In particular, the mathematical programmingproblem is not treated here.Since the aim has been to give a problem-solving person some tools with whichto work, the mathematical detail in the pages that follow has been mostly confinedto that required for explanatory purposes. No claim is made to rigour in any‘proof’, though a sincere effort has been made to ensure that the statement oftheorems is correct and precise.1.2. MACHINE CHARACTERISTICSIn the first edition, a ‘small computer’ was taken to have about 6000 characters ofmain memory to hold both programs and data.
This logical machine, which might bea part of a larger physical computer, reflected the reality facing a quite large group ofusers in the mid- to late-1970s.A more up-to-date definition of ‘small computer’ could be stated, but it is notreally necessary. Users of this book are likely to be those scientists, engineers, andstatisticians who must, for reasons of financial or administrative necessity orconvenience, carry out their computations in environments where the programs4Compact numerical methods for computerscannot be acquired simply and must, therefore, be written in-house. There are also anumber of portable computers now available. This text is being entered on a TandyRadio Shack TRS-80 Model 100, which is only the size of a large book and ispowered by four penlight batteries.Users of the various types of machines mentioned above often do not have muchchoice as to the programming tools available.
On ‘borrowed’ computers, one has toput up with the compilers and interpreters provided by the user who has paid for theresource. On portables, the choices may be limited by the decisions of the manufacturer. In practice, I have, until recently, mostly programmed in BASIC, despite itslimitations, since it has at least been workable on most of the machines available tome.Another group of users of the material in this book is likely to be softwaredevelopers. Some scenarios which have occurred are:—software is being developed in a particular computing environment (e.g. LISPfor artificial intelligence) and a calculation is required for which suitable off-the-shelfroutines are not available;—standard routines exist but when linked into the package cause the executablecode to be too large for the intended disk storage or memory;—standard routines exist, but the coding is incompatible with the compiler orinterpreter at hand.It is worth mentioning that programming language standards have undergoneconsiderable development in the last decade.
Nevertheless, the goal of portablesource codes of numerical methods is still only won by careful and conservativeprogramming practices.Because of the emphasis on the needs of the user to program the methods, there isconsiderable concern in this book to keep the length and complexity of thealgorithms to a minimum. Ideally, I like program codes for the algorithms to be nolonger than a page of typed material, and at worse, less than three pages. This makesit feasible to implement the algorithms in a single work session. However, even thislevel of effort is likely to be considered tedious and it is unnecessary if the code can beprovided in a suitable form. Here we provide source code in Turbo Pascal for thealgorithms in this book and for the driver and support routines to run them (underTurbo Pascal version 3.01a).The philosophy behind this book remains one of working with available toolsrather than complaining that better ones exist, albeit not easily accessible.
Thisshould not lead to complacency in dealing with the machine but rather to an activewariness of any and every feature of the system. A number of these can and should bechecked by using programming devices which force the system to reveal itself in spiteof the declarations in the manual(s). Others will have to be determined by exploringevery error possibility when a program fails to produce expected results. In mostcases programmer error is to blame, but I have encountered at least one system errorin each of the systems I have used seriously.
For instance, trigonometric functions areusually computed by power series approximation. However, these approximationshave validity over specified domains, usually [0, π /2] or [0, π /2] (see Abramowitz andStegun 1965, p 76). Thus the argument of the function must first be transformed toA starting point5bring it into the appropriate range.
For exampleorsin( π – φ ) = sin(1.1)sin( π /2 – φ ) = cos(1.2)Unless this range reduction is done very carefully the results may be quiteunexpected. On one system, hosted by a Data General NOVA, I have observedthat the sine of an angle near π and the cosine of an angle near π/2 were bothcomputed as unity instead of a small value, due to this type of error. Similarly, onsome early models of Hewlett- Packard pocket calculators, the rectangular-to-polarcoordinate transformation may give a vector 180° from the correct direction.
(Thisappears to have been corrected now.)Testing the quality of the floating-point arithmetic and special functions istechnically difficult and tedious. However, some developments which aid the userhave been made by public-spirited scientists. Of these, I consider the most worthyexample to be PARANOIA, a program to examine the floating-point arithmeticprovided by a programming language translator.
Devised originally by Professor WKahan of the University of California, Berkeley, it has been developed and distributed in a number of programming languages (Karpinski 1985). Its output isdidactic, so that one does not have to be a numerical analyst to interpret the results. Ihave used the BASIC, FORTRAN, Pascal and C versions of PARANOIA, and have seenreports of Modula-2 and ADA®† versions.In the area of special functions, Cody and Waite (1980) have devised software toboth calculate and test a number of the commonly used transcendental functions(sin, cos, tan, log, exp, sqrt, xy).
The ELEFUNT testing software is available in theirbook, written in FORTRAN. A considerable effort would be needed to translate it intoother programming languages.An example from our own work is the program DUBLTEST, which is designed todetermine the precision to which the standard special functions in BASIC arecomputed (Nash and Walker-Smith 1987). On the IBM PC, many versions ofMicrosoft BASIC (GWBASIC, BASICA) would only compute such functions in singleprecision, even if the variables involved as arguments or results were doubleprecision. For some nonlinear parameter estimation problems, the resulting lowprecision results caused unsatisfactory operation of our software.Since most algorithms are in some sense iterative, it is necessary that one hassome criterion for deciding when sufficient progress has been made that theexecution of a program can be halted. While, in general, I avoid tests whichrequire knowledge of the machine, preferring to use the criterion that no progresshas been made in an iteration, it is sometimes convenient or even necessary toemploy tests involving tolerances related to the structure of the computing deviceat hand.The most useful property of a system which can be determined systematically isthe machine precision.
This is the smallest number, eps, such that1+eps>1† ADA is a registered name of the US Department of Defense.(1.3)6Compact numerical methods for computerswithin the arithmetic of the system. Two programs in FORTRAN for determining themachine precision, the radix or base of the arithmetic, and machine rounding ortruncating properties have been given by Malcolm (1972). The reader is cautionedthat, since these programs make use of tests of conditions like (1.3), they may befrustrated by optimising compilers which are able to note that (1.3) in exactarithmetic is equivalent toeps>0.(1.4)Condition (1.4) is not meaningful in the present context.
The Univac compilershave acquired some notoriety in this regard, but they are by no means the onlyoffenders.To find the machine precision and radix by using arithmetic of the computeritself, it is first necessary to find a number q such that (1 + q ) and q arerepresented identically, that is, the representation of 1 having the same exponentas q has a digit in the (t + 1)th radix position where t is the number of radix digitsin the floating-point mantissa. As an example, consider a four decimal digitmachine.
If q = 10,000 or greater, then q is represented as (say)0.1 * 1E5while 1 is represented as0·00001 * 1E5.The action of storing the five-digit sum0·10001 * 1E5in a four-digit word causes the last digit to be dropped. In the example,q = 10 000 is the smallest number which causes (1 + q) and q to be representedidentically, but any numberq > 9999will have the same property. If the machine under consideration has radix R, thenanyq > Rt(1.5)t +1will have the desired property.
If, moreover, q and R are represented so thatthenq < Rt+1q+R>q.(1.6)(1.7)In our example, R = 10 and t = 4 so the largest q consistent with (1.6) isq = 105-10 = 99 990 = 0·9999 * 1E5and99 990 + 10 = 100 000 = 0·1000 * 1E6 > q.Starting with a trial value, say q = 1, successive doubling will give some numberq = 2kA starting point7such that (q + 1) and q are represented identically. By then setting r to successiveintegers 2, 3, 4 , . . . , a value such thatq+r>q(1.8)will be found.
On a machine which truncates, r is then the radix R. However, ifthe machine rounds in some fashion, the condition (1.8) may be satisfied for r < R.Nevertheless, the representations of q and (q + r) will differ by R. In the example,doubling will produce q = 16 384 which will be represented as0·1638 * 1E5so q + r is represented as0·1639 * 1E5for some r 10. Then subtraction of these gives0·0001 * 1E5 = 10.Unfortunately, it is possible to foresee situations where this will not work.Suppose that q = 99 990, then we have0·9999 * 1E5 + 10 = 0·1000 * 1E6and0·1000 * 1E6–0·9999 * 1E5 = R'.But if the second number in this subtraction is first transformed to0·0999 * 1E6then R´ is assigned the value 100. Successive doubling should not, unless themachine arithmetic is extremely unusual, give q this close to the upper bound of(1.6).Suppose that R has been found and that it is greater than two.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
