Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 46
Текст из файла (страница 46)
The subject has something of the sameentertaining, if impractical, flavor as that of fast matrix multiplication, discussedin §2.11.Turn now to algebraic manipulations. You multiply a polynomial of degreen − 1 (array of range [0..n-1]) by a monomial factor x − a by a bit of codelike the following,c[n]=c[n-1];for (j=n-1;j>=1;j--) c[j]=c[j-1]-c[j]*a;c[0] *= (-a);Likewise, you divide a polynomial of degree n by a monomial factor x − a(synthetic division again) usingrem=c[n];c[n]=0.0;for(i=n-1;i>=0;i--) {swap=c[i];c[i]=rem;rem=swap+rem*a;}which leaves you with a new polynomial array and a numerical remainder rem.Multiplication of two general polynomials involves straightforward summingof the products, each involving one coefficient from each polynomial.
Divisionof two general polynomials, while it can be done awkwardly in the fashion taughtusing pencil and paper, is susceptible to a good deal of streamlining. Witness thefollowing routine based on the algorithm in [3].void poldiv(float u[], int n, float v[], int nv, float q[], float r[])Given the n+1 coefficients of a polynomial of degree n in u[0..n], and the nv+1 coefficientsof another polynomial of degree nv in v[0..nv], divide the polynomial u by the polynomialv (“u”/“v”) giving a quotient polynomial whose coefficients are returned in q[0..n], and aremainder polynomial whose coefficients are returned in r[0..n]. The elements r[nv..n]and q[n-nv+1..n] are returned as zero.{int k,j;for (j=0;j<=n;j++) {r[j]=u[j];q[j]=0.0;}for (k=n-nv;k>=0;k--) {q[k]=r[nv+k]/v[nv];for (j=nv+k-1;j>=k;j--) r[j] -= q[k]*v[j-k];}for (j=nv;j<=n;j++) r[j]=0.0;}176Chapter 5.Evaluation of FunctionsRational FunctionsYou evaluate a rational function likeR(x) =p 0 + p1 x + · · · + pµ x µPµ (x)=Qν (x)q 0 + q 1 x + · · · + q ν xν(5.3.4)in the obvious way, namely as two separate polynomials followed by a divide.
As amatter of convention one usually chooses q0 = 1, obtained by dividing numeratorand denominator by any other q0 . It is often convenient to have both sets ofcoefficients stored in a single array, and to have a standard function available fordoing the evaluation:double ratval(double x, double cof[], int mm, int kk)Given mm, kk, and cof[0..mm+kk], evaluate and return the rational function (cof[0] +cof[1]x + · · · + cof[mm]xmm )/(1 + cof[mm+1]x + · · · + cof[mm+kk]xkk ).{int j;double sumd,sumn;Note precision! Change to float if desired.for (sumn=cof[mm],j=mm-1;j>=0;j--) sumn=sumn*x+cof[j];for (sumd=0.0,j=mm+kk;j>=mm+1;j--) sumd=(sumd+cof[j])*x;return sumn/(1.0+sumd);}CITED REFERENCES AND FURTHER READING:Acton, F.S.
1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), pp. 183, 190. [1]Mathews, J., and Walker, R.L. 1970, Mathematical Methods of Physics, 2nd ed. (Reading, MA:W.A. Benjamin/Addison-Wesley), pp. 361–363.
[2]Knuth, D.E. 1981, Seminumerical Algorithms, 2nd ed., vol. 2 of The Art of Computer Programming(Reading, MA: Addison-Wesley), §4.6. [3]Fike, C.T. 1968, Computer Evaluation of Mathematical Functions (Englewood Cliffs, NJ: PrenticeHall), Chapter 4.Winograd, S. 1970, Communications on Pure and Applied Mathematics, vol. 23, pp.
165–179. [4]Kronsjö, L. 1987, Algorithms: Their Complexity and Efficiency, 2nd ed. (New York: Wiley). [5]5.4 Complex ArithmeticAs we mentioned in §1.2, the lack of built-in complex arithmetic in C is anuisance for numerical work. Even in languages like FORTRAN that have complexdata types, it is disconcertingly common to encounter complex operations thatproduce overflows or underflows when both the complex operands and the complexresult are perfectly representable.
This occurs, we think, because software companiesassign inexperienced programmers to what they believe to be the perfectly trivialtask of implementing complex arithmetic.1775.4 Complex ArithmeticActually, complex arithmetic is not quite trivial. Addition and subtractionare done in the obvious way, performing the operation separately on the real andimaginary parts of the operands. Multiplication can also be done in the obvious way,with 4 multiplications, one addition, and one subtraction,(a + ib)(c + id) = (ac − bd) + i(bc + ad)(5.4.1)(the addition before the i doesn’t count; it just separates the real and imaginary partsnotationally).
But it is sometimes faster to multiply via(a + ib)(c + id) = (ac − bd) + i[(a + b)(c + d) − ac − bd](5.4.2)which has only three multiplications (ac, bd, (a + b)(c + d)), plus two additions andthree subtractions. The total operations count is higher by two, but multiplicationis a slow operation on some machines.While it is true that intermediate results in equations (5.4.1) and (5.4.2) canoverflow even when the final result is representable, this happens only when the finalanswer is on the edge of representability. Not so for the complex modulus, if youare misguided enough to compute it as|a + ib| =(a2 + b 2(bad!)(5.4.3)whose intermediate result will overflow if either a or b is as large as the squareroot of the largest representable number (e.g., 1019 as compared to 1038 ).
The rightway to do the calculation is+|a + ib| =(|a|( 1 + (b/a)2|b| 1 + (a/b)2|a| ≥ |b||a| < |b|(5.4.4)Complex division should use a similar trick to prevent avoidable overflows,underflow, or loss of precision,[a + b(d/c)] + i[b − a(d/c)]a + ib c + d(d/c)= [a(c/d) + b] + i[b(c/d) − a]c + id c(c/d) + d|c| ≥ |d|(5.4.5)|c| < |d|Of course you should calculate repeated subexpressions, like c/d or d/c, only once.Complex square root is even more complicated, since we must both guardintermediate results, and also enforce a chosen branch cut (here taken to be thenegative real axis). To take the square root of c + id, first compute0 *(2( |c| 1 + 1 + (d/c)2w≡*((|c/d| + 1 + (c/d)2 |d|2c=d=0|c| ≥ |d||c| < |d|(5.4.6)178Chapter 5.Evaluation of FunctionsThen the answer is√0 w=0dw= 0, c ≥ 0w+i2wc + id =|d|+ iw2w |d|− iw2ww = 0, c < 0, d ≥ 0(5.4.7)w = 0, c < 0, d < 0Routines implementing these algorithms are listed in Appendix C.CITED REFERENCES AND FURTHER READING:Midy, P., and Yakovlev, Y.
1991, Mathematics and Computers in Simulation, vol. 33, pp. 33–49.Knuth, D.E. 1981, Seminumerical Algorithms, 2nd ed., vol. 2 of The Art of Computer Programming(Reading, MA: Addison-Wesley) [see solutions to exercises 4.2.1.16 and 4.6.4.41].5.5 Recurrence Relations and Clenshaw’sRecurrence FormulaMany useful functions satisfy recurrence relations, e.g.,(n + 1)Pn+1 (x) = (2n + 1)xPn (x) − nPn−1 (x)Jn+1 (x) =2nJn (x) − Jn−1 (x)x(5.5.1)(5.5.2)nEn+1 (x) = e−x − xEn (x)(5.5.3)cos nθ = 2 cos θ cos(n − 1)θ − cos(n − 2)θ(5.5.4)sin nθ = 2 cos θ sin(n − 1)θ − sin(n − 2)θ(5.5.5)where the first three functions are Legendre polynomials, Bessel functions of the firstkind, and exponential integrals, respectively. (For notation see [1].) These relationsare useful for extending computational methods from two successive values of n toother values, either larger or smaller.Equations (5.5.4) and (5.5.5) motivate us to say a few words about trigonometricfunctions.
If your program’s running time is dominated by evaluating trigonometricfunctions, you are probably doing something wrong. Trig functions whose argumentsform a linear sequence θ = θ0 + nδ, n = 0, 1, 2, . . ., are efficiently calculated bythe following recurrence,cos(θ + δ) = cos θ − [α cos θ + β sin θ]sin(θ + δ) = sin θ − [α sin θ − β cos θ](5.5.6)5.5 Recurrence Relations and Clenshaw’s Recurrence Formulawhere α and β are the precomputed coefficients δβ ≡ sin δα ≡ 2 sin22179(5.5.7)The reason for doing things this way, rather than with the standard (and equivalent)identities for sums of angles, is that here α and β do not lose significance if theincremental δ is small.
Likewise, the adds in equation (5.5.6) should be done inthe order indicated by square brackets. We will use (5.5.6) repeatedly in Chapter12, when we deal with Fourier transforms.Another trick, occasionally useful, is to note that both sin θ and cos θ can becalculated via a single call to tan: θt ≡ tan2cos θ =1 − t21 + t2sin θ =2t1 + t2(5.5.8)The cost of getting both sin and cos, if you need them, is thus the cost of tan plus2 multiplies, 2 divides, and 2 adds. On machines with slow trig functions, this canbe a savings. However, note that special treatment is required if θ → ±π.
And alsonote that many modern machines have very fast trig functions; so you should notassume that equation (5.5.8) is faster without testing.Stability of RecurrencesYou need to be aware that recurrence relations are not necessarily stableagainst roundoff error in the direction that you propose to go (either increasing n ordecreasing n). A three-term linear recurrence relationyn+1 + an yn + bn yn−1 = 0,n = 1, 2, .
. .(5.5.9)has two linearly independent solutions, fn and gn say. Only one of these correspondsto the sequence of functions fn that you are trying to generate. The other one gnmay be exponentially growing in the direction that you want to go, or exponentiallydamped, or exponentially neutral (growing or dying as some power law, for example).If it is exponentially growing, then the recurrence relation is of little or no practicaluse in that direction.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
