f20-6

906 Chapter20. Less-NumericalAlgorithmsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).SUBROUTINE arcsum(iin,iout,ja,nwk,nrad,nc) INTEGER ja,nc,nrad,nwk,iin(*),iout(*) Usedby arcode. Addtheinteger jatotheradix nradmultiple-precisioninteger iin(nc..nwk) . Return the result in iout(nc..nwk) . INTEGER j,jtmp,karrykarry=0do 11j=nwk,nc+1,-1 jtmp=ja ja=ja/nradiout(j)=iin(j)+(jtmp-ja*nrad)+karryif (iout(j).ge.nrad) then iout(j)=iout(j)-nrad karry=1 else karry=0 endif enddo 11 iout(nc)=iin(nc)+ja+karry return END If radix-changing, rather than compression, is your primary aim (for example to convert an arbitrary file into printable characters) then you are of course free to set all the components of nfreqequal, say, to 1. CITED REFERENCES AND FURTHER READING: Bell,T.C.,Cleary,J.G.,andWitten,I.H.1990, TextCompression (EnglewoodCliffs,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 (arithmeticgeometric mean) method,which also finds applicationto the calculation of elliptic integrals (cf. §6.11)and in advancedimplementationsof 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 ) 20.6ArithmeticatArbitraryPrecision 907Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).and then, for i=0,1,..., repeats the iteration Xi+1=1 2/parenleftbigg/radicalbig Xi+1√Xi/parenrightbigg πi+1=πi/parenleftbiggXi+1+1 Yi+1/parenrightbigg Yi+1=Yi/radicalbig Xi+1+1/radicalbig 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-levellanguagelike FORTRAN,anaturalchoiceistoworkinradix(base)256,so thatcharacterarrayscanbedirectlyinterpretedasstringsofdigits. Attheveryendof ourcalculation,wewillwanttoconvertouranswertoradix10,butthatisessentially a frill for the benefit of human ears, accustomed to the familiar chant, “three pointonefouronefivenine ....” Foranyless frivolouscalculation,wewouldlikelynever leavebase256(orthethencetriviallyreachablehexadecimal,octal,orbinarybases). We will adopt the convention of storing digit strings in the “human” ordering, that is, with the first stored digit in an array being most significant, the last stored digit being least significant. The opposite convention would, of course, also bepossible. “Carries,” where we need to partition a number larger than 255 into a low-order byte and a high-order carry, present a minor programming annoyance, solved, in the routines below, by the use of FORTRAN’sEQUIVALENCE facility, and some initial testing of the order in which bytes are stored in a FORTRAN integer. It is easy at this point, following Knuth [2], to write a routine for the “fast” arithmetic operations: short addition (adding a single byte to a string), addition, subtraction, short multiplication (multiplying a string by a single byte), short division,ones-complementnegation;andacoupleofutilityoperations,copyingandleft-shifting strings. SUBROUTINE mpops(w,u,v) CHARACTER*1 w(*),u(*),v(*) Multipleprecision arithmetic operations done oncharacter strings, interpreted as radix256numbers. This routine collects the simpler operations. INTEGER i,ireg,j,n,ir,is,iv,ii1,ii2 CHARACTER*1 creg(4) SAVE ii1,ii2EQUIVALENCE (ireg,creg) Itisassumedthatwiththeaboveequivalence, creg(ii1) addresses thelow-order byteof ireg,a n d creg(ii2) addresses the next higher order byte. The values ii1andii2are set by an initial call to mpinit. ENTRY mpinit ireg=256*ichar(’2’)+ichar(’1’) do11j=1,4 Figure out the byte ordering. if (creg(j).eq.’1’) ii1=j if (creg(j).eq.’2’) ii2=j enddo 11 returnENTRY mpadd(w,u,v,n) Adds the unsigned radix 256 integers u(1:n)andv(1:n)yielding the unsigned integer w(1:n+1). ireg=0 do12j=n,1,-1 908 Chapter20. Less-NumericalAlgorithmsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).ireg=ichar(u(j))+ichar(v(j))+ichar(creg(ii2)) w(j+1)=creg(ii1) enddo 12 w(1)=creg(ii2) return ENTRY mpsub(is,w,u,v,n) Subtractstheunsignedradix256integer v(1:n)from u(1:n)yieldingtheunsignedinteger w(1:n). If the result is negative (wraps around), isis returned as −1;o t h e r wi s ei ti s returned as 0. ireg=256 do13j=n,1,-1 ireg=255+ichar(u(j))-ichar(v(j))+ichar(creg(ii2))w(j)=creg(ii1) enddo 13 is=ichar(creg(ii2))-1 return ENTRY mpsad(w,u,n,iv) Short addition: the integer iv(inthe range 0≤iv≤255) is added to the unsigned radix 256 integer u(1:n), yielding w(1:n+1). ireg=256*iv do14j=n,1,-1 ireg=ichar(u(j))+ichar(creg(ii2))w(j+1)=creg(ii1) enddo 14 w(1)=creg(ii2) returnENTRY mpsmu(w,u,n,iv) Shortmultiplication: theunsignedradix256integer u(1:n)ismultipliedbytheinteger iv (in the range 0≤iv≤255), yielding w(1:n+1). ireg=0 do15j=n,1,-1 ireg=ichar(u(j))*iv+ichar(creg(ii2))w(j+1)=creg(ii1) enddo 15 w(1)=creg(ii2) returnENTRY mpsdv(w,u,n,iv,ir) Short division: the unsigned radix 256 integer u(1:n)is divided by the integer iv(in the range 0≤iv≤255),yieldingaquotient w(1:n)andaremainder ir(with 0≤ir≤255). ir=0do 16j=1,n i=256*ir+ichar(u(j)) w(j)=char(i/iv)ir=mod(i,iv) enddo 16 returnENTRY mpneg(u,n) Ones-complement negate the unsigned radix 256 integer u(1:n). ireg=256 do17j=n,1,-1 ireg=255-ichar(u(j))+ichar(creg(ii2))u(j)=creg(ii1) enddo 17 return ENTRY mpmov(u,v,n) Move v(1:n)onto u(1:n). do18j=1,n u(j)=v(j) enddo 18 returnENTRY mplsh(u,n) Left shift u(2..n+1) onto u(1:n). do19j=1,n u(j)=u(j+1) 20.6ArithmeticatArbitraryPrecision 909Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).enddo 19 return END Full multiplicationof two digit strings, if done by the traditionalhand method, is not a fast operation: In multiplying two strings of length N, the multiplicand would be short-multiplied in turn by each byte of the multiplier, requiring O(N2) operationsinall. Wewillsee,however,that allthearithmeticoperationsonnumbers of length Ncan in fact be done in O(N×logN×log log N)operations. The trick is to recognizethat multiplicationis essentially a convolution (§13.1) of the digits of the multiplicand and multiplier, followed by some kind of carryoperation. Consider,forexample,two ways ofwritingthe calculation 456×789: 456 ×789 4104 3648 3192 359784456 ×789 36 45 54 32 40 48 28 35 42 28 67 118 93 54 359 784 The tableau on the left shows the conventional method of multiplication, in which three separate short multiplications of the full multiplicand (by 9, 8, and 7) are added to obtain the final result. The tableau on the right shows a different method (sometimes taught for mental arithmetic), where the single-digit cross products are all computed (e.g. 8×6=4 8), then added in columns to obtain an incompletely carried result (here, the list 28,67,118,93,54). The final step is a single pass from right to left, recordingthe single least-significant digit and carryingthe higher digit or digits into the total to the left (e.g. 93 + 5 = 98 , record the 8, carry 9). You can see immediately that the column sums in the right-hand method are componentsof the convolutionof the digit strings, for example 118 = 4 ×9+5× 8+6×7.I n§13.1 we learned how to compute the convolution of two vectors by the fast Fouriertransform(FFT): Each vectoris FFT’d, the two complextransforms are multiplied, and the result is inverse-FFT’d. Since the transforms are done with floating arithmetic, we need sufficient precision so that the exact integer value ofeach component of the result is discernible in the presence of roundoff error. We should therefore allow a (conservative) few times log 2(log2N)bits for roundoff in the FFT. A number of length Nbytes in radix 256 can generate convolution components as large as the order of (256)2N, thus requiring 16 + log2Nbits of precision for exact storage. If itis the number of bits in the floating mantissa (cf.§20.1), we obtain the condition 16 + log2N+few×log2log2N< it (20.6.3 ) We see that single precision, say with it =2 4, is inadequate for any interesting value of N, while double precision, say with it =5 3, allows Nto be greater than 106, corresponding to some millions of decimal digits. The following routine 910 Chapter20. Less-NumericalAlgorithmsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).thereforepresumesdoubleprecisionversionsof realft(§12.3)and four1(§12.2), here called drealft anddfour1. (These routines are included on the Numerical Recipesdiskettes.) SUBROUTINE mpmul(w,u,v,n,m) INTEGER m,n,NMAXCHARACTER*1 w(n+m),u(n),v(m)DOUBLE PRECISION RX PARAMETER (NMAX=8192,RX=256.D0) C USES drealft DOUBLE PRECISION version of realft. Uses Fast Fourier Transform to multiply the unsigned radix 256 integers u(1:n)and v(1:m), yielding a product w(1:n+m). INTEGER j,mn,nn DOUBLE PRECISION cy,t,a(NMAX),b(NMAX)mn=max(m,n) nn=1 Find the smallestuseable power oftwo for the transform. 1 if(nn.lt.mn) then nn=nn+nn goto 1 endif nn=nn+nnif(nn.gt.NMAX)pause ’NMAX too small in fftmul’ do 11j=1,n Move Uto a double precision floating array. a(j)=ichar(u(j)) enddo 11 do12j=n+1,nn a(j)=0.D0 enddo 12 do13j=1,m Move Vto a double precision floating array. b(j)=ichar(v(j)) enddo 13 do14j=m+1,nn b(j)=0.D0 enddo 14 Perform the convolution: First, the two Fourier transforms. call drealft(a,nn,1)call drealft(b,nn,1)b(1)=b(1)*a(1) Thenmultiplythecomplexresults(realandimaginaryparts). b(2)=b(2)*a(2) do 15j=3,nn,2 t=b(j) b(j)=t*a(j)-b(j+1)*a(j+1) b(j+1)=t*a(j+1)+b(j+1)*a(j) enddo 15 call drealft(b,nn,-1) Then do the inverse Fourier transform. cy=0. M a k eafi n a lp a s st od oa l lt h ec a r r i e s . do16j=nn,1,-1 t=b(j)/(nn/2)+cy+0.5D0 The 0.5 allows for roundoff error. b(j)=mod(t,RX) cy=int(t/RX) enddo 16 if (cy.ge.RX) pause ’cannot happen in fftmul’ w(1)=char(int(cy)) Copy answer to output. do17j=2,n+m w(j)=char(int(b(j-1))) enddo 17 return END With multiplication thus a “fast” operation, division is best performed by multiplying the dividend by the reciprocal of the divisor. The reciprocal of a value 20.6ArithmeticatArbitraryPrecision 911Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Vis calculated by iteration of Newton’s rule, Ui+1=Ui(2−VU i)( 20.6.4 ) which results in the quadratic convergence of U∞to1/V, as you can easily prove. (Many supercomputers and RISC machines actually use this iteration toperformdivisions.) We can now see where the operationscount NlogNlog log N, mentionedabove,originates: NlogNis in the Fouriertransform,with the iteration to converge Newton’s rule giving an additional factor of log log N. SUBROUTINE mpinv(u,v,n,m) INTEGER m,n,MF,NMAXCHARACTER*1 u(n),v(m) REAL BI PARAMETER (MF=4,BI=1./256.,NMAX=8192) Character string v(1:m)is interpreted as a radix 256 number with the radix point after (nonzero) v(1);u(1:n)issettothemostsignificantdigitsofitsreciprocal,withtheradix point after u(1). C USES mpmov,mpmul,mpneg INTEGER i,j,mmREAL fu,fv CHARACTER*1 rr(2*NMAX+1),s(NMAX)if(max(n,m).gt.NMAX)pause ’NMAX too small in mpinv’mm=min(MF,m) fv=ichar(v(mm)) Use ordinary floating arithmetic to get an initial ap- proximation. do 11j=mm-1,1,-1 fv=fv*BI+ichar(v(j)) enddo 11 fu=1./fv do12j=1,n i=int(fu) u(j)=char(i) fu=256.*(fu-i) enddo 12 1 continue Iterate Newton’s rule to convergence. call mpmul(rr,u,v,n,m) Construct 2−UVinS. call mpmov(s,rr(2),n)call mpneg(s,n) s(1)=char(ichar(s(1))-254) Multiply SUintoU. call mpmul(rr,s,u,n,n)call mpmov(u,rr(2),n)do 13j=2,n-1 Iffractionalpartof Sisnotzero,ithasnotconverged to1. if(ichar(s(j)).ne.0)goto 1 enddo 13 continue return END Divisionnowfollowsasasimplecorollary,withonlythenecessityofcalculating the reciprocal to sufficient accuracy to get an exact quotient and remainder. SUBROUTINE mpdiv(q,r,u,v,n,m) INTEGER m,n,NMAX,MACC CHARACTER*1 q(n-m+1),r(m),u(n),v(m) PARAMETER (NMAX=8192,MACC=6) Divides unsigned radix 256 integers u(1:n)byv(1:m)(with m≤nrequired), yielding a quotient q(1:n-m+1) and a remainder r(1:m). C USES mpinv,mpmov,mpmul,mpsad,mpsub INTEGER isCHARACTER*1 rr(2*NMAX),s(2*NMAX) if(n+MACC.gt.NMAX)pause ’NMAX too small in mpdiv’ 912 Chapter20. Less-NumericalAlgorithmsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).call mpinv(s,v,n+MACC,m) SetS=1/V. call mpmul(rr,s,u,n+MACC,n) SetQ=SU. call mpsad(s,rr,n+MACC-1,1) call mpmov(q,s(3),n-m+1)call mpmul(rr,q,v,n-m+1,m) Multiply and subtract to get the remainder. call mpsub(is,rr(2),u,rr(2),n) if (is.ne.0) pause ’MACC too small in mpdiv’call mpmov(r,rr(n-m+2),m)return END Square roots are calculated by a Newton’s rule much like division. If Ui+1=1 2Ui(3−VU2 i)( 20.6.5 ) thenU∞convergesquadraticallyto 1/√ V. A finalmultiplicationby Vgives√ V. SUBROUTINE mpsqrt(w,u,v,n,m) INTEGER m,n,NMAX,MF CHARACTER*1 w(*),u(*),v(*)REAL BIPARAMETER (NMAX=2048,MF=3,BI=1./256.) C USES mplsh,mpmov,mpmul,mpneg,mpsdv Character string v(1:m)is interpreted as a radix 256 number with the radix point after v(1);w(1:n)isset toitssquare root (radix point after w(1)), and u(1:n)isset tothe reciprocal thereof (radix point before u(1)).wanduneed not be distinct, in which case they are set to the square root. INTEGER i,ir,j,mmREAL fu,fv CHARACTER*1 r(NMAX),s(NMAX) if(2*n+1.gt.NMAX)pause ’NMAX too small in mpsqrt’mm=min(m,MF) fv=ichar(v(mm)) Use ordinary floating arithmetic to get an initial approx- imation. do 11j=mm-1,1,-1 fv=BI*fv+ichar(v(j)) enddo 11 fu=1./sqrt(fv)do 12j=1,n i=int(fu)u(j)=char(i) fu=256.*(fu-i) enddo 12 1 continue Iterate Newton’s rule to convergence. call mpmul(r,u,u,n,n) Construct S=( 3−VU2)/2. call mplsh(r,n)call mpmul(s,r,v,n,m)call mplsh(s,n) call mpneg(s,n) s(1)=char(ichar(s(1))-253)call mpsdv(s,s,n,2,ir) do 13j=2,n-1 If fractional part of Sis not zero, it has not converged to1. if(ichar(s(j)).ne.0)goto 2 enddo 13 call mpmul(r,u,v,n,m) Get square root from reciprocal and return. call mpmov(w,r(2),n) return 2 continue call mpmul(r,s,u,n,n) Replace UbySU. call mpmov(u,r(2),n) goto 1END 20.6ArithmeticatArbitraryPrecision 913Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).We already mentioned that radix conversion to decimal is a merely cosmetic operationthatshouldnormallybeomitted. Thesimplestwaytoconvertafractiontodecimalis tomultiplyit repeatedlyby 10,pickingoff(andsubtracting)theresulting integerpart. This, has an operationscount of O(N 2), however,since each liberated decimal digit takes an O(N)operation. It ispossible to do the radix conversion as a fast operation by a “divide and conquer” strategy, in which the fraction is (fast) multiplied by a large power of 10, enough to move about half the desired digits to the left of the radix point. The integer and fractional pieces are now processed independently, each further subdivided. If our goal were a few billion digits of π, instead of a few thousand, we would need to implement this scheme. For presentpurposes, the following lazy routine is adequate: SUBROUTINE mp2dfr(a,s,n,m) INTEGER m,n,IAZ CHARACTER*1 a(*),s(*) PARAMETER (IAZ=48) C USES mplsh,mpsmu Converts a radix 256 fraction a(1:n)(radix point before a(1)) to a decimal fraction representedasanasciistring s(1:m),where misareturnedvalue. Theinputarray a(1:n) isdestroyed. NOTE:Forsimplicity,thisroutineimplementsaslow( ∝N2)algorithm. Fast (∝NlnN), more complicated, radix conversion algorithms do exist. INTEGER j m=2.408*ndo 11j=1,m call mpsmu(a,a,n,10) s(j)=char(ichar(a(1))+IAZ)call mplsh(a,n) enddo 11 returnEND Finally,then,wearriveataroutineimplementingequations(20.6.1)and(20.6.2): SUBROUTINE mppi(n)INTEGER n,IAOFF,NMAX PARAMETER (IAOFF=48,NMAX=8192) C USES mpinit,mp2dfr,mpadd,mpinv,mplsh,mpmov,mpmul,mpsdv,mpsqrt Demonstrate multiple precision routines by calculating and printing the first nbytes of π. INTEGER ir,j,m CHARACTER*1 x(NMAX),y(NMAX),sx(NMAX),sxi(NMAX),t(NMAX),s(3*NMAX), * pi(NMAX) call mpinit t(1)=char(2) SetT=2. do11j=2,n t(j)=char(0) enddo 11 call mpsqrt(x,x,t,n,n) SetX0=√ 2. call mpadd(pi,t,x,n) Setπ0=2 +√ 2. call mplsh(pi,n) call mpsqrt(sx,sxi,x,n,n) SetY0=21/4. call mpmov(y,sx,n) 1 continue call mpadd(x,sx,sxi,n) SetXi+1 =(X1/2 i+X−1/2 i)/2. call mpsdv(x,x(2),n,2,ir) call mpsqrt(sx,sxi,x,n,n) Form the temporary T=YiX1/2 i+1+X−1/2 i+1. call mpmul(t,y,sx,n,n)call mpadd(t(2),t(2),sxi,n) x(1)=char(ichar(x(1))+1) Increment X i+1andYiby 1. 914 Chapter20. Less-NumericalAlgorithmsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).3.1415926535897932384626433832795028841971693993751058209749445923078164062 862089986280348253421170679821480865132823066470938446095505822317253594081 284811174502841027019385211055596446229489549303819644288109756659334461284 756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185 480744623799627495673518857527248912279381830119491298336733624406566430860 213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465 495853710507922796892589235420199561121290219608640344181598136297747713099 605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130 019278766111959092164201989380952572010654858632788659361533818279682303019 520353018529689957736225994138912497217752834791315155748572424541506959508295331168617278558890750983817546374649393192550604009277016711390098488240 128583616035637076601047101819429555961989467678374494482553797747268471040 475346462080466842590694912933136770289891521047521620569660240580381501935112533824300355876402474964732639141992726042699227967823547816360093417216412199245863150302861829745557067498385054945885869269956909272107975093029 553211653449872027559602364806654991198818347977535663698074265425278625518 184175746728909777727938000816470600161452491921732172147723501414419735685481613611573525521334757418494684385233239073941433345477624168625189835694 855620992192221842725502542568876717904946016534668049886272327917860857843 838279679766814541009538837863609506800642251252051173929848960841284886269456042419652850222106611863067442786220391949450471237137869609563643719172874677646575739624138908658326459958133904780275900994657640789512694683983 525957098258226205224894077267194782684826014769909026401363944374553050682 034962524517493996514314298091906592509372216964615157098583874105978859597729754989301617539284681382686838689427741559918559252459539594310499725246808459872736446958486538367362226260991246080512438843904512441365497627807 977156914359977001296160894416948685558484063534220722258284886481584560285 Figure 20.6.1. The first 2398 decimal digits of π, computed by the routines in this section. y(1)=char(ichar(y(1))+1) call mpinv(s,y,n,n) SetY i+1 =T/(Yi+1 ). call mpmul(y,t(3),s,n,n)call mplsh(y,n)call mpmul(t,x,s,n,n) Form temporary T=(X i+1 +1 )/(Yi+1 ). continue IfT=1t h e nweh a v ec o n v e r g e d . m=mod(255+ichar(t(2)),256)do 12j=3,n if(ichar(t(j)).ne.m)goto 2 enddo 12 if (abs(ichar(t(n+1))-m).gt.1)goto 2write (*,*) ’pi=’ s(1)=char(ichar(pi(1))+IAOFF) s(2)=’.’call mp2dfr(pi(2),s(3),n-1,m) Converttodecimalforprinting. NOTE:Theconversionroutine,forthisdemonstra- tion only, is a slow( ∝N 2) algorithm. Fast ( ∝NlnN), more complicated, radix conversion algorithms do exist. write (*,’(1x,64a1)’) (s(j),j=1,m+1) return 2 continue call mpmul(s,pi,t(2),n,n) Setπi+1 =Tπ i. call mpmov(pi,s(2),n) goto 1 END 20.6ArithmeticatArbitraryPrecision 915Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Figure 20.6.1 gives the result, computed with n= 1000. As an exercise, you might enjoy checking the first hundreddigits of the figure against the first 12 termsof Ramanujan’s celebrated identity [3] 1 π=√ 8 9801∞/summationdisplay n=0(4n)! (1103 + 26390 n) (n! 396n)4(20.6.6 ) using the above routines. You might also use the routines to verify that the number 2512+1is not a prime, but has factors 2,424,833 and 7,455,602,825,647,884,208,337,395,736,200,454,918,783,366,342,657 (which are in fact prime; the remaining prime factor being about 7.416×1098)[4]. CITED REFERENCES AND FURTHER READING: Borwein,J.M.,andBorwein,P.B.1987, PiandtheAGM:AStudyinAnalyticNumberTheoryand Computational Complexity (New York: Wiley). [1] Knuth,D.E.1981, SeminumericalAlgorithms ,2nded.,vol.2of TheArtofComputerProgramming (Reading, MA: Addison-Wesley), §4.3. [2] Ramanujan, S. 1927, Collected Papers of Srinivasa Ramanujan , G.H. Hardy, P.V. Seshu Aiyar, and B.M. Wilson,eds. (Cambridge, U.K.: Cambridge University Press), pp. 23–39. [3] Kolata, G. 1990, June 20, The New York Times . [4] Kronsj¨o, L. 1987, Algorithms: Their Complexity and Efficiency , 2nd ed. (New York: Wiley).