20.6 Arithmetic at Arbitrary Precision 921 void mpmul(unsigned char w[],unsigned char u[],unsigned char v[],int n, int m); void mpsad(unsigned char v[],unsigned char u[],int n,int iv); void mpsub(int *is,unsigned char w[],unsigned char u[],unsigned char v[], int n); int is; unsigned char*rr,米s; rr=cvector(1,(n+MACC)<<1); s=cvector(1,(n+MACC)<<1); mpinv(s,v,n+MACC,m); Set S=1/V mpmul(rr,s,u,n+MACC,n); Set Q=SU. mpsad(s,rr,n+MACC-1,1) mpmov (q,&s [2],n-m+1); mpmul(rr,q,v,n-m+1,m); Multiply and subtract to get the remainder mpsub(&is,&rr[1],u,&rr[1],n); granted for 18881892 if (is)nrerror("MACC too small in mpdiv"); mpmov(r,&rr[n-m+1],m); 1800 free_cvector(s,1,(n+MACC)<<1): free_cvector(rr,1,(n+MACC)<<1); from NUMERICAL RECIPES I Square roots are calculated by a Newton's rule much like division.If server (Nort U+1=2,3-VW） (20.6.5) America computer Press. ART then U converges quadratically to 1/VV.A final multiplication by V gives vV. 9 #include <math.h> Programs #include "nrutil.h" #define MF 3 #define BI (1.0/256) void mpsqrt(unsigned char w],unsigned char un,unsigned char v,int n, to dir int m) Character string v[1..m]is interpreted as a radix 256 number with the radix point after v[1]: w[1..n]is set to its square root (radix point after w[1]),and u[1..n]is set to the reciprocal 1920 SCIENTIFIC COMPUTING(ISBN thereof (radix point before u[1]).w and u need not be distinct,in which case they are set to the square root. void mplsh(unsigned char u[▣，intn); void mpmov(unsigned char u[],unsigned char v,int n); void mpmul(unsigned char w[],unsigned char u[],unsigned char v[],int n, idge.org Numerical Recipes 0621 -43108 int m); void mpneg(unsigned char u[],int n); void mpsdv(unsigned char w[],unsigned char u],int n,int iv,int *ir); int i,ir,i,mm; (outside float fu,fv; North Software. unsigned char r,*s; r=cvector(1,n<<1); s=cvector(1,n<<1); mm=IMIN(m,MF) fv=(float)v[mm]; Use ordinary floating arithmetic to get an initial ap- or（灯=mm-1j>=1;j-)( proximation fV米=BI; fv +v[j]; fu=1.0/sqrt(fv); for(j=1;j<=n:j++){ i=(int)fu; u[j]=(unsigned char)i;20.6 Arithmetic at Arbitrary Precision 921 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). void mpmul(unsigned char w[], unsigned char u[], unsigned char v[], int n, int m); void mpsad(unsigned char w[], unsigned char u[], int n, int iv); void mpsub(int *is, unsigned char w[], unsigned char u[], unsigned char v[], int n); int is; unsigned char *rr,*s; rr=cvector(1,(n+MACC)<<1); s=cvector(1,(n+MACC)<<1); mpinv(s,v,n+MACC,m); Set S = 1/V . mpmul(rr,s,u,n+MACC,n); Set Q = SU. mpsad(s,rr,n+MACC-1,1); mpmov(q,&s[2],n-m+1); mpmul(rr,q,v,n-m+1,m); Multiply and subtract to get the remainder. mpsub(&is,&rr[1],u,&rr[1],n); if (is) nrerror("MACC too small in mpdiv"); mpmov(r,&rr[n-m+1],m); free_cvector(s,1,(n+MACC)<<1); free_cvector(rr,1,(n+MACC)<<1); } Square roots are calculated by a Newton’s rule much like division. If Ui+1 = 1 2 Ui(3 − V U2 i ) (20.6.5) then U∞ converges quadratically to 1/ √ V . A final multiplication by V gives √ V . #include <math.h> #include "nrutil.h" #define MF 3 #define BI (1.0/256) void mpsqrt(unsigned char w[], unsigned char u[], unsigned char v[], int n, int m) Character string v[1..m] is interpreted as a radix 256 number with the radix point after v[1]; w[1..n] is set to its square root (radix point after w[1]), and u[1..n] is set to the reciprocal thereof (radix point before u[1]). w and u need not be distinct, in which case they are set to the square root. { void mplsh(unsigned char u[], int n); void mpmov(unsigned char u[], unsigned char v[], int n); void mpmul(unsigned char w[], unsigned char u[], unsigned char v[], int n, int m); void mpneg(unsigned char u[], int n); void mpsdv(unsigned char w[], unsigned char u[], int n, int iv, int *ir); int i,ir,j,mm; float fu,fv; unsigned char *r,*s; r=cvector(1,n<<1); s=cvector(1,n<<1); mm=IMIN(m,MF); fv=(float) v[mm]; Use ordinary floating arithmetic to get an initial ap￾for (j=mm-1;j>=1;j--) { proximation. fv *= BI; fv += v[j]; } fu=1.0/sqrt(fv); for (j=1;j<=n;j++) { i=(int) fu; u[j]=(unsigned char) i;