Arndt - Algorithms for Programmers, страница 52

ARITHMETICAL ALGORITHMS301where ∆ is the difference operator.Note that the mathematically equivalent formx∗k=2xk xk+2 − xk+1xk+2 − 2 xk+1 + xk(13.170)sometimes given should be avoided with numerical computations due to possible cancellation.PkIf xk = i=0 ai and the ratio of consecutive summands ai is approximately constant (that is, a is close toa geometric series) then x∗ converges significantly faster to x∞ than x. Let us partly rewrite the formulausing xk − xk−1 = ak :2x∗k= xk+2 −(ak+2 )ak+2 − ak+1(13.171)Then for a geometric series with ak+1 /ak = qx∗k¡¢22a0 q k+2(ak+2 )= xk+2 −= xk+2 −ak+2 − ak+1a0 (q k+2 − q k+1 )¢1 − q k+3q k+2a0 ¡= a0+ a0 q k+2 · k+1=1 − q k+3 + q k+3k+21−qq−q1−qa0=1−q(13.172a)(13.172b)(13.172c)which is the exact sum.Why do we meet the delta squared here? Consider the sequencex0 ,x1 = Φ(x0 ),x2 = Φ(x1 ) = Φ(Φ(x0 )),...(13.173)of better and better approximations to some root r of a function f .

Think of the xk as partial sums of aseries whose sum is the root r. Apply the idea to define an improved iteration Φ∗ from a given Φ:2Φ∗ (x)=Φ(Φ(x)) −(Φ(Φ(x)) − Φ(x))Φ(Φ(x)) − 2 Φ(x) + x(13.174)The good news is that Φ∗ will give quadratic convergence even√if Φ only gives linear convergence. Asan example let us take f (x) = (x2 − d)2 , forget that its root d is a double root and happily defineΦ(x) = x − f (x)/f 0 (x) = x − (x2 − d)/(4x). Convergence is only linear:¶µ√√e e2+ O(e3 )(13.175)Φ( d · (1 + e)) =d· 1+ +24Then tryΦ∗ (x) =d (7 x2 + d)x (3 x2 + 5 d)(13.176)and find that it offers quadratic convergence√Φ( d · (1 + e))=√¶µe3e24++ O(e )d· 1−416(13.177)In general, if Φ has convergence of order n > 1 then Φ∗ will be of order 2n − 1, linear convergence (n = 1)is turned into second order.

(See [10]).CHAPTER 13. ARITHMETICAL ALGORITHMS13.8.6302Improvements of the delta squared process *Let us rewrite the partial sums of equation 13.168 as values of a function F : F (k) := xk . The deltasquared process extrapolated the function F :F (∞)≈ F (k) −F 0 (k)2F 00 (k)(13.178)if we approximate the n-th derivatives using the backward difference operator ∆ asF 0 (k) ≈F 00 (k) ≈F (3) (k) ≈F (n)≈∆ F (k) = F (k) − F (k − 1)∆2 F (k) = F (k) − 2 F (k − 1) + F (k − 2)∆3 F (k) = F (k) − 3 F (k − 1) + 3 F (k − 2) − F (k − 3)µ ¶kXni k∆ F (k) =(−1)F (k − i)ii=0(13.179a)(13.179b)(13.179c)(13.179d)and identifyF (0)F (1)F (2)F (3)F (k)=====x0f (x) = x1f (f (x)) = x2f (f (f (x))) = x3f ◦k (x) = xk(13.180a)(13.180b)(13.180c)(13.180d)(13.180e)where f ◦k (.) denotes the k-fold composition of f (.).

