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

To get N bits of precision onehas to add proportional N terms of the sum, each term involves one (length-N ) short division (and oneaddition). Therefore the total work is proportional N 2 , which makes it impossible to compute billions ofdigits from linearly convergent series even if they are as ‘good’ as Chudnovsky’s famous series for π:¶µ¶∞ µ16541681608 X 13591409(6k)!(−1)k= √+k(13.281)3π(k!)3 (3k)! 6403203k640320 k=0 545140134=√∞X126403203(−1)kk=0Here is an alternative way to evaluate a sumof consecutive terms:(6k)!13591409 + k 5451401343(k!) (3k)!(640320)3kPN −1k=0(13.282)ak of rational summands: One looks at the ratios rkakak−1(13.283)=: r0 (1 + r1 (1 + r2 (1 + r3 (1 + . .

. (1 + rN −1 ) . . . ))))(13.284)rk:=(set a−1 := 1 to avoid a special case for k = 0)That isN−1Xakk=013 e.g.arc-cotangent series with arguments > 1CHAPTER 13. ARITHMETICAL ALGORITHMS317Now definerm,nrm,m:= rm (1 + rm+1 (. . . (1 + rn ) . . . )):= rmwhere m < n(13.285)(13.286)thenrm,n=nX1am−1ak(13.287)k=mand especiallyr0,n=nXak(13.288)k=0Withrm,n= rm + rm · rm+1 + rm · rm+1 · rm+2 + . . .· · · + rm · · · · · rx + rm · · · · · rx · [rx+1 + · · · + rx+1 · · · · · rn ]xY= rm,x +rk · rx+1,n(13.289)(13.290)k=mThe product telescopes, one getsrm,n=rm,x +ax· rx+1,nam−1(13.291)(where m ≤ x < n).Now we can formulate the binary splitting algorithm by giving a binsplit function r:function r(function a, int m, int n){rational ret;if m==n then{ret := a(m)/a(m-1)}else{x := floor( (m+n)/2 )ret := r(a,m,x) + a(x) / a(m-1) * r(a,x+1,n)}print( "r:", m, n, "=", ret )return ret}Here a(k) must be a function that returns the k-th term of the series we wish to compute, in additionone must have a(-1)=1.

A trivial example: to compute arctan(1/10) one would usefunction a(int k){if k<0 then return 1elsereturn (-1)^k/((2*k+1)*10^(2*k+1))}Calling r(a,0,N) returnsPNk=0ak .In case the programming language used does not provide rational numbers one needs to rewrite formulaUm,n13.291 in separate parts for denominator and numerator. With ai = pqii , p−1 = q−1 = 1 and rm,n =: Vm,none getsUm,n=pm−1 qx Um,x Vx+1,n + px qm−1 Ux+1,n Vm,x(13.292)Vm,n=pm−1 qx Vm,x Vx+1,n(13.293)CHAPTER 13.

ARITHMETICAL ALGORITHMS318The reason why binary splitting is better than the straight forward way is that the involved work is onlyO((log N )2 M (N )), where M (N ) is the complexity of one N -bit multiplication (see [42]). This meansthat sums of linear but sufficient convergence are again candidates for high precision computations.In addition, the ratio r0,N −1 (i.e.

the sum of the first N terms) can be reused if one wants to evaluatethe sum to a higher precision than before. To get twice the precision user0,2 N −1=r0,N −1 + aN −1 · rN,2 N −1(13.294)(this is formula 13.291 with m = 0, x = N − 1, n = 2N − 1). With explicit rational arithmetic:U0,2N −1V0,2N −1==qN −1 U0,N −1 VN,2N −1 + pN −1 UN,2N −1 V0,N −1qN −1 V0,N −1 VN,2N −1(13.295)(13.296)Thereby with the appearance of some new computer that can multiply two length 2·N numbers14 one onlyneeds to combine the two ratios r0,N −1 and rN,2N −1 that had been precomputed by the last generationof computers.

This costs only a few full-size multiplications on your new and expensive supercomputer(instead of several hundreds for the iterative schemes), which means that one can improve on priorcomputations at low cost.If one wants to stare at zillions of decimal digits of the floating point expansion then one division is alsoneeded which costs not more than 4 multiplications (cf. section 13.3).Note that this algorithm can trivially be extended (or rather simplified) to infinite products, e.g. matrixproducts as Bellard’s"#·¸∞2 (k− 21 ) (k+2)Y0 π+6102427 (k+ 3 ) (k+ 3 )=(13.297)0101k=0Cf.

[42] and [48].13.14The magic sumalt algorithmThe following algorithm is due to Cohen, Villegas and Zagier, see [49].P∞Pseudo code to compute an estimate of k=0 xk using the first n summands. The xk summands areexpected in x[0,1, ...,n-1].function sumalt(x[], n){d := (3+sqrt(8))^nd := (d+1/d)/2b := 1c := ds := 0for k:=0 to n-1{c := c - bs := s + c * x[k]b := b * (2*(n+k)*(n-k)) / ((2*k+1)*(k+1))}return s/d}With alternating sums the accuracy of the estimate will be (3 +√8)−n ≈ 5.82−n .As an example let us explicitely write down the estimate for the 4 · arctan(1) using the first 8 termsµ¶1 1 1 1 1111π ≈ 4·− + − + −+−= 3.017 .

. .(13.298)1 3 5 7 9 11 13 1514 assumingone could multiply length-N numbers beforeCHAPTER 13. ARITHMETICAL ALGORITHMS319The sumalt-massaged estimate isµ665856 665728 663040 641536π ≈ 4·−+−+1357¶557056 376832 163840 32768+−+−/6658579111315= 4 · 3365266048/4284789795 = 3.141592665 . . .(13.299)and already gives 7 correct digits of π. The linear but impressive growth of the successive sumalt estimateswith n, the number of terms used, is illustrated by the following table:n1234567891011121314151617181920sumalt(n)2.6666666666666666666663.1372549019607843137253.1407407407407407407403.1416357184121482215073.1415865464036739683483.1415933442156594036603.1415925649375401220153.1415926652243158640173.1415926520088119516193.1415926538097315693183.1415926535585787555133.1415926535942963384703.1415926535891345805173.1415926535898907186253.1415926535897786643753.1415926535897954367753.1415926535897929042853.1415926535897932896143.1415926535897932305843.141592653589793239682sumalt(n)−π0.4749259869231265717950.0043377516290089247370.000851912849052497721-0.0000430648223549830440.000006107186119270114-0.0000006906258661651970.000000088652253116447-0.0000000116345226255550.000000001580981286843-0.0000000002199383308560.000000000031214482948-0.0000000000045031000070.000000000000658657944-0.0000000000000974801630.000000000000014574087-0.0000000000000021983120.000000000000000334177-0.0000000000000000511510.000000000000000007877-0.000000000000000001220Therefore even slowly converging series likeπClog(2)===4·∞X(−1)k2k + 1k=0∞Xk=0∞Xk=0ζ(s)== 4 · arctan(1)(−1)k(2 k + 1)2(−1)kk+1= 0.9159655941772190 .

. .= 0.6931471805599453 . . .∞X1(−1)k1 − 21−sks(13.300)(13.301)(13.302)(13.303)k=1can be used to compute estimates that are correct up to thousands of digits. The algorithm scales like n2if the series terms in x[] are small rational values and like n3 · log(n) if they are full precision (rationalor float) values.All values ck and bk occurring in the computation are integers.

In fact, the bk in the computation with nterms are the coefficients of the expanded n-th Chebychev polynomial (of the first kind) with argument1 + 2x:CHAPTER 13. ARITHMETICAL ALGORITHMSk012345678bk320ck6658576658566657286630406415365570563768321638403276811282688215048448018022421299213107232768T8 (1 + 2x) = 1 + 128x + 2688x2 + 21504x3 + 84480x4 +√+180224x5 + 212992x6 + 131072x7 + 32768x8 = T16 ( 1 + x)T16 (x) = 1 − 128x2 + 2688x4 − 21504x6 + 84480x8 −−180224x10 + 212992x12 − 131072x14 + 32768x16(13.304)(13.305)Now observe that one has always cn = bn =√22n−1 in a length-n sumalt computation. Obviously, ‘goingbackwards’ avoids the computation of (3 + 8)n :function sumalt(x[], n){b := 2**(2*n-1)c := bs := 0for k:=n-1 to 0 step -1{s := s + c * x[k]b := b * ((2*k+1)*(k+1)) / (2*(n+k)*(n-k))c := c + b}return s/c}The bk and ck occurring in a length-n sumalt computation can be given explicitely asµ¶nn + i 2kbk =2n + k 2kµ¶knXXn + i 2inck = dn −ck =2n+i2ii=0(13.306)(13.307)i=k+1P∞To compute an estimate of k=0 xk using the first n partial sums use the following pseudo code (thePkpartial sums pk = j=0 xj are expected in p[0,1,...,n-1]):function sumalt_partial(p[], n){d := (3+sqrt(8))^nd := (d+1/d)/2b := 1c := ds := 0for k:=0 to n-1{s := s + b * p[k]b := b * (2*(n+k)*(n-k)) / ((2*k+1)*(k+1))}return s/d}The backward variant is:CHAPTER 13.

ARITHMETICAL ALGORITHMS321function sumalt_partial(p[], n){b := 2**(2*n-1)c := bs := 0for k:=n-1 to 0 step -1{s := s + b * p[k]b := b * ((2*k+1)*(k+1)) / (2*(n+k)*(n-k))c := c + b}return s/c}Cf. [hfloat: src/hf/sumalt.cc]For series of already geometrical rate of convergence [49] givesfunction sumalt_partial(p[], n, e){d := ( 2*e + 1 + 2*sqrt(e*(e+1)) )^nd := (d+1/d)/2b := 1c := ds := 0for k:=0 to n-1{s := s + b * p[k]b := b * (2*(n+k)*(n-k)) / ((2*k+1)*(k+1)) * e}return s/d}pfor series where |ak /ak+1 | ≈ e. Convergence is improved from ∼ e−n to ∼ (2e + 1 + 2 e(e + 1))−n ≈(4e + 2)−n .

This algorithm specializes to the original one for e = 1.13.15Chebyshev polynomials *The Chebychev polynomials of the first and second kind can be defined by the functionsTn (x)=cos(n arccos(x))sin((n + 1) arccos(x))√1 − x2Un (x) =(13.308)(13.309)For integral n both of them are polynomials. We haveT−n (x)= Tn (x)T−1 (x) = xT0 (x) = 1T1 (x)T2 (x)(13.310)(13.311)(13.312)= x= 2 x2 − 1(13.313)(13.314)T3 (x) = 4 x3 − 3 xT4 (x) = 8 x4 − 8 x2 + 1(13.315)(13.316)T5 (x) = 16 x5 − 20 x3 + 5 xT6 (x) = 32 x6 − 48 x4 + 18 x2 − 1(13.317)(13.318)T7 (x)=64 x7 − 112 x5 + 56 x3 − 7 x(13.319)CHAPTER 13. ARITHMETICAL ALGORITHMS322bn/2cn X (−1)k (n − k − 1)! (2x)n−2k2k! (n − 2k)!Tn (x) =(13.320)k=0bn/2c µ¶n(2x)n−2k (x2 − 1)k2k(13.321)−Un−2 (x)012x4 x2 − 18 x3 − 4 x16 x4 − 12 x3 + 132 x5 − 32 x3 + 6 x64 x6 − 80 x4 + 24 x2 − 1128 x7 − 192 x5 + 80 x3 − 8 x(13.322)(13.323)(13.324)(13.325)(13.326)(13.327)(13.328)(13.329)(13.330)(13.331)X=k=0andU−n (x)U−1 (x)U0 (x)U1 (x)U2 (x)U3 (x)U4 (x)U5 (x)U6 (x)U7 (x)==========bn/2cX (−1)k (n − k)! (2x)n−2kk! (n − 2k)!Un (x) =(13.332)k=0bn/2+1c µX=k=0¶n+1xn−2k (x2 − 1)k2k + 1(13.333)Both obey the same recurrence (omitting argument x)TnUn= 2 x Tn−1 − Tn−2= 2 x Un−1 − Un−2(13.334)(13.335)Their generating functions are1 − xtt2 − 2 x t + 1=1t2 − 2 x t + 1=∞Xn=0∞Xtn Tn (x)(13.336)tn Un (x)(13.337)n=0Composition is multiplication of indices:Tn (Tm (x))=Tn m (x)(13.338)For example,T2n (x)= T2 (Tn (x)) = 2 Tn2 (x) − 1(13.339)Some relations between T and U are1(Un − Un−2 )2Tn=Un − x Un−1 = x Un−1 − Un−2 =(13.340)Tn+1=x Tn − (1 − x2 ) Un−1(13.341)U2n−1=2 Tn Un−1(13.342)CHAPTER 13.