20.6 Arithmetic at Arbitrary Precision 915 CITED REFERENCES AND FURTHER READING: Bell,T.C.,Cleary,J.G.,and Witten,I.H.1990,Text Compression (Englewood Cliffs,NJ:Prentice- Hall). Nelson,M.1991.The Data Compression Book(Redwood City.CA:M&T Books). Witten,I.H.,Neal,R.M.,and Cleary,J.G.1987,Communications of the ACM,vol.30,pp.520- 540.[1] 20.6 Arithmetic at Arbitrary Precision 83g nted for Let's compute the number a to a couple of thousand decimal places.In doing so,we'll learn some things about multiple precision arithmetic on computers and Cam meet quite an unusual application of the fast Fourier transform(FFT).We'll also develop a set of routines that you can use for other calculations at any desired level of arithmetic precision. RECIPESI To start with,we need an analytic algorithm for m.Useful algorithms are quadratically convergent,i.e.,they double the number of significant digits at each iteration.Quadratically convergent algorithms for are based on the AGM (arithmetic geometric mean)method,which also finds application to the calculation of elliptic integrals(cf.$6.11)and in advanced implementations of the ADI method for elliptic partial differential equations($19.5).Borwein and Borwein [1]treat this subject,which is beyond our scope here.One of their algorithms for starts with 。多2 the initializations IENTIFIC X0=V2 61 T0=2+V2 (20.6.1) %=迈 (ISBN and then,fori=0,1,...,repeats the iteration Numerica 10.621 (+左) Recipes 43108 Xi+1= Xi+1+1 (outside Ti+1= Y+1/ (20.6.2) North Software. YivXi+1+- Y+1= Xi+1 Y+1 The value m emerges as the limit Too. Now,to the question of how to do arithmetic to arbitrary precision:In a high-level language like C,a natural choice is to work in radix(base)256,so that character arrays can be directly interpreted as strings of digits.At the very end of our calculation,we will want to convert our answer to radix 10,but that is essentially a frill for the benefit of human ears,accustomed to the familiar chant,"three point20.6 Arithmetic at Arbitrary Precision 915 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). CITED REFERENCES AND FURTHER READING: Bell, T.C., Cleary, J.G., and Witten, I.H. 1990, Text Compression (Englewood Cliffs, NJ: Prentice￾Hall). Nelson, M. 1991, The Data Compression Book (Redwood City, CA: M&T Books). Witten, I.H., Neal, R.M., and Cleary, J.G. 1987, Communications of the ACM, vol. 30, pp. 520– 540. [1] 20.6 Arithmetic at Arbitrary Precision Let’s compute the number π to a couple of thousand decimal places. In doing so, we’ll learn some things about multiple precision arithmetic on computers and meet quite an unusual application of the fast Fourier transform (FFT). We’ll also develop a set of routines that you can use for other calculations at any desired level of arithmetic precision. To start with, we need an analytic algorithm for π. Useful algorithms are quadratically convergent, i.e., they double the number of significant digits at each iteration. Quadratically convergent algorithms for π are based on the AGM (arithmetic geometric mean) method, which also finds application to the calculation of elliptic integrals (cf. §6.11) and in advanced implementations of the ADI method for elliptic partial differential equations (§19.5). Borwein and Borwein [1] treat this subject, which is beyond our scope here. One of their algorithms for π starts with the initializations X0 = √ 2 π0 =2+ √ 2 Y0 = √4 2 (20.6.1) and then, for i = 0, 1,... , repeats the iteration Xi+1 = 1 2
Xi + 1 √Xi πi+1 = πi
Xi+1 + 1 Yi + 1 Yi+1 = Yi Xi+1 + 1 Xi+1 Yi + 1 (20.6.2) The value π emerges as the limit π∞. Now, to the question of how to do arithmetic to arbitrary precision: In a high-level language like C, a natural choice is to work in radix (base) 256, so that character arrays can be directly interpreted as strings of digits. At the very end of our calculation, we will want to convert our answer to radix 10, but that is essentially a frill for the benefit of human ears, accustomed to the familiar chant, “three point