c20-6 (779628), страница 2

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 2)

Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Full multiplication of two digit strings, if done by the traditional hand method,is not a fast operation: In multiplying two strings of length N , the multiplicandwould be short-multiplied in turn by each byte of the multiplier, requiring O(N 2 )operations in all. We will see, however, that all the arithmetic operations on numbersof length N can in fact be done in O(N × log N × log log N ) operations.The trick is to recognize that multiplication is essentially a convolution (§13.1)of the digits of the multiplicand and multiplier, followed by some kind of carryoperation.

Consider, for example, two ways of writing the calculation 456 × 789:20.6 Arithmetic at Arbitrary Precision919#include "nrutil.h"#define RX 256.0mn=IMAX(m,n);while (nn < mn) nn <<= 1;Find the smallest usable power of two for the transnn <<= 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 transdrealft(b,nn,1);forms.b[1] *= a[1];Then multiply the complex results (real and imagib[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 bymultiplying the dividend by the reciprocal of the divisor.

The reciprocal of a valueV 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 easilyprove. (Many supercomputers and RISC machines actually use this iteration toperform 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 iterationto converge Newton’s rule giving an additional factor of log log N .Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.

To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).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;920Chapter 20.Less-Numerical Algorithms#include "nrutil.h"#define MF 4#define BI (1.0/256)maxmn=IMAX(n,m);rr=cvector(1,1+(maxmn<<1));s=cvector(1,maxmn);mm=IMIN(MF,m);fv=(float) v[mm];Use ordinary floating arithmetic to get an initial apfor (j=mm-1;j>=1;j--) {proximation.fv *= BI;fv += v[j];}fu=1.0/fv;for (j=1;j<=n;j++) {i=(int) fu;u[j]=(unsigned char) i;fu=256.0*(fu-i);}for (;;) {Iterate Newton’s rule to convergence.mpmul(rr,u,v,n,m);Construct 2 − U V in S.mpmov(s,&rr[1],n);mpneg(s,n);s[1] -= 254;Multiply SU into U .mpmul(rr,s,u,n,n);mpmov(u,&rr[1],n);for (j=2;j<n;j++)If fractional part of S is not zero, it has not convergedif (s[j]) break;to 1.if (j==n) {free_cvector(s,1,maxmn);free_cvector(rr,1,1+(maxmn<<1));return;}}}Division now follows as a simple corollary, with only the necessity of calculatingthe reciprocal to sufficient accuracy to get an exact quotient and remainder.#include "nrutil.h"#define MACC 6void mpdiv(unsigned char q[], unsigned char r[], unsigned char u[],unsigned char v[], int n, int m)Divides unsigned radix 256 integers u[1..n] by v[1..m] (with m ≤ n required), yielding aquotient q[1..n-m+1] and a remainder r[1..m].{void mpinv(unsigned char u[], unsigned char v[], int n, int m);void mpmov(unsigned char u[], unsigned char v[], int n);Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).void mpinv(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(nonzero) v[1]; u[1..n] is set to the most significant digits of its reciprocal, with the radixpoint after u[1].{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);unsigned char *rr,*s;int i,j,maxmn,mm;float fu,fv;20.6 Arithmetic at Arbitrary Precision921void 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;Set S = 1/V .Set Q = SU .Multiply and subtract to get the remainder.mpdiv");}Square roots are calculated by a Newton’s rule much like division.

If1Ui (3 − V Ui2 )(20.6.5)2√√converges quadratically to 1/ V . A final multiplication by V gives V .Ui+1 =then U∞#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 reciprocalthereof (radix point before u[1]). w and u need not be distinct, in which case they are setto 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];for (j=mm-1;j>=1;j--) {fv *= BI;fv += v[j];}fu=1.0/sqrt(fv);for (j=1;j<=n;j++) {i=(int) fu;u[j]=(unsigned char) i;Use ordinary floating arithmetic to get an initial approximation.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).rr=cvector(1,(n+MACC)<<1);s=cvector(1,(n+MACC)<<1);mpinv(s,v,n+MACC,m);mpmul(rr,s,u,n+MACC,n);mpsad(s,rr,n+n+MACC/2,1);mpmov(q,&s[2],n-m+1);mpmul(rr,q,v,n-m+1,m);mpsub(&is,&rr[1],u,&rr[1],n);if (is) nrerror("MACC too small inmpmov(r,&rr[n-m+1],m);free_cvector(s,1,(n+MACC)<<1);free_cvector(rr,1,(n+MACC)<<1);922Chapter 20.Iterate Newton’s rule to convergence.Construct S = (3 − V U 2 )/2.If fractional part of S is not zero, it has not convergedto 1.Replace U by SU .Get square root from reciprocal and return.}We already mentioned that radix conversion to decimal is a merely cosmeticoperation that should normally be omitted.

Характеристики

Список файлов книги