20.6 Arithmetic at Arbitrary Precision 919 #include "nrutil.h" #define RX 256.0 void mpmul(unsigned char w[],unsigned char u[],unsigned char v[],int n, int m) 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]. void drealft(double data[],unsigned long n,int isign); double version of realft. int j,mn,nn=1; double cy,t,*a,*b; mn=IMAX (m,n); while (nn mn)nn <<1; Find the smallest usable power of two for the trans- nn<<=1; form. a=dvector(1,nn); 19881992 b=dvector(1,nn); for (j=1;j<=n;j++) Move U to a double precision floating array. a[j]=(double)u[j]; for(j=n+1;j<nn;j++)a[j]=0.0; for （j=1;j<=m;j+） Move V to a double precision floating array. from NUMERICAL RECIPES I b[j]=(double)v[]; for(j=m+1;j<=nn;j++)b[j]=0.0; drealft(a,nn,1); Perform the convolution:First,the two Fourier trans- (Nort server drealft(b,nn,1); forms. b[1]*=a[1]: Then multiply the complex results (real and imagi- b[2]*=a[2]; nary parts). Americ computer, for(j=3;j<=nn;j+=2){ ART b[j]=(t=b[j])*a[j]-b[j+1]*a[j+1]; b[j+1]=t*a[j+1]+b[j+1]*a[j]; Programs drealft(b,nn,-1); Then do the inverse Fourier transform. cy=0.0; Make a final pass to do all the carries. for (i=nn;i>=1;--){ 1CIYP ic t=b[]/(nn>>1)+cy+0.5; The 0.5 allows for roundoff error. cy=(unsigned long)(t/RX); to dir b[j]=t-cy*RX; if (cy >=RX)nrerror("cannot happen in fftmul"); OF SCIENTIFIC COMPUTING(ISBN w[1]=(unsigned char)cy; Copy answer to output for (j=2;j<=n+m;j++) w[j]=(unsigned char)b[j-1]; free_dvector(b,1,nn); v@cambri 1988-1992 by Numerical Recipes 10-:6211 free_dvector(a,1,nn); 43108 With multiplication thus a "fast"operation,division is best performed by (outside multiplying the dividend by the reciprocal of the divisor.The reciprocal of a value V is calculated by iteration of Newton's rule, North Software. U:+1=U(2-VU) (20.6.4) visit website machine which results in the quadratic convergence of 0 to 1/V,as you can easily prove.(Many supercomputers and RISC machines actually use this iteration to perform divisions.)We can now see where the operations count Nlog N log log N. mentioned above,originates:Nlog N is in the Fourier transform,with the iteration to converge Newton's rule giving an additional factor of log log N.20.6 Arithmetic at Arbitrary Precision 919 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). #include "nrutil.h" #define RX 256.0 void mpmul(unsigned char w[], unsigned char u[], unsigned char v[], int n, int m) 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]. { void drealft(double data[], unsigned long n, int isign); double version of realft. int j,mn,nn=1; double cy,t,*a,*b; mn=IMAX(m,n); while (nn < mn) nn <<= 1; Find the smallest usable power of two for the trans￾nn <<= 1; form. a=dvector(1,nn); b=dvector(1,nn); for (j=1;j<=n;j++) Move U to a double precision floating array. a[j]=(double)u[j]; for (j=n+1;j<=nn;j++) a[j]=0.0; for (j=1;j<=m;j++) Move V to a double precision floating array. b[j]=(double)v[j]; for (j=m+1;j<=nn;j++) b[j]=0.0; drealft(a,nn,1); Perform the convolution: First, the two Fourier trans￾drealft(b,nn,1); forms. b[1] *= a[1]; Then multiply the complex results (real and imagi￾b[2] *= a[2]; nary parts). for (j=3;j<=nn;j+=2) { b[j]=(t=b[j])*a[j]-b[j+1]*a[j+1]; b[j+1]=t*a[j+1]+b[j+1]*a[j]; } drealft(b,nn,-1); Then do the inverse Fourier transform. cy=0.0; Make a final pass to do all the carries. for (j=nn;j>=1;j--) { t=b[j]/(nn>>1)+cy+0.5; The 0.5 allows for roundoff error. cy=(unsigned long) (t/RX); b[j]=t-cy*RX; } if (cy >= RX) nrerror("cannot happen in fftmul"); w[1]=(unsigned char) cy; Copy answer to output. for (j=2;j<=n+m;j++) w[j]=(unsigned char) b[j-1]; free_dvector(b,1,nn); free_dvector(a,1,nn); } 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 V is calculated by iteration of Newton’s rule, Ui+1 = Ui(2 − V Ui) (20.6.4) which results in the quadratic convergence of U∞ to 1/V , as you can easily prove. (Many supercomputers and RISC machines actually use this iteration to perform divisions.) We can now see where the operations count N log N log log N, mentioned above, originates: N log N is in the Fourier transform, with the iteration to converge Newton’s rule giving an additional factor of log log N