924 Chapter 20.Less-Numerical Algorithms mpsqrt(x,x,t,n,n); Set Xo v2. mpadd(pi,t,x,n); Set To =2+v2 mplsh(pi,n); mpsqrt(sx,sxi,x,n,n); Set Yo =21/4 mpmov(y,sx,n); for (;;) mpadd(x,sx,sxi,n); SetX+1=(X/2+X-1/2)/2 mpsdv(x,&x[1],n,2,&ir); mpsqrt(sx,sxi,x,n,n); Form the temporaryT=Y mpmul(t,y,sx,n,n); mpadd(&t[1],&t[1],sxi,n); x[1]++; Increment Xi+1 and Yi by 1. y[1]+; mpinv(s,y,n,n); SetY+1=T/(Y:+1) mpmul(y,&t [2],s,n,n); 19881992 mplsh(y,n); mpmul(t,x,s,n,n); Form temporary T=(Xi+1+1)/(Yi+1). mm=t[2]-1; If T =1 then we have converged. 100 for(j=3;j<=n;j++)[ if (t [j]!mm)break; Cambridge from NUMERICAL RECIPESI m=t [n+1]-mm 1f(j<=n11m>1Ilm<-1){ mpmul(s,pi,&t [1],n,n); Set不i+1=TT mpmov (pi,&s [1],n); continue; server computer, (North America University Press. THE 2 printf("pi=\n"); ART s[1]=p1[1]+IA0FF; s[2]=’.’: Programs m=mm: strictly proh mp2dfr(pi[1],&s[2],n-1,&m); Convert to decimal for printing.NOTE:The conversion routine,for this demonstration only,is a slow N2)algorithm.Fast (o N In N),more complicated,radix conversion algorithms do exist. s[m+3]=0; to dir printf("%64s\n",&s[1]); free_cvector(pi,1,n+1); free_cvector(s,1,3*n); ectcustser 18881920 OF SCIENTIFIC COMPUTING(ISBN free_cvector(t,1,n<<1); free_cvector(sxi,1,n); free_cvector(sx,1,n); free_cvector(y,1,n<<1); Numerical Recipes 10-621 free_cvector(x,1,n+1); return; 43108 (outside Figure 20.6.1 gives the result,computed with n =1000.As an exercise,you North Software. might enjoy checking the first hundred digits of the figure against the first 12 terms of Ramanujan's celebrated identity [3] Ame ⑧、 visit website 1 (4n)1(1103+26390m) machine (20.6.6) 9801 (n!396n)4 n=0 using the above routines.You might also use the routines to verify that the number 2512+1 is not a prime,but has factors 2,424,833 and 7.455.602,825,647.884.208.337,395.736,200.454.918.783366,342.657(which are in fact prime;the remaining prime factor being about 7.416 x 1098)[4].924 Chapter 20. Less-Numerical Algorithms 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). mpsqrt(x,x,t,n,n); Set X0 = √2. mpadd(pi,t,x,n); Set π0 = 2+ √2. mplsh(pi,n); mpsqrt(sx,sxi,x,n,n); Set Y0 = 21/4. mpmov(y,sx,n); for (;;) { mpadd(x,sx,sxi,n); Set Xi+1 = (X1/2 i + X−1/2 i )/2. mpsdv(x,&x[1],n,2,&ir); mpsqrt(sx,sxi,x,n,n); Form the temporary T = YiX1/2 i+1 + X−1/2 i+1 . mpmul(t,y,sx,n,n); mpadd(&t[1],&t[1],sxi,n); x[1]++; Increment Xi+1 and Yi by 1. y[1]++; mpinv(s,y,n,n); Set Yi+1 = T /(Yi + 1). mpmul(y,&t[2],s,n,n); mplsh(y,n); mpmul(t,x,s,n,n); Form temporary T = (Xi+1 + 1)/(Yi + 1). mm=t[2]-1; If T = 1 then we have converged. for (j=3;j<=n;j++) { if (t[j] != mm) break; } m=t[n+1]-mm; if (j <= n || m > 1 || m < -1) { mpmul(s,pi,&t[1],n,n); Set πi+1 = T πi. mpmov(pi,&s[1],n); continue; } printf("pi=\n"); s[1]=pi[1]+IAOFF; s[2]=’.’; m=mm; mp2dfr(&pi[1],&s[2],n-1,&m); Convert to decimal for printing. NOTE: The conversion routine, for this demonstration only, is a slow (∝ N2) algorithm. Fast (∝ N ln N), more complicated, radix conversion algorithms do exist. s[m+3]=0; printf(" %64s\n",&s[1]); free_cvector(pi,1,n+1); free_cvector(s,1,3*n); free_cvector(t,1,n<<1); free_cvector(sxi,1,n); free_cvector(sx,1,n); free_cvector(y,1,n<<1); free_cvector(x,1,n+1); return; } } Figure 20.6.1 gives the result, computed with n = 1000. As an exercise, you might enjoy checking the first hundred digits of the figure against the first 12 terms of Ramanujan’s celebrated identity [3] 1 π = √8 9801 ∞ n=0 (4n)! (1103 + 26390n) (n! 396n)4 (20.6.6) using the above routines. You might also use the routines to verify that the number 2512 + 1 is 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 × 10 98) [4]