Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C, страница 12

That isanother reason for our being rigorous about using the ANSI prototype mechanism,and a good reason for you to use an ANSI-compatible compiler.Some older C compilers do provide an optional compilation mode in whichthe implicit conversion of float to double is suppressed. Use this if you can.In this book, when we write float, we mean float; when we write double,we mean double, i.e., there is a good algorithmic reason for having higherprecision. Our routines all can tolerate the traditional implicit conversion rules,but they are more efficient without them.

Of course, if your application actuallyrequires double precision, you can change our declarations from float to doublewithout difficulty. (The brute force approach is to add a preprocessor statement#define float double !)A Few WrinklesWe like to keep code compact, avoiding unnecessary spaces unless they addimmediate clarity. We usually don’t put space around the assignment operator “=”.Through a quirk of history, however, some C compilers recognize the (nonexistent)operator “=-” as being equivalent to the subtractive assignment operator “-=”, and“=*” as being the same as the multiplicative assignment operator “*=”. That is whyyou will see us write y= -10.0; or y=(-10.0);, and y= *a; or y=(*a);.We have the same viewpoint regarding unnecessary parentheses.

You can’t write(or read) C effectively unless you memorize its operator precedence and associativityrules. Please study the accompanying table while you brush your teeth every night.We never use the register storage class specifier. Good optimizing compilersare quite sophisticated in making their own decisions about what to keep in registers,and the best choices are sometimes rather counter-intuitive.Different compilers use different methods of distinguishing between definingand referencing declarations of the same external name in several files. We followthe most common scheme, which is also the ANSI standard.

The storage classextern is explicitly included on all referencing top-level declarations. The storageclass is omitted from the single defining declaration for each external variable. Wehave commented these declarations, so that if your compiler uses a different schemeyou can change the code. The various schemes are discussed in §4.8 of [1].26Chapter 1.PreliminariesOperator Precedence and Associativity Rules in Cfunction callarray elementstructure or union memberpointer reference to structureleft-to-rightlogical notbitwise complementunary minusincrementdecrementaddress ofcontents ofcast to typesize in bytesright-to-left*/%multiplydivideremainderleft-to-right+-addsubtractleft-to-right<<>>bitwise left shiftbitwise right shiftleft-to-right<><=>=arithmetic less thanarithmetic greater thanarithmetic less than or equal toarithmetic greater than or equal toleft-to-right==!=arithmetic equalarithmetic not equalleft-to-right&bitwise andleft-to-right^bitwise exclusive orleft-to-right|bitwise orleft-to-right&&logical andleft-to-right||logical orleft-to-rightconditional expressionright-to-left()[].->!~++-&*(type)sizeof?:=assignment operatoralso += -= *= /= %=<<= >>= &= ^= |=,sequential expressionright-to-leftleft-to-rightWe have already alluded to the problem of computing small integer powers ofnumbers, most notably the square and cube.

The omission of this operation from Cis perhaps the language’s most galling insult to the scientific programmer. All goodFORTRAN compilers recognize expressions like (A+B)**4 and produce in-line code,in this case with only one add and two multiplies. It is typical for constant integerpowers up to 12 to be thus recognized.1.2 Some C Conventions for Scientific Computing27In C, the mere problem of squaring is hard enough! Some people “macro-ize”the operation as#define SQR(a) ((a)*(a))However, this is likely to produce code where SQR(sin(x)) results in two calls tothe sine routine! You might be tempted to avoid this by storing the argument of thesquaring function in a temporary variable:static float sqrarg;#define SQR(a) (sqrarg=(a),sqrarg*sqrarg)The global variable sqrarg now has (and needs to keep) scope over the wholemodule, which is a little dangerous.

Also, one needs a completely different macro tosquare expressions of type int. More seriously, this macro can fail if there are twoSQR operations in a single expression. Since in C the order of evaluation of pieces ofthe expression is at the compiler’s discretion, the value of sqrarg in one evaluationof SQR can be that from the other evaluation in the same expression, producingnonsensical results. When we need a guaranteed-correct SQR macro, we use thefollowing, which exploits the guaranteed complete evaluation of subexpressions ina conditional expression:static float sqrarg;#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)A collection of macros for other simple operations is included in the file nrutil.h(see Appendix B) and used by many of our programs. Here are the synopses:SQR(a)DSQR(a)FMAX(a,b)FMIN(a,b)DMAX(a,b)DMIN(a,b)IMAX(a,b)IMIN(a,b)LMAX(a,b)LMIN(a,b)SIGN(a,b)Square a float value.Square a double value.Maximum of two float values.Minimum of two float values.Maximum of two double values.Minimum of two double values.Maximum of two int values.Minimum of two int values.Maximum of two long values.Minimum of two long values.Magnitude of a times sign of b.Scientific programming in C may someday become a bed of roses; for now,watch out for the thorns!CITED REFERENCES AND FURTHER READING:Harbison, S.P., and Steele, G.L., Jr.

1991, C: A Reference Manual, 3rd ed. (Englewood Cliffs,NJ: Prentice-Hall). [1]AT&T Bell Laboratories 1985, The C Programmer’s Handbook (Englewood Cliffs, NJ: PrenticeHall).Kernighan, B., and Ritchie, D. 1978, The C Programming Language (Englewood Cliffs, NJ:Prentice-Hall). [Reference for K&R “traditional” C. Later editions of this book conform tothe ANSI C standard.]Hogan, T. 1984, The C Programmer’s Handbook (Bowie, MD: Brady Communications).28Chapter 1.Preliminaries1.3 Error, Accuracy, and StabilityAlthough we assume no prior training of the reader in formal numerical analysis,we will need to presume a common understanding of a few key concepts.

We willdefine these briefly in this section.Computers store numbers not with infinite precision but rather in some approximation that can be packed into a fixed number of bits (binary digits) or bytes (groupsof 8 bits). Almost all computers allow the programmer a choice among severaldifferent such representations or data types. Data types can differ in the number ofbits utilized (the wordlength), but also in the more fundamental respect of whetherthe stored number is represented in fixed-point (int or long) or floating-point(float or double) format.A number in integer representation is exact. Arithmetic between numbers ininteger representation is also exact, with the provisos that (i) the answer is not outsidethe range of (usually, signed) integers that can be represented, and (ii) that divisionis interpreted as producing an integer result, throwing away any integer remainder.In floating-point representation, a number is represented internally by a sign bits (interpreted as plus or minus), an exact integer exponent e, and an exact positiveinteger mantissa M .

Taken together these represent the numbers × M × B e−E(1.3.1)where B is the base of the representation (usually B = 2, but sometimes B = 16),and E is the bias of the exponent, a fixed integer constant for any given machineand representation. An example is shown in Figure 1.3.1.Several floating-point bit patterns can represent the same number. If B = 2,for example, a mantissa with leading (high-order) zero bits can be left-shifted, i.e.,multiplied by a power of 2, if the exponent is decreased by a compensating amount.Bit patterns that are “as left-shifted as they can be” are termed normalized.

Mostcomputers always produce normalized results, since these don’t waste any bits ofthe mantissa and thus allow a greater accuracy of the representation. Since thehigh-order bit of a properly normalized mantissa (when B = 2) is always one, somecomputers don’t store this bit at all, giving one extra bit of significance.Arithmetic among numbers in floating-point representation is not exact, even ifthe operands happen to be exactly represented (i.e., have exact values in the form ofequation 1.3.1). For example, two floating numbers are added by first right-shifting(dividing by two) the mantissa of the smaller (in magnitude) one, simultaneouslyincreasing its exponent, until the two operands have the same exponent.

Low-order(least significant) bits of the smaller operand are lost by this shifting. If the twooperands differ too greatly in magnitude, then the smaller operand is effectivelyreplaced by zero, since it is right-shifted to oblivion.The smallest (in magnitude) floating-point number which, when added to thefloating-point number 1.0, produces a floating-point result different from 1.0 istermed the machine accuracy m . A typical computer with B = 2 and a 32-bitwordlength has m around 3 × 10−8 . (A more detailed discussion of machinecharacteristics, and a program to determine them, is given in §20.1.) Roughly29aantiss-bitm23thbe is b“p it cha ount ldom”8bitexpsignbitonent1.3 Error, Accuracy, and Stability= 01000000010000000000000000000000(a)3 = 01000001011000000000000000000000(b)1⁄ 2= 00111111110000000000000000000000(c)10 − 7 = 001101001(d)= 0100000101 1 0 ...10110101111111001010000000000000000000000003 + 10 − 7 = 01000001011000000000000000000000(f )1⁄ 4(e)Figure 1.3.1.

Floating point representations of numbers in a typical 32-bit (4-byte) format. (a) Thenumber 1/2 (note the bias in the exponent); (b) the number 3; (c) the number 1/4; (d) the number10−7 , represented to machine accuracy; (e) the same number 10−7 , but shifted so as to have the sameexponent as the number 3; with this shifting, all significance is lost and 10−7 becomes zero; shifting toa common exponent must occur before two numbers can be added; (f) sum of the numbers 3 + 10−7,which equals 3 to machine accuracy.

Even though 10−7 can be represented accurately by itself, it cannotaccurately be added to a much larger number.speaking, the machine accuracy m is the fractional accuracy to which floating-pointnumbers are represented, corresponding to a change of one in the least significantbit of the mantissa. Pretty much any arithmetic operation among floating numbersshould be thought of as introducing an additional fractional error of at least m .

Thistype of error is called roundoff error.It is important to understand that m is not the smallest floating-point numberthat can be represented on a machine. That number depends on how many bits thereare in the exponent, while m depends on how many bits there are in the mantissa.Roundoff errors accumulate with increasing amounts of calculation. If, in thecourse of obtaining a calculated value, you perform N such arithmetic operations,√you might be so lucky as to have a total roundoff error on the order of N m , ifthe roundoff errors come in randomly up or down.