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

The key to fast computations is1−EK= 1−R=∞1 X n 022 cn2 n=0(13.237)(see [7] p.8) and soR0 (k)E(k) =:=K 0 (k)=E 0 (k)R0 (k) K(k) ="∞1 X n 21−2 cn2 n=02 AGM (1, 1 −k2 )#−1πP∞· (1 − n=0 2n−1 c2n )(13.238)(13.239)Similar as for K 0 one definesE 0 (k):=E(k 0 ) = E(1 − k 2 )(13.240)CHAPTER 13. ARITHMETICAL ALGORITHMS309An ellipse with major axis pa and minor axis b has the circumference (or arclength) L = 4 a E 0 (b/a) =4 a E(e). The quantity e = 1 − a2 /b2 is the eccentricity of the ellipse.Legendre’s relation between K and E is (arguments omitted for readability, choose your favorite form):E0E+ 0 −1 =KKE K0 + E0 K − K K0π2 K K0π2=(13.241)(13.242)Equivalently,AGM (1, k)For k =√12RE/K=(1 − E 0 /K 0 )(1 − R0 )=(13.243)=: s we have k = k 0 , thereby K = K 0 and E = E 0 , soK(s)πµ2 E(s) K(s)−ππ¶=12πAs expressions 13.224 and 13.239 provide a fast AGM based computation ofcan be used to compute π (cf. [8]).(13.244)KπandEπthe above formula2Using E − K = k k 0 dKdk − k K one can express the derivative of K in terms of E and K and therebycompute that quantity fast:dKdk13.10E − k 02 Kk k 02=(13.245)Computation of π/ log(q)For the computation of the natural logarithm one can use the relationlog(m rx )=log(m) + x log(r)(13.246)where m is the mantissa and r the radix of the floating point numbers.There is a nice way to compute the value of log(r) if the value of π has been precomputed.

One definesθ2 (q) =θ3 (q) =θ4 (q) =∞Xn=−∞∞Xn=−∞∞Xq (n+1/2)qn22(13.247)(13.248)2(−1)n q n(13.249)n=−∞Thenlog1q= − log q = πK0K(13.250)(so q = exp(−π K 0 /K)) and (cf. [8] p.225)πlog(1/q)=−¡¢π= AGM θ3 (q)2 , θ2 (q)2log(q)(13.251)CHAPTER 13. ARITHMETICAL ALGORITHMS310Computing θ3 (q) is easy when q = 1/r:θ3 (q)=∞X1+22q n = 2 (1 +n=1∞X2qn ) − 1(13.252)n=1However, the computation of θ2 (q) suggests to choose q = 1/r4 =: b4 :θ2 (q) ==0+22b∞X2q (n+1/2) = 2n=0∞Xn2 +nq= 2 b (1 +n=0∞Xn=0∞X2b4nqn2+4n+1+nwhere q = b4)(13.253)(13.254)n=1[hfloat: src/tz/pilogq.cc]One has (see [22])Θ22 (q) =2kKπΘ23 (q) =Θ22 (q)Θ23 (q)k =πlog(q)Θ24 (q) =k0 =2 k0 KπΘ24 (q)Θ23 (q)AGM (1, k)AGM (1, k 0 )¢¡AGM Θ23 (q), Θ24 (q)=−1 =13.112Kπ(13.255)(13.256)(13.257)(13.258)Computation of q = exp(−π K 0 /K)We use (cf. [7] p.35/36)q=µ¶K0exp −πK(13.259)Following the usual convention, we writeK0Kwhere k 0 = b0 and k = b00 ==AGM (1, k 0 )AGM (1, b0 )=AGM (1, k)AGM (1, b0 0 )(13.260)p1 − b20 and use ([7] p.38, note the missing ‘4’ there)π AGM (1, b0 )2 AGM (1, b0 0 )=limn→∞14 anlogn2cn(13.261)thereby¶µ¶4 an14 an1= lim exp − n−1 logexp −2 lim n logn→∞n→∞ 2cn2cnµ¶−1/2n−1µ¶−1/(2n−1 )4 an4 anlim exp log= limn→∞n→∞cncnµq==(13.262)(13.263)soµq=limn→∞cn4 an¶1/(2n−1 )(13.264)CHAPTER 13.

ARITHMETICAL ALGORITHMS311One obtains an algorithm for exp(−x) by first solving for k, k 0 so that x = π K 0 /K (precomputed π) andapplying the last relation that implies the computation of a 2n−1 -th root. Note that the quantity c shouldc2nbe computed via cn+1 = 4 an+1throughout the AGM computation in order to preserve its accuracy. (cf.also [8] p.227)√For k = 1/ 2 =: s one has k = k 0 and so q = exp(−π). Thus the calculation of exp(−π) =0.0432139182637 . .

. can directly be done via a single AGM computation as (cn /(4an ))N where N =1/2(n−1) . The quantity ii = exp(−π/2) = 0.2078795763507 . . . can be obtained using N = 1/2n . (cf. also[7] p.13)13.12Iterations for high precision computations of πIn this section various iterations for computing π with at least second order convergence are given.The number of full precision multiplications (FPM) are an indication of the efficiency of the algorithm.The approximate number of FPMs that were counted with a computation of π to 4 million decimaldigits12 is indicated like this: #FPM=123.4.AGM as in [hfloat: src/pi/piagm.cc], #FPM=98.4 (#FPM=149.3 for the quartic variant):a0= 1(13.265a)b0=1√2pn=2 a2Pnn+1 k 21 − k=0 2 ckπ − pn=π 2 2n+4 e−π 2AGM 2 (a0 , b0 )(13.265b)→π(13.265c)n+1(13.265d)A fourth order version uses 13.200a, cf.

also [hfloat: src/pi/piagm.cc].AGM variant as in [hfloat: src/pi/piagm3.cc], #FPM=99.5 (#FPM=155.3 for the quartic variant):a0= 1√b0=pn=√√π − pn<6+4√(13.266a)2(13.266b)2 a2n+1Pn3 (1 − k=0 2k c2k ) − 1√3 π 2 2n+4 e− 3 π 2AGM 2 (a0 , b0 )→π(13.266c)n+1(13.266d)AGM variant as in [hfloat: src/pi/piagm3.cc], #FPM=108.2 (#FPM=169.5 for the quartic variant):a0= 1√b0=pn=π − pn12 usingradix 10, 000 and 1 million LIMBs.<√6−4√(13.267a)2(13.267b)6 a2n+1Pn3 (1 − k=0 2k c2k ) + 1√13π 2 2n+4 e→π(13.267c)− √13 π 2n+1AGM (a0 , b0 )2(13.267d)CHAPTER 13. ARITHMETICAL ALGORITHMS312Borwein’s quartic (fourth order) iteration, variant r = 4 as in [hfloat: src/pi/pi4th.cc], #FPM=170.5:y0a0yk+1=2−1√= 6−4 2==ak+10√=1 − (1 − yk4 )1/41 + (1 − yk4 )1/4(1 − yk4 )−1/4 − 1(1 − yk4 )−1/4 + 1(13.268a)(13.268b)→0+(13.268c)(13.268d)2)ak (1 + yk+1 )4 − 22k+3 yk+1 (1 + yk+1 + yk+1→1π= ak ((1 + yk+1 )2 )2 − 22k+3 yk+1 ((1 + yk+1 )2 − yk+1 )n< ak − π −1 ≤ 16 · 4n 2 e−4 2 π(13.268e)(13.268f)(13.268g)Identities 13.268d and 13.268f show how to save operations.Borwein’s quartic (fourth order) iteration, variant r = 16 as in [hfloat: src/pi/pi4th.cc], #FPM=164.4:1 − 2−1/41 + 2−1/4√8/ 2 − 2(2−1/4 + 1)4y0=a0=yk+1=(1 − yk4 )−1/4 − 1(1 − yk4 )−1/4 + 1ak+1=2ak (1 + yk+1 )4 − 22k+4 yk+1 (1 + yk+1 + yk+1)0(13.269a)(13.269b)→0+< ak − π −1 ≤ 16 · 4n 4 e−4n(13.269c)→1π4π(13.269d)(13.269e)Same operation count as before, but this variant gives approximately twice as much precision after thesame number of steps.The general form of the quartic iterations (13.268a and 13.269a) ispy0 =λ∗ (r)a0 = α(r)yk+1=(1 − yk4 )−1/4 − 1(1 − yk4 )−1/4 + 1→0+√2= ak (1 + yk+1 )4 − 22k+2 r yk (1 + yk+1 + yk+1)√√n0 < ak − π −1 ≤ 16 · 4n r e−4 r πak+1Cf.

[8], p.170f.(13.270a)(13.270b)(13.270c)→1π(13.270d)(13.270e)CHAPTER 13. ARITHMETICAL ALGORITHMS313Derived AGM iteration (second order) as in [hfloat: src/pi/pideriv.cc], #FPM=276.2:√x0 =2√p0 = 2 + 2y1 = 21/4µ¶1 √1xk+1 =xk + √(k ≥ 0) → 1 +2xk√yk xk + √1xkyk+1 =(k ≥ 1) → 1 +yk + 1xk + 1pk+1 = pk(k ≥ 1) → π +yk + 1pk − π= 10−2k+1(13.271a)(13.271b)(13.271c)(13.271d)(13.271e)(13.271f)(13.271g)Cubic AGM from [46], as in [hfloat: src/pi/picubagm.cc], #FPM=182.7:a0=b0=an+1=bn+1=pn=1√(13.272a)3−12an + 2 bnr 3223 bn (an + an bn + bn )33 a2nPn1 − k=0 3k (a2k − a2k+1 )(13.272b)(13.272c)(13.272d)(13.272e)Second order iteration, as in [hfloat: src/pi/pi2nd.cc], #FPM=255.7:y0=a0=yk+1==ak+1ak − π −11√2121 − (1 − yk2 )1/21 + (1 − yk2 )1/2(13.273a)(13.273b)→0+(13.273c)(1 − yk2 )−1/2 − 1(1 − yk2 )−1/2 + 1= ak (1 + yk+1 )2 − 2k+1 yk+1≤ 16 · 2k+1 e−2k+1π13.273d shows how to save 1 multiplication per step (cf.

section 13.3).(13.273d)→1π(13.273e)(13.273f)CHAPTER 13. ARITHMETICAL ALGORITHMS314Quintic (5th order) iteration from the article [43], as in [hfloat: src/pi/pi5th.cc], #FPM=353.2:√s0 = 5( 5 − 2)(13.274a)1a0 =(13.274b)225sn+1 =→1(13.274c)sn (z + x/z + 1)25where x =−1 →4(13.274d)snand y = (x − 1)2 + 7 → 16(13.274e)³x ³´´1/5pand z =y + y 2 − 4x3→2(13.274f)2µ 2¶sn − 5 p1an+1 = s2n an − 5n+ sn (s2n − 2sn + 5)→(13.274g)2πn1an −< 16 · 5n e−π 5(13.274h)πCubic (third order) iteration from [44], as in [hfloat: src/pi/pi3rd.cc], #FPM=200.3:13√a0=(13.275a)s0=rk+1=sk+1=ak+122= rk+1ak − 3k (rk+1− 1)3−12(13.275b)31 + 2 (1 − s3k )1/3rk+1 − 12(13.275c)(13.275d)→1π(13.275e)Nonic (9th order) iteration from [44], as in [hfloat: src/pi/pi9th.cc], #FPM=273.7:a0=r0=s0tuvm13√(13.276a)3−12= (1 − r03 )1/3= 1 + 2 rk¡¢1/3= 9 rk (1 + rk + rk2 )2(13.276b)(13.276c)(13.276d)(13.276e)2= t + tu + u27 (1 + sk + s2k )=v(13.276f)(13.276g)1πak+1= m ak + 32 k−1 (1 − m)sk+1=(1 − rk )3(t + 2 u) v(13.276i)rk+1= (1 − s3k )1/3(13.276j)Summary of operation count vs.

algorithms:→(13.276h)CHAPTER 13. ARITHMETICAL ALGORITHMS#FPM78.42498.42499.510108.241149.324155.265164.359169.544170.519182.710200.261255.699273.763276.221353.202-315algorithm name in hfloatpi_agm_sch()pi_agm()pi_agm3(fast variant)pi_agm3(slow variant)pi_agm(quartic)pi_agm3(quartic, fast variant)pi_4th_order(r=16 variant)pi_agm3(quartic, slow variant)pi_4th_order(r=4 variant)pi_cubic_agm()pi_3rd_order()pi_2nd_order()pi_9th_order()pi_derived_agm()pi_5th_order()TBD: notes: discontin.TBD: slow quartic, slow quart.AGMTBD: other quant: num of variablesMore iterations for πThese are not (yet) implemented in hfloat.A third order algorithm from [45]:v0= 2−1/8v1= 2−7/8w0α0β0= 1= 1= 0vn+1wn+1αn+1βn+1πn³1/2(1 − 3)2−1/2+31/4´n£¤1/3 o1/2= vn3 − vn6 + 4vn2 (1 − vn8 )+ vn−1¡ 2¢322vn + vn+1 3vn+1 vn − 1¡ 2¢ wn=32vn+1− vn 3vn+1vn2 − 1¶µ 32vn+1+ 1 αn=vnµ 3¶v 2 αn2vn+1=+ 1 βn + (6wn+1 vn − 2vn+1 wn ) n+12vnvn=8 · 21/8αn βn→π(13.277a)(13.277b)(13.277c)(13.277d)(13.277e)(13.277f)(13.277g)(13.277h)(13.277i)(13.277j)CHAPTER 13.

ARITHMETICAL ALGORITHMS316A second order algorithm from [47]:α0m0mn+1αn+1= 1/3= 2=(13.278a)(13.278b)4p1+(13.278c)(4 − mn ) (2 + mn )2n1(1 − mn ) →= mn αn +3π(13.278d)Another second order algorithm from [47]:α0s12∗ 2(sn ) + (sn )(1 + 3 sn+1 ) (1 + 3 s∗n )αn+1= 1/3= 1/3= 1= 4(13.279a)(13.279b)(13.279c)(13.279d)(1 + 3 sn+1 )αn − 2n sn+1=→1π(13.279e)A fourth order algorithm from [47]:α0s1∗ 4(sn ) + (sn )3 sn+1 ) (1 + 3 s∗n )4(1 +αn+113.13= 1/3√=2−1= 1= 2(13.280a)(13.280b)(13.280c)(13.280d)= (1 + sn+1 )4 αn +4n+14(1 − (1 + sn+1 ) )3→1π(13.280e)The binary splitting algorithm for rational seriesThe straight forward computation of a series for which each term adds a constant amount of precision13to a precision of N digits involves the summation of proportional N terms.