Arndt - Algorithms for Programmers (523138), страница 3
Текст из файла (страница 3)
. . . . . . . . . . . . . . . . . . . . . . . . 28313.3.2 Square root extraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28413.3.3 Cube root extraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28513.4 Square root extraction for rationals . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . 28513.5 A general procedure for the inverse n-th root . . . . . . . . . . . . . . . . . . . . . . . . . 28713.6 Re-orthogonalization of matrices . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . 29013.7 n-th root by Goldschmidt’s algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29213.8 Iterations for the inversion of a function . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29413.8.1 Householder’s formula . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29413.8.2 Schröder’s formula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29513.8.3 Dealing with multiple roots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29613.8.4 A general scheme .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29813.8.5 Improvements by the delta squared process . . . . . . . . . . . . . . . . . . . . . . 30013.8.6 Improvements of the delta squared process * . . . . .
. . . . . . . . . . . . . . . . 30213.9 Transcendental functions & the AGM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30313.9.1 The AGM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30313.9.2 log . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30513.9.3 exp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30613.9.4 sin, cos and tan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . 30713.9.5 Elliptic K . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30713.9.6 Elliptic E . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30813.10Computation of π/ log(q) . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . 30913.11Computation of q = exp(−π K 0 /K) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31013.12Iterations for high precision computations of π . . . . . . . . . . . . . . . . . . . . .
. . . 311CONTENTS713.13The binary splitting algorithm for rational series . . . . . . . . . . . . . . . . . . . . . . . 31613.14The magic sumalt algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31813.15Chebyshev polynomials * . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32113.16Continued fractions *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32313.17Some hypergeometric identities * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32413.17.1 Definition .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32413.17.2 Transformations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32513.17.3 Examples: elementary functions . . . . . . . . . . . . . . . . .
. . . . . . . . . . . 32813.17.4 Elliptic K and E . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 330A List of important Symbols332B The pseudo language Sprache334Bibliography336Some important remarks. . . about this document.This is a draft of what is intended to turn into a book about selected algorithms. The audience in mindare programmers who are interested in the treated algorithms and actually want to create and understandworking and reasonably optimized code.The printable full version will always stay online for free download. The referenced sources are online aspart of FXT (fast transforms and low level routines) or hfloat (arithmetical algorithms).The reader is welcome to criticize and make suggestions.
Thanks go to those1 who helped to improvethis document so far! Thanks also to the people who share their ideas or source code on the net. I try togive due references to original sources and authors wherever I can. However, I am in no way an expertfor history of algorithms and I pretty sure will never be one. So if you feel that a reference is missingsomewhere, let me know.New sections appear as soon as they contain anything useful, sometimes just listings or remarks outliningwhat is to appear. A ”TBD: something to be done” is a reminder to myself to fill in something that ismissing or would be nice to have.The style varies from chapter to chapter which I do not consider bad per se: while some topics (as fastFourier transforms) need a clear and explicit introduction others (like the bit wizardry chapter) seem tobe best presented by basically showing the code with just a few comments.
Still other parts (like thechapter about sorting and searching) are presented elsewhere extremely well so I only the basic ideas areintroduced shortly (accompanied by code) and references for further studies are given.The pseudo language Sprache is used when I see a clear advantage to do so, mainly when the corresponding C++ does not appear to be self explanatory. Larger pieces of code are presented inC++. A tiny starter about C++ (some good reasons in favor of C++ and some of the very basicsof classes/overloading/templates) might be included.
C programmers do not need to be shocked bythe ‘++’: only an rather minimal set of the C++ features is used.Enjoy reading !1 inparticular André Piotrowski and Edith Parzefall.”Why make things difficult, when it is possible to make them crypticand totally illogic, with just a little bit more effort?”– Aksel Peter JørgensenChapter 1The Fourier transform1.1The discrete Fourier transformThe discrete Fourier transform (DFT or simply FT) of a complex sequence a of length n is defined ascck=:=F [a]1√n(1.1)n−1Xax z +x kwhere z = e± 2 π i/n(1.2)x=0z is an n-th root of unity: z n = 1.Back-transform (or inverse discrete Fourier transform IDFT or simply IFT) is thenF −1 [c]a =ax1√n=n−1X(1.3)ck z −x k(1.4)k=0To see this, consider element y of the IFT of the FT of a:F −1 [F [a]]y==n−1n−11 X 1 X√√(ax z x k ) z −y knn x=0k=0X1 Xax(z x−y )kn x(1.5)(1.6)kPAs k (z x−y )k = n for x = y and zero else (because z is an n-th root of unity). Therefore the wholeexpression is equal to1 Xax δx,ynnx=ay(1.7)1 (x = y)0 (x =6 y)(1.8)where½δx,y=Here we will call the FT with the plus in the exponent the forward transform.
The choice is actuallyarbitrary1 .1 Electricalengineers prefer the minus for the forward transform, mathematicians the plus.10CHAPTER 1. THE FOURIER TRANSFORM11The FT is a linear transform, i.e. for α, β ∈ CF [α a + β b] =α F [a] + β F [b](1.9)For the FT Parseval’s equation holds, let c = F [a], thenn−1Xa2x=x=0n−1Xc2k(1.10)k=0The normalization factor √1n in front of the FT sums is sometimes replaced by a single n1 in front of theinverse FT sum which is often convenient in computation.
Then, of course, Parseval’s equation has to bemodified accordingly.A straight forward implementation of the discrete Fourier transform, i.e. the computation of n sums eachof length n requires ∼ n2 operations:void slow_ft(Complex *f, long n, int is){Complex h[n];const double ph0 = is*2.0*M_PI/n;for (long w=0; w<n; ++w){Complex t = 0.0;for (long k=0; k<n; ++k){t += f[k] * SinCos(ph0*k*w);}h[w] = t;}copy(h, f, n);}[FXT: slow ft in slow/slowft.cc] is must be +1 (forward transform) or −1 (backward transform),SinCos(x) returns a Complex(cos(x), sin(x)).A fast FourierPm transform (FFT) algorithm is an algorithm that improves the operation count to proportional n k=1 (pk − 1), where n = p1 p2 · · · pm is a factorization of n. In case of a power n = pm thevalue computes to n (p − 1) logp (n).
In the special case p = 2 even n/2 log2 (n) (complex) multiplicationssuffice. There are several different FFT algorithms with many variants.1.2Symmetries of the Fourier transformA bit of notation turns out to be useful:Let a be the sequence a (length n) reversed around element with index n/2:a0an/2ak:= a0:= an/2:= an−kif n even(1.11)(1.12)(1.13)Let aS , aA be the symmetric, antisymmetric part of the sequence a, respectively:aSaA:= a + a:= a − a(1.14)(1.15)(The elements with indices 0 and n/2 of aA are zero). Now let a ∈ R (meaning that each element of a is∈ R), thenF [aS ]∈F [aS ] =F [aA ] ∈F [aA ] =R(1.16)F [aS ]iR−F [aA ](1.17)(1.18)(1.19)CHAPTER 1.
THE FOURIER TRANSFORM12i.e. the FT of a real symmetric sequence is real and symmetric and the FT of a real antisymmetricsequence is purely imaginary and antisymmetric. Thereby the FT of a general real sequence is thecomplex conjugate of its reversed:F [a]=F [a]∗f ora∈R(1.20)Similarly, for a purely imaginary sequence b ∈ iR:F [bS ] ∈F [bS ] =F [bA ] ∈F [bA ] =iRF [bS ]R−F [bA ](1.21)(1.22)(1.23)(1.24)The FT of a complex symmetric/antisymmetric sequence is symmetric/antisymmetric, respectively.1.3Summary of definitions of Fourier transforms *This section summarizes the definitions of the continuous, semi-continuous and the discrete Fouriertransform.The continuous Fourier transformThe (continuous) Fourier transform (FT) of a function f : Cn → Cn , ~x 7→ f (~x) is defined byZ1F (~ω ) := ¡√ ¢nf (~x) eσ i ~x ω~ dn x2πCn(1.25)where σ = ±1. The FT is is a unitary transform.Its inverse (‘back-transform’) isf (~x) =√Z12πnF (~ω ) e−σ ~x ω~ dn ω(1.26)Cni.e.
the complex conjugate transform.For the 1-dimensional case one hasF (ω) =f (x) =Z +∞1√f (x)eσ x ω dx2 π −∞Z +∞1√F (ω) e−σ x ω dω2 π −∞(1.27)(1.28)The ‘frequency’-form isZfˆ(ν) =+∞−∞Z +∞f (x) =−∞f (x)eσ 2 π i x ν dx(1.29)fˆ(ν) e−σ 2 π i x ν dν(1.30)CHAPTER 1. THE FOURIER TRANSFORM13The semi-continuous Fourier transformFor periodic functions defined on a interval L ∈ R, f : L → R, x 7→ f (x) one has the semi-continuousFourier transform:Z1ck := √f (x) eσ 2 π i k x/L dx(1.31)L LThen½k=+∞f (x)if f continuous at x1 X√ck e−σ 2 π i k x/L =(1.32)f (x+0)+f (x−0)elseL2k=−∞Another (equivalent) form is given byak:=bk:=f (x)=Z12πkx√f (x) cosdx,k = 0, 1, 2, . . .LL LZ12πkx√f (x) sindx,k = 1, 2, . .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