We had f = Φ before. Write the Taylor expansionasF (k + D) =F (k) + F 0 (k) D +F 00 (k) 2 F 000 (k) 3D +D + ...26(13.181)Use the second order iteration for an extremum of FΦ {F 0 }2 (k) =X−F 0 (k)F 00 (k)(13.182)set X := F (k)Φ {F 0 }2 (k) = F (k) −F 0 (k)= F (k) + DF 00 (k)(13.183)and insert D = Φ − F into equation 13.181 up to first order:F (k + D) =F (k) + F 0 (k) D = F (k) + F 0 (k)F 0 (k)F 00 (k)(13.184)This is Aitken’s delta squared method.Note that for the numerical process it is not necessary to know F algebraically, one just needs the valuesof F (k), the partial sums.There were two degrees of freedom in our derivation: the order a after which the Taylor series in D wastruncated and the order b of the iteration used for D = Φ {F 0 } − F .

Define an extrapolation operatorT[a,b] (F ):=whereF (k) + F 0 (k) D +0D := Φ {F }b − FF (a) aF 00 (k) 2D + ... +D2a!(13.185)(13.186)CHAPTER 13. ARITHMETICAL ALGORITHMS303Φb can be any iteration of order b, we’ll use Householder’s variant in what follows. T[1,2] is the deltasquared operator which ‘improves’ given partial sums (or sequence of functional iterates) of order n toorder 2n − 1.T[2,2] is a bit of a surprise:µµ¶¶2F0F 00F01 F 02− 00T[2,2] (F ) = F + F − 00 +=F−F2F2 F 000(13.187)and gets order 2n2 − n out of n. A linearly convergent iteration will not be made super-linear by T[2,2] .In general, T[a,2]T[a,2] (F )µ¶iaXF (i)F0= F+− 00i!Fi=1=F−(13.188)1 F 02F 000 F 03F 0000 F 04−+± ...2 F 006 F 00324 F 004(13.189)takes order n to na−1 (2n − 1) = 2na − na−1 .T[a,b] also (only) takes order n to na−1 (2n − 1) = 2na − na−1 .

Linear convergence is turned into secondorder for a = b − 1.Here is what we get:a=1a=2a=3a=4a=5an=2361224482a+1 − 2a−1n=35154513512152 3a − 3a−1n=472811244817922 4a − 4a−1n2n − 12n2 − n2n3 − n22n4 − n32n5 − n42na − na−1Using a Padé approximant for the expansion in D and the general Padé type iteration Φ[] one can definethe extrapolation operator³´P[a,a0 ][b,b0 ] (F ):=Padé[a,a0 ] F (x + D)(13.190)whereD := Φ {F 0 }[b,b0 ] − F(13.191)P[1,0][2,0] is Aitken’s delta squared process.

For a = 1, a0 = 1, b = 2, b0 = 1 one getsP[1,1][2,1] (F ) ==2F 02 D=D − 2F 02F 02 F 00F−3F 002 − F 0 F 000F−F 00with D =2F 0 F 00− 2F 002F 0 F 000(13.192)(13.193)which makes linear convergence quadratic and takes super-linear convergence to 2n2 − n.13.9Transcendental functions & the AGM13.9.1The AGMThe AGM (arithmetic geometric mean) plays a central role in the high precision computation of logarithmsand π.CHAPTER 13. ARITHMETICAL ALGORITHMS304The AGM (a, b) is defined as the limit of the iteration AGM iteration, cf.13.194a :ak+1=bk+1=ak + bkp 2ak bk(13.194a)(13.194b)starting with a0 = a and b0 = b. Both of the values converge quadratically to a common limit. Therelated quantity ck (used in many AGM based computations) is defined asc2k= a2k − b2k= (ak−1 − ak )2(13.195)(13.196)One further defines (cf.

sections 13.9.5 and 13.9.6)"R0 (k):=∞1 X n 21−2 cn2 n=0#−1(13.197)where c2n := a2n − b2n corresponding to AGM (1, k).An alternative way for the computation for the AGM iteration isck+1=ak+1=bk+1=ak − bk2ak + bkq 2a2k+1 − c2k+1(13.198a)(13.198b)(13.198c)Schönhage gives the most economic variant of the AGM, which, apart from the square root, only needsone squaring per step:A0B0t0Skbkak+1Ak+1Bk+1c2k+1tk+1= a20= b20= 1 − (A0 − B0 )Ak + Bk=p 4=Bksquare root computationak + bk=2= a2k+1squaring√ ¶2√µ√Ak + BkAk + BkAk B k=+=242= 2 (Ak+1 − SK )= b2k+1= Ak+1 − Bk+1= tk −2k+1 c2k+1Starting with a0 = A0 = 1, B0 = 1/2 one has π ≈ (2 a2n )/tn .(13.199a)(13.199b)(13.199c)(13.199d)(13.199e)(13.199f)(13.199g)(13.199h)(13.199i)(13.199j)(13.199k)CHAPTER 13.

ARITHMETICAL ALGORITHMS305Combining two steps of the AGM iteration leads to the 4th order AGM iteration:√α0 =a0pβ0 =b0αk + βkαk+1 =2µ¶1/4αk βk (αk2 + βk2 )βk+1 =2γk4(Note that αk =√a2k and βk =√= αk4 − βk4(13.200b)(13.200c)(13.200d)(13.200e)b2k .) and"0= c2k/2(13.200a)R (k) =1−∞XÃn4αn4µ−n=0αn2 + βn22¶2 !#−1(13.201)√corresponding to AGM 4(1, k) (cf. [8] p.17).An alternative formulation of the 4th order AGM iteration is:γk+1αk+1βk+1c2k/2 + 2 c2k/2+113.9.2αk − βk2αk + βk=2¡ 4¢1/44= αk+1 − γk+1¡¢24= αk−1− αk2 − γk2=(13.202a)(13.202b)(13.202c)(13.202d)logThe (natural) logarithm can be computed using the following relation (cf.

[8] p.221)¯¯¯log(x) − R0 (10−n ) + R0 (10−n x)¯log(x)n102(n−1)≈ R0 (10−n ) − R0 (10−n x)≤(13.203)(13.204)that holds for n ≥ 3 and x ∈] 12 , 1[. Note that the first term on the rhs. is constant and might be storedfor subsequent log-computations. See also section 13.10.[hfloat: src/tz/log.cc]If one has some efficient algorithm for exp() one can compute log() from exp() usingylog(d):= 1 − d e−x= x + log(1 − y)¡¢= x + log 1 − (1 − d e−x ) = x + log (e−x d) = x + (−x + log(d))(13.205)(13.206)(13.207)Thenlog(d)¶µy3y2++ ...= x + log (1 − y) = x − y +23Truncation of the series after the n-th power of y gives an iteration of order n + 1:¶µy3y n−1y2++ ··· +xk+1 = Φn (xk ) := x − y +23n−1(13.208)(13.209)CHAPTER 13. ARITHMETICAL ALGORITHMS306Padé series P[i,j] (z) of log (1 − z) at z = 0 produce (order i + j + 2) iterations. For i = j we get[i, j] →7[0, 0] →7[1, 1] 7→[2, 2] 7→[4, 4] 7→x + P[i,j] (z = 1 − d e−x )x−z6−zx−z·6 − 4z30 − 21z + z 2x−z·30 − 36z + 9z 23780 − 6510z + 3360z 2 − 505z 3 + 6z 4x−z·3780 − 8400z + 6300z 2 − 1800z 3 + 150z 4(13.210a)(13.210b)(13.210c)(13.210d)(13.210e)Compared to the power series based iteration one needs one additional long division but saves half of theexponentiations.

This can be a substantial saving for high order iterations.13.9.3expThe exponential function can be computed using the iteration that is obtained as follows:exp(d) = x exp (d − log(x))= x exp(y) where y := d − log(x)µ¶y2y3= x 1+y+++ ...23!(13.211)(13.212)(13.213)The corresponding n-th oder iteration isµxk+1=Φn (xk ) := xk1+y+y3y n−1y2++ ...23!(n − 1)!¶(13.214)As the computation of logarithms is expensive one should use a higher (e.g.

8th) order iteration.[hfloat: src/tz/itexp.cc]Padé series P[i,j] (z) of exp (z) at z = 0 produce (order i + j + 1) iterations. For i = j we get[i, j] 7→[1, 1] 7→[2, 2] 7→[4, 4] 7→x P[i,j] (z = d − log x)z+2x·z−212 + 6z + z 2x·12 − 6z + z 21680 + 840z + 180z 2 + 20z 3 + z 4x·1680 − 840z + 180z 2 − 20z 3 + z 4The [i, j]-th Padé approximant of exp(z) is( i ¡ ¢) ( j ¡ ¢)X i zkX j (−z)kkk¡i+j ¢¡i+j ¢P[i,j] (z) =/k!k!kkk=0(13.215a)(13.215b)(13.215c)(13.215d)(13.216)k=0The numerator for i = j (multiplied by (2i)!/i! in order to avoid rational coefficients) is=¡i¢i(2i)! X k z k¡2 i¢·i!k!kk=0(13.217)CHAPTER 13. ARITHMETICAL ALGORITHMS13.9.4307sin, cos and tanFor arcsin, arccos and arctan use the complex analogue of the AGM.

For sin, cos and tan use the expiteration above think complex.13.9.5Elliptic KThe function K can be defined asZK(k) =π/2K(k)==21 − k 2 sin ϑ0One hasZdϑp1=0dtp2(1 − t ) (1 − k 2 t2 )µ1 1¯ ¶π, ¯F 2 2 ¯ k221¶2∞ µπ X (2 i − 1)!!2=π2=π2(13.218)(13.219)k 2ii i!2i=0Ã!µ ¶2µ¶2µ¶211·31·3·52461+k +k +k + ...22·42·4·6µ¶19 425 61225 83969 101 + k2 +k +k +k +k + ...4642561638465536(13.220)(13.221)(13.222)A special value isπ2(13.223)ππ=2 AGM (1, k 0 )2 AGM (1, 1 − k 2 )(13.224)π2 AGM (1, k)(13.225)K(0) =The computational interesting form isK(k)=One defines k 0 = 1 − k 2 and K 0 asK 0 (k) :=K(k 0 ) = K(1 − k 2 ) =[hfloat: src/tz/elliptick.cc]Product forms for K and K 0 that are also candidates for fast computations are2 0K (k0 ) =π2 0K (k0 ) =π∞Y∞Y2=1 + kn01+knn=0n=1√2 knwhere kn+1 :=,1 + kn∞Y1√knn=00 < k0 ≤ 11 + kn√ ,2 kn0 < k0 ≤ 1wherekn+1 :=(13.226)(13.227)These follow directly from 13.225 (and AGM (a, b) = a AGM (1, b/a)):³11 + k √ ´−1=AGM (, k)AGM (1, k)2√³1 + k2 k ´−1AGM (1,)=(first form)21+k³√1 + k ´−1=k AGM ( √ , 1)(second form)2 k(13.228)CHAPTER 13.

ARITHMETICAL ALGORITHMS308The second form seems computationally more effective. Similarly,2K(k0 ) =π∞Y∞Y2=1 + kn1 + kn0n=0n=1where kn+1(13.229)p1 − 1 − kn21 − kn0p=:=,1 + kn01 + 1 − kn20 < k0 ≤ 1With an efficient algorithm for K the logarithm can be computed using¯¯¯ 0¯¯K (k) − log ( 4 )¯ ≤ 4k 2 (8 + | log k|) where 0 < k ≤ 1¯k ¯13.9.6(13.230)Elliptic EThe function E can be defined asZπ/2E(k) =Zp221 − k sin ϑ d ϑ =010√1 − k 2 t2√dt1 − t2(13.231)One hasE(k) ====µ 1 1¯ ¶π−2, 2 ¯ 2F¯k21Ã!¶2∞ µX(2 i − 1)!!k 2iπ1−22i i!2i − 1i=0Ã!µ ¶2µ¶2 4 µ¶2 611·3k1·3·5kπ21−k −−− ...222·432·4·65µ¶π13 45 6175 8441 101 − k2 −k −k −k −k − ...24642561638465536(13.232)(13.233)(13.234)(13.235)Special values areE(0) =π2E(1) = 1(13.236)(The latter leads to a series for 2/π).