Python代碼實現了fft與逆fft。
參考資料:
https://walkccc.github.io/CLRS/Chap30/30.1/
https://sites.math.rutgers.edu/~ajl213/CLRS/Ch30.pdf

30 Polynomials and the FFT

30.1 Representing polynomials

1.Multiply the polynomials?and??using equations (30.1)?and (30.2).

2.Another way to evaluate a polynomial?A(x) of degree-bound?n?at a given point??is to divide?A(x) by the polynomial?, obtaining a quotient polynomial q(x)?of degree-bound?n - 1 and a remainder?r, such that

Clearly,?. Show how to compute the remainder?r?and the coefficients of?q(x) in time? from???and the coefficients of?A.

Let?A?be the matrix with?1's on the diagonal,??'s on the super diagonal, and?0's everywhere else. Let?q?be the vector?. If??then we need to solve the matrix equation?Aq = a to compute the remainder and coefficients. Since?A?is tridiagonal, Problem 28-1 (e) tells us how to solve this equation in linear time.

3.Derive a point-value representation for??from a point-value representation for, assuming that none of the points is?0.

4.Prove that?n?distinct point-value pairs are necessary to uniquely specify a polynomial of degree-bound?nn, that is, if fewer than?nn?distinct point-value pairs are given, they fail to specify a unique polynomial of degree-bound?n. Hint:?Using Theorem 30.1, what can you say about a set of?n - 1 point-value pairs to which you add one more arbitrarily chosen point-value pair?)

5.Show how to use equation (30.5)?to interpolate in time. Hint:?First compute the coefficient representation of the polynomial?and then divide by?as necessary for the numerator of each term; see Exercise 30.1-2. You can compute each of the?nn?denominators in time?O(n).

6.Explain what is wrong with the "obvious" approach to polynomial division using a point-value representation, i.e., dividing the corresponding?yy?values. Discuss separately the case in which the division comes out exactly and the case in which it doesn't.

If we wish to compute P/Q?but Q?takes on the value zero at some of these points, then we can't carry out the "obvious" method. However, as long as all point value pairs??we choose for?Q?are such that?, then the approach comes out exactly as we would like.

7.Consider two sets?A?and?B, each having?n?integers in the range from?0?to 10n. We wish to compute the?Cartesian sum?of?A?and?B, defined by

Note that the integers in?C?are in the range from?0?to 20n. We want to find the elements of?C?and the number of times each element of?C?is realized as a sum of elements in?A?and?B. Show how to solve the problem in O(nlgn)?time.Hint:?Represent?A?and?B?as polynomials of degree at most?10n.

30.2 The DFT and FFT

1.Prove Corollary 30.4.

根據引理30.3消去引理，這里設d=n/2，k=2，就可以得出：

2.Compute the DFT?of the vector (0,1,2,3).

3.Do Exercise 30.1-1 by using the Θ(nlgn)-time scheme.

4.Write pseudocode to compute??Θ(nlgn)?time.

Python代碼：

1).DFT:

import numpy as npdef recursive_fft(a):n = len(a)if n == 1:return aw_n = np.exp(complex(0,2*np.pi/n))w = 1a0 = a[0:n:2]a1 = a[1:n:2]y0 = recursive_fft(a0)y1 = recursive_fft(a1)y = np.zeros(n,dtype='complex')for k in np.arange(0,n//2):y[k] = y0[k] + w*y1[k]y[k + n//2] = y0[k] - w*y1[k]w = w*w_nreturn yprint(recursive_fft([0,1,2,3]))

示例計算了本節習題2，輸出為：[ 6.+0.j -2.-2.j -2.+0.j -2.+2.j]

2).逆DFT：

import numpy as npdef recursive_inverse_fft_base(y):n = len(y)if n == 1:return yw_n = np.exp(complex(0,-2*np.pi/n))w = 1y0 = y[0:n:2]y1 = y[1:n:2]a0 = recursive_inverse_fft_base(y0)a1 = recursive_inverse_fft_base(y1)a = np.zeros(n,dtype='complex')for k in np.arange(0,n//2):a[k] = (a0[k] + w*a1[k])a[k + n//2] = (a0[k] - w*a1[k])w = w*w_nreturn adef recursive_inverse_fft(y):a = recursive_inverse_fft_base(y)return a / len(y)print(recursive_inverse_fft([6,-2-2j,-2,-2+2j]))

這是對本節第2題求逆，結果應該是[0,1,2,3].代碼運行結果為：

[0.+0.000000e+00j 1.-6.123234e-17j 2.+0.000000e+00j 3.+6.123234e-17j]

第2項和第4項有一個極小的j，這是計算浮點數的誤差導致，可以認為是0，結果就是[0,1,2,3].

3).結合fft的代碼計算本節第3題：

a = np.array([7,-1,1,-10,0,0,0,0]) b = np.array([8,0,-6,3,0,0,0,0]) ya = recursive_fft(a) yb = recursive_fft(b) yc = ya*yb c = recursive_inverse_fft(yc) print(c)

輸出結果：

[ 56.-3.55271368e-15j ?-8.+8.75870054e-15j -34.+0.00000000e+00j
?-53.-6.98234370e-15j ?-9.+3.55271368e-15j ?63.-5.45215418e-15j
?-30.+0.00000000e+00j ? 0.+3.67579734e-15j]

因為我們做了后補0，因此需要把輸出中最后一項0去掉，系數從高到低就是[56,-8,-34,-53,-9,63,-30]

5.Describe the generalization of theFFT?procedure to the case in which?nn?is a power of?3. Give a recurrence for the running time, and solve the recurrence.

6.Suppose that instead of performing an?nn-element FFT?over the field of complex numbers (where?nn?is even), we use the ring Zm??of integers modulo?m, where??and?t?is an arbitrary positive integer. Use??instead of??as a principal nth root of unity, modulo?mm. Prove that the DFT?and the inverse DFT?are well defined in this system.

7.Given a list of values??(possibly with repetitions), show how to find the coefficients of a polynomial P(x)?of degree-bound n+1?that has zeros only at?(possibly with repetitions). Your procedure should run in time?. (Hint:?The polynomial? P(x)?has a zero at??if and only if P(x)?is a multiple of

8.The?chirp transform?of a vector?is the vector?, whereand?z?is any complex number. The?DFT?is therefore a special case of the chirp transform, obtained by taking??. Show how to evaluate the chirp transform in time??for any complex number?z.? Hint:?Use the equation

to view the chirp transform as a convolution.)

30.3 Efficient FFT implementations

1.Show how ITERATIVE-FFT?computes the?DFT?of the input vector (0,2,3,?1,4,5,7,9).

2.Show how to implement an FFT?algorithm with the bit-reversal permutation occurring at the end, rather than at the beginning, of the computation. Hint:?Consider the inverse DFT.

We can consider ITERATIVE-FFT as working in two steps, the copy step (line 1) and the iteration step (lines 4 through 14). To implement in inverse
iterative FFT algorithm, we would need to first reverse the iteration step, then reverse the copy step. We can do this by implementing the INVERSEINTERATIVE-FFT algorithm exactly as we did with the FFT algorithm, but this time creating the iterative version of RECURSIVE-FFT-INV as given in exercise 30.2-4.

3.How many times does ITERATIVE-FFT?compute twiddle factors in each stage? Rewrite ITERATIVE-FFT?to compute twiddle factors only??times in stage?s.

4.Suppose that the adders within the butterfly operations of the FFT?circuit sometimes fail in such a manner that they always produce a zero output, independent of their inputs. Suppose that exactly one adder has failed, but that you don't know which one. Describe how you can identify the failed adder by supplying inputs to the overall FFT?circuit and observing the outputs. How efficient is your method?

First observe that if the input to a butterfly operation is all zeros, then so is the output. Our initial input will be () where if 0 ≤ i ≤ n/2 - 1, and i otherwise. By examining the diagram on page 919, we?see that the zeros will propogate as half the input to each butterfly operation in the final stage. Thus, if any of the output values are 0 we know that one of
the ’s is zero, where k ≥ n/2 - 1. If this is the case, then the faulty butterfly added occurs along the wire connected to yk. If none of the output values are zero, the faulty butterfly adder must occur as one which only deals with the first n/2 - 1 wires. In this case, we rerun the test with input such that ?if 0 ≤ i ≤ n=2 - 1 and 0 otherwise. This time, we will know for sure which of the first n/2-1 wires touches the faulty butterfly adder. Since only lg n adders touch a given wire, we have reduced the problem of size n to a problem of size lg n. Thus, at most ??tests are required.
?

Problems

1.Divide-and-conquer multiplication:

a.?Show how to multiply two linear polynomials?ax + b? and?cx + d? using only three multiplications. Hint:?One of the multiplications is

b.?Give two divide-and-conquer algorithms for multiplying two polynomials of degree-bound?nn?inΘ(nlg3)?time. The first algorithm should divide the input polynomial coefficients into a high half and a low half, and the second algorithm should divide them according to whether their index is odd or even.

c.?Show how to multiply two?nn-bit integers in O(nlg3)?steps, where each step operates on at most a constant number of?1-bit values.

2.A?Toeplitz matrix?is an? n×n?matrix??such that???for??and?.

a.?Is the sum of two Toeplitz matrices necessarily Toeplitz? What about the product?

b.?Describe how to represent a Toeplitz matrix so that you can add two n×n?Toeplitz matrices in? O(n)?time.

c.?Give an? O(nlgn)-time algorithm for multiplying an? n×n?Toeplitz matrix by a vector of length?n. Use your representation from part (b).

d.?Give an efficient algorithm for multiplying two n×n?Toeplitz matrices. Analyze its running time.

3. Multidimensional fast Fourier transform

We can generalize the?1-dimensional discrete Fourier transform defined by equation(30.8)?to?d?dimensions. The input is a?d-dimensional array??whose dimensions are??, where?. We define the?d-dimensional discrete Fourier transform by the equation

for?

a.?Show that we can compute a?dd-dimensional DFT?by computing 1-dimensional DFTs on each dimension in turn. That is, we first compute ???separate?1-dimensional DFTs along dimension?1. Then, using the result of the DFTs along dimension?1?as the input, we compute ?separate?1-dimensional DFTs along dimension?2. Using this result as the input, we compute???separate?1-dimensional DFTs along dimension?3, and so on, through dimension?d.

b.?Show that the ordering of dimensions does not matter, so that we can compute a?d-dimensional DFT?by computing the?1-dimensional DFTs in any order of the?d?dimensions.

c.?Show that if we compute each?1-dimensional DFT?by computing the fast Fourier transform, the total time to compute a?d-dimensional DFT?is O(nlgn), independent of?d.

4.Evaluating all derivatives of a polynomial at a point

Given a polynomial? A(x)?of degree-bound?n, we define its?tth derivative by

From the coefficient representation??of?A(x) and a given point??, we wish to determine??for?.

a.?Given coefficients??such that

show how to compute??for??in?O(n) time.

b.?Explain how to find???in??time, given??for?.

c.?Prove that

where??and

d.?Explain how to evaluate???for???in O(nlgn)?time. Conclude that we can evaluate all nontrivial derivatives of A(x)?at??in O(nlgn)?time.

5. Polynomial evaluation at multiple points

We have seen how to evaluate a polynomial of degree-bound?nn?at a single point in O(n)?time using Horner's rule. We have also discovered how to evaluate such a polynomial at all?n?complex roots of unity in O(nlgn)?time using the FFT. We shall now show how to evaluate a polynomial of degree-bound?n?at?n?arbitrary points in??time.

To do so, we shall assume that we can compute the polynomial remainder when one such polynomial is divided by another inO(nlgn)?time, a result that we state without proof. For example, the remainder of???when divided by??is

.

Given the coefficient representation of a polynomial??and?n?points??, we wish to compute the?n?values?. For?, define the polynomials??and?. Note that??has degree at most?j - i.

a.?Prove that??for any point?zz.

b.?Prove that??and that?.

c.?Prove that for?, we have??and?

d.?Give an?-time algorithm to evaluate?

6. FFT using modular arithmetic

As defined, the discrete Fourier transform requires us to compute with complex numbers, which can result in a loss of precision due to round-off errors. For some problems, the answer is known to contain only integers, and by using a variant of the FFT?based on modular arithmetic, we can guarantee that the answer is calculated exactly. An example of such a problem is that of multiplying two polynomials with integer coefficients. Exercise 30.2-6 gives one approach, using a modulus of length Ω(n)?bits to handle a DFT?on?n?points. This problem gives another approach, which uses a modulus of the more reasonable lengthO(lgn); it requires that you understand the material of Chapter 31. Let?n?be a power of?2.

a.?Suppose that we search for the smallest?k?such thatp=kn+1?is prime. Give a simple heuristic argument why we might expect?k?to be approximately lnn. (The value of?k?might be much larger or smaller, but we can reasonably expect to examineO(lgn)?candidate values of?k?on average.) How does the expected length of?pp?compare to the length of?n?

Let?g?be a generator of, and let?.

b.?Argue that the DFT?and the inverseDFT?are well-defined inverse operations modulo?p, where?w?is used as a principal?nth root of unity.

c.?Show how to make the FFT?and its inverse work modulo?p?in time O(nlgn), where operations on words of O(lgn)?bits take unit time. Assume that the algorithm is given?p?and?w.

d.?Compute the DFT?modulo p=17?of the vector (0,5,3,7,7,2,1,6). Note that g=3?is a generator of??.

b,c,d omitted.

?

?

?

總結

以上是生活随笔為你收集整理的算法导论第三版 第30章习题答案的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。