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

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

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

The simplest way to convert a fraction todecimal is to multiply it repeatedly by 10, picking off (and subtracting) the resultinginteger part. This, has an operations count of O(N 2 ), however, since each liberateddecimal digit takes an O(N ) operation. It is possible to do the radix conversion asa 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 digitsto the left of the radix point.

The integer and fractional pieces are now processedindependently, 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 presentpurposes, the following lazy routine is adequate:#define IAZ 48void 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 represented 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 (∝ N 2) 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):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).fu=256.0*(fu-i);}for (;;) {mpmul(r,u,u,n,n);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 (s[j]) {mpmul(r,s,u,n,n);mpmov(u,&r[1],n);break;}}if (j<n) continue;mpmul(r,u,v,n,IMIN(m,n));mpmov(w,&r[1],n);free_cvector(s,1,n<<1);free_cvector(r,1,n<<1);return;}Less-Numerical Algorithms20.6 Arithmetic at Arbitrary Precision923Figure 20.6.1. The first 2398 decimal digits of π, computed by the routines in this section.#include <stdio.h>#include "nrutil.h"#define IAOFF 48void mppi(int n)Demonstrate multiple precision routines by calculating and printing the first n bytes of π.{void mp2dfr(unsigned char a[], unsigned char s[], int n, int *m);void mpadd(unsigned char w[], unsigned char u[], unsigned char v[], int n);void mpinv(unsigned char u[], unsigned char v[], int n, int m);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 mpsdv(unsigned char w[], unsigned char u[], int n, int iv, int *ir);void mpsqrt(unsigned char w[], unsigned char u[], unsigned char v[], int n,int m);int ir,j,m;unsigned char mm,*x,*y,*sx,*sxi,*t,*s,*pi;x=cvector(1,n+1);y=cvector(1,n<<1);sx=cvector(1,n);sxi=cvector(1,n);t=cvector(1,n<<1);s=cvector(1,3*n);pi=cvector(1,n+1);t[1]=2;for (j=2;j<=n;j++) t[j]=0;Set T = 2.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).3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989380952572010654858632788659361533818279682303019520353018529689957736225994138912497217752834791315155748572424541506959508295331168617278558890750983817546374649393192550604009277016711390098488240128583616035637076601047101819429555961989467678374494482553797747268471040475346462080466842590694912933136770289891521047521620569660240580381501935112533824300355876402474964732639141992726042699227967823547816360093417216412199245863150302861829745557067498385054945885869269956909272107975093029553211653449872027559602364806654991198818347977535663698074265425278625518184175746728909777727938000816470600161452491921732172147723501414419735685481613611573525521334757418494684385233239073941433345477624168625189835694855620992192221842725502542568876717904946016534668049886272327917860857843838279679766814541009538837863609506800642251252051173929848960841284886269456042419652850222106611863067442786220391949450471237137869609563643719172874677646575739624138908658326459958133904780275900994657640789512694683983525957098258226205224894077267194782684826014769909026401363944374553050682034962524517493996514314298091906592509372216964615157098583874105978859597729754989301617539284681382686838689427741559918559252459539594310499725246808459872736446958486538367362226260991246080512438843904512441365497627807977156914359977001296160894416948685558484063534220722258284886481584560285924Chapter 20.Less-Numerical Algorithms}Figure 20.6.1 gives the result, computed with n = 1000.

As an exercise, youmight enjoy checking the first hundred digits of the figure against the first 12 termsof Ramanujan’s celebrated identity [3]√ ∞18 X (4n)! (1103 + 26390n)=π9801 n=0(n! 396n)4(20.6.6)using the above routines. You might also use the routines to verify that thenumber 2512 + 1 is not a prime, but has factors 2,424,833 and7,455,602,825,647,884,208,337,395,736,200,454,918,783,366,342,657 (which arein fact prime; the remaining prime factor being about 7.416 × 1098 ) [4].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).√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 (;;) {1/2−1/2mpadd(x,sx,sxi,n);Set Xi+1 = (Xi + Xi)/2.mpsdv(x,&x[1],n,2,&ir);1/2−1/2mpsqrt(sx,sxi,x,n,n);Form the temporary T = Yi Xi+1 + Xi+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.

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

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