922 Chapter 20. Less-Numerical Algorithms fu=256.0*(fu-1); for(;;)[ Iterate Newton's rule to convergence. mpmul(r,u,u,n,n) Construct S=(3-VU2)/2. mplsh(r,n); mpmul(s,r,v,n,IMIN(m,n)) mplsh(s,n); mpneg(s,n); s[1]-=253; mpsdv(s,s,n,2,&ir); for (j=2;j<nij++) If fractional part of S is not zero,it has not converged if (s[j]) to l mpmul(r,s,u,n,n); Replace U by SU. mpmov(u,&r[1],n); break; 83g granted for 19881992 if (j<n)continue; 1.800 mpmul(r,u,v,n,IMIN(m,n)); Get square root from reciprocal and return mpmov(w,&r[1],n); free_cvector(s,1,n<<1); free_cvector(r,1,n<<1); return; (Nort We already mentioned that radix conversion to decimal is a merely cosmetic America server computer, operation that should normally be omitted.The simplest way to convert a fraction to ART decimal is to multiply it repeatedly by 10,picking off(and subtracting)the resulting integer part.This,has an operations count of O(N2),however,since each liberated Programs decimal digit takes an O(N)operation.It is possible 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 to dir independently,each further subdivided.If our goal were a few billion digits of m, instead of a few thousand,we would need to implement this scheme.For present OF SCIENTIFIC COMPUTING(ISBN purposes,the following lazy routine is adequate: 1988-19920 #define IAZ 48 10-621 void mp2dfr(unsigned char a[],unsigned char s[],int n,int *m) Numerical Recipes 43108 Converts a radix 256 fraction a[1..n](radix point before a[1])to a decimal fraction rep- resented as an ascii string s[1..m],where m is a returned value.The input array a[1..n] is destroyed.NOTE:For simplicity,this routine implements a slow N2)algorithm.Fast NIn N),more complicated,radix conversion algorithms do exist. (outside Software. void mplsh(unsigned char u[],int n); void mpsmu(unsigned char w[],unsigned char u[],int n,int iv); int ji ying of *m=(int)(2.408n): for(j=1;j<=(*m);j+)[ mpsmu(a,a,n,10); s[j]=a[1]+IAZ; mplsh(a,n); Finally,then,we arrive at a routine implementing equations(20.6.1)and(20.6.2):922 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). fu=256.0*(fu-i); } for (;;) { Iterate Newton’s rule to convergence. mpmul(r,u,u,n,n); Construct S = (3 − V U2)/2. mplsh(r,n); mpmul(s,r,v,n,IMIN(m,n)); mplsh(s,n); mpneg(s,n); s[1] -= 253; mpsdv(s,s,n,2,&ir); for (j=2;j<n;j++) { If fractional part of S is not zero, it has not converged if (s[j]) { to 1. mpmul(r,s,u,n,n); Replace U by SU. mpmov(u,&r[1],n); break; } } if (j<n) continue; mpmul(r,u,v,n,IMIN(m,n)); Get square root from reciprocal and return. mpmov(w,&r[1],n); free_cvector(s,1,n<<1); free_cvector(r,1,n<<1); return; } } We already mentioned that radix conversion to decimal is a merely cosmetic operation that should normally be omitted. The simplest way to convert a fraction to decimal is to multiply it repeatedly by 10, picking off (and subtracting) the resulting integer part. This, has an operations count of O(N 2), however, since each liberated decimal digit takes an O(N) operation. It is possible 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 present purposes, the following lazy routine is adequate: #define IAZ 48 void mp2dfr(unsigned char a[], unsigned char s[], int n, int *m) Converts a radix 256 fraction a[1..n] (radix point before a[1]) to a decimal fraction rep￾resented as an ascii string s[1..m], where m is a returned value. The input array a[1..n] is destroyed. NOTE: For simplicity, this routine implements a slow (∝ N2) algorithm. Fast (∝ N ln N), more complicated, radix conversion algorithms do exist. { void mplsh(unsigned char u[], int n); void mpsmu(unsigned char w[], unsigned char u[], int n, int iv); int j; *m=(int) (2.408*n); for (j=1;j<=(*m);j++) { mpsmu(a,a,n,10); s[j]=a[1]+IAZ; mplsh(a,n); } } Finally, then, we arrive at a routine implementing equations (20.6.1) and (20.6.2):