c20-6 (779628)
Текст из файла
20.6 Arithmetic at Arbitrary Precision915CITED REFERENCES AND FURTHER READING:Bell, T.C., Cleary, J.G., and Witten, I.H. 1990, Text Compression (Englewood Cliffs, NJ: PrenticeHall).Nelson, M. 1991, The Data Compression Book (Redwood City, CA: M&T Books).Witten, I.H., Neal, R.M., and Cleary, J.G. 1987, Communications of the ACM, vol. 30, pp. 520–540. [1]Let’s compute the number π to a couple of thousand decimal places.
In doingso, we’ll learn some things about multiple precision arithmetic on computers andmeet quite an unusual application of the fast Fourier transform (FFT). We’ll alsodevelop a set of routines that you can use for other calculations at any desired levelof arithmetic precision.To start with, we need an analytic algorithm for π. Useful algorithms arequadratically convergent, i.e., they double the number of significant digits ateach iteration. Quadratically convergent algorithms for π are based on the AGM(arithmetic geometric mean) method, which also finds application to the calculationof elliptic integrals (cf.
§6.11) and in advanced implementations of the ADI methodfor elliptic partial differential equations (§19.5). Borwein and Borwein [1] treat thissubject, which is beyond our scope here. One of their algorithms for π starts withthe initializations√X0 = 2√π0 = 2 + 2(20.6.1)√4Y0 = 2and then, for i = 0, 1, . . . , repeats the iterationXi+1πi+1Yi+1p1Xi + √XiXi+1 + 1= πiYi + 1pYi Xi+1 + p 1Xi+1=Yi + 11=2(20.6.2)The value π emerges as the limit π∞ .Now, to the question of how to do arithmetic to arbitrary precision: In ahigh-level language like C, a natural choice is to work in radix (base) 256, so thatcharacter arrays can be directly interpreted as strings of digits.
At the very end ofour calculation, we will want to convert our answer to radix 10, but that is essentiallya frill for the benefit of human ears, accustomed to the familiar chant, “three pointSample 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).20.6 Arithmetic at Arbitrary Precision916Chapter 20.Less-Numerical Algorithms#define LOBYTE(x) ((unsigned char) ((x) & 0xff))#define HIBYTE(x) ((unsigned char) ((x) >> 8 & 0xff))Multiple precision arithmetic operations done on character strings, interpreted as radix 256numbers.
This set of routines collects the simpler operations.void mpadd(unsigned char w[], unsigned char u[], unsigned char v[], int n)Adds the unsigned radix 256 integers u[1..n] and v[1..n] yielding the unsigned integerw[1..n+1].{int j;unsigned short ireg=0;for (j=n;j>=1;j--) {ireg=u[j]+v[j]+HIBYTE(ireg);w[j+1]=LOBYTE(ireg);}w[1]=HIBYTE(ireg);}void mpsub(int *is, unsigned char w[], unsigned char u[], unsigned char v[],int n)Subtracts the unsigned radix 256 integer v[1..n] from u[1..n] yielding the unsigned integerw[1..n].
If the result is negative (wraps around), is is returned as −1; otherwise it is returnedas 0.{int j;unsigned short ireg=256;for (j=n;j>=1;j--) {ireg=255+u[j]-v[j]+HIBYTE(ireg);w[j]=LOBYTE(ireg);}*is=HIBYTE(ireg)-1;}void mpsad(unsigned char w[], unsigned char u[], int n, int iv)Short addition: the integer iv (in the range 0 ≤ iv ≤ 255) is added to the unsigned radix 256integer u[1..n], yielding w[1..n+1].{int j;unsigned short ireg;ireg=256*iv;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).one four one five nine.
. . .” For any less frivolous calculation, we would likely neverleave base 256 (or the thence trivially reachable hexadecimal, octal, or binary bases).We will adopt the convention of storing digit strings in the “human” ordering,that is, with the first stored digit in an array being most significant, the last stored digitbeing least significant. The opposite convention would, of course, also be possible.“Carries,” where we need to partition a number larger than 255 into a low-orderbyte and a high-order carry, present a minor programming annoyance, solved, in theroutines below, by the use of the macros LOBYTE and HIBYTE.It is easy at this point, following Knuth [2], to write a routine for the “fast”arithmetic operations: short addition (adding a single byte to a string), addition,subtraction, short multiplication (multiplying a string by a single byte), shortdivision, ones-complement negation; and a couple of utility operations, copyingand left-shifting strings.
(On the diskette, these functions are all in the singlefile mpops.c.)20.6 Arithmetic at Arbitrary Precision917for (j=n;j>=1;j--) {ireg=u[j]+HIBYTE(ireg);w[j+1]=LOBYTE(ireg);}w[1]=HIBYTE(ireg);}for (j=n;j>=1;j--) {ireg=u[j]*iv+HIBYTE(ireg);w[j+1]=LOBYTE(ireg);}w[1]=HIBYTE(ireg);}void mpsdv(unsigned char w[], unsigned char u[], int n, int iv, int *ir)Short division: the unsigned radix 256 integer u[1..n] is divided by the integer iv (in therange 0 ≤ iv ≤ 255), yielding a quotient w[1..n] and a remainder ir (with 0 ≤ ir ≤ 255).{int i,j;*ir=0;for (j=1;j<=n;j++) {i=256*(*ir)+u[j];w[j]=(unsigned char) (i/iv);*ir=i % iv;}}void mpneg(unsigned char u[], int n)Ones-complement negate the unsigned radix 256 integer u[1..n].{int j;unsigned short ireg=256;for (j=n;j>=1;j--) {ireg=255-u[j]+HIBYTE(ireg);u[j]=LOBYTE(ireg);}}void mpmov(unsigned char u[], unsigned char v[], int n)Move v[1..n] onto u[1..n].{int j;for (j=1;j<=n;j++) u[j]=v[j];}void mplsh(unsigned char u[], int n)Left shift u(2..n+1) onto u[1..n].{int j;for (j=1;j<=n;j++) u[j]=u[j+1];}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 mpsmu(unsigned char w[], unsigned char u[], int n, int iv)Short multiplication: the unsigned radix 256 integer u[1..n] is multiplied by the integer iv(in the range 0 ≤ iv ≤ 255), yielding w[1..n+1].{int j;unsigned short ireg=0;918Chapter 20.Less-Numerical Algorithms456× 7894104364831923597844× 73632 4028 35 4228 67 1183 5 9 75 68 945 544893 548 4The tableau on the left shows the conventional method of multiplication, in whichthree separate short multiplications of the full multiplicand (by 9, 8, and 7) areadded to obtain the final result.
The tableau on the right shows a different method(sometimes taught for mental arithmetic), where the single-digit cross products areall computed (e.g. 8 × 6 = 48), then added in columns to obtain an incompletelycarried result (here, the list 28, 67, 118, 93, 54).
The final step is a single pass fromright to left, recording the single least-significant digit and carrying the higher digitor digits into the total to the left (e.g. 93 + 5 = 98, record the 8, carry 9).You can see immediately that the column sums in the right-hand method arecomponents of the convolution of the digit strings, for example 118 = 4 × 9 + 5 ×8 + 6 × 7. In §13.1 we learned how to compute the convolution of two vectors bythe fast Fourier transform (FFT): Each vector is FFT’d, the two complex transformsare multiplied, and the result is inverse-FFT’d.
Since the transforms are done withfloating arithmetic, we need sufficient precision so that the exact integer value ofeach component of the result is discernible in the presence of roundoff error. Weshould therefore allow a (conservative) few times log2 (log2 N ) bits for roundoffin the FFT. A number of length N bytes in radix 256 can generate convolutioncomponents as large as the order of (256)2 N , thus requiring 16 + log2 N bits ofprecision for exact storage. If it is the number of bits in the floating mantissa(cf. §20.1), we obtain the condition16 + log2 N + few × log2 log2 N < it(20.6.3)We see that single precision, say with it = 24, is inadequate for any interestingvalue of N , while double precision, say with it = 53, allows N to be greaterthan 106 , corresponding to some millions of decimal digits.
The following routinetherefore presumes double precision versions of realft (§12.3) and four1 (§12.2),here called drealft and dfour1. (These routines are included on the NumericalRecipes diskettes.)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.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
