chap32.htm

程序设计的艺术,此书是英文版的,经久不衰,非常的经典

<img src="784_c.gif"><P>
form a group under multiplication (see Section 33.3). This group has the same structure as the additive group (<B>Z</B><I><SUB>n</SUB>,</I> +) modulo <I>n</I>, since <img src="784_d.gif"> implies that <img src="784_e.gif">. Similarly, <img src="784_f.gif">. Essential properties of the complex <I>n</I>th roots of unity are given in the following lemmas.<P>
<a name="0992_1b17">Lemma 32.3<a name="0992_1b17"><P>
<a name="0992_1b13">For any integers <I>n</I> <img src="../images/gteq.gif"> 0, <I>k</I> <img src="../images/gteq.gif"> 0, and <I>d</I> &gt; 0,<P>
<img src="785_a.gif"><P>
<h4><a name="0992_1b18">(32.7)<a name="0992_1b18"></sub></sup></h4><P>
<I><B>Proof</I></B>     The lemma follows directly from equation (32.6), since<P>
<img src="785_b.gif"><P>
<a name="0992_1b19">Corollary 32.4<a name="0992_1b19"><P>
For any even integer <I>n</I> &gt; 0,<P>
<img src="785_c.gif"><P>
<I><B>Proof</I></B>     The proof is left as Exercise 32.2-1.      <P>
<a name="0992_1b1a"><a name="0992_1b14">Lemma 32.5<a name="0992_1b1a"><P>
If <I>n</I> &gt; 0 is even, then the squares of the <I>n</I> complex <I>n</I>th roots of unity are the <I>n</I>/2 complex (<I>n</I>/2)th roots of unity.<P>
<I><B>Proof</I></B>     By the cancellation lemma, we have <img src="785_d.gif">, for any non-negative integer <I>k</I>. Note that if we square all of the complex <I>n</I>th roots of unity, then each (<I>n</I>/2)th root of unity is obtained exactly twice, since<P>
<img src="785_e.gif"><P>
Thus, <img src="785_f.gif"> have the same square. This property can also be proved using Corollary 32.4, since <img src="785_g.gif"> implies<I><SUP> </I></SUP><img src="785_h.gif">, and thus <img src="785_i.gif">.      <P>
As we shall see, the halving lemma is essential to our divide-and-conquer approach for converting between coefficient and point-value representations of polynomials, since it guarantees that the recursive subproblems are only half as large.<P>
<a name="0992_1b1b">Lemma 32.6<a name="0992_1b1b"><P>
<a name="0992_1b15">For any integer <I>n</I> <img src="../images/gteq.gif"> 1 and nonnegative integer <I>k</I> not divisible by <I>n</I>,<P>
<img src="786_a.gif"><P>
<I><B>Proof</I></B>     Because equation (3.3) applies to complex values,<P>
<img src="786_b.gif"><P>
Requiring that <I>k</I> not be divisible by <I>n</I> ensures that the denominator is not 0, since <img src="786_c.gif"> only when <I>k</I> is divisible by <I>n</I>.      <P>
<P>







<h2>The DFT</h2><P>
Recall that we wish to evaluate a polynomial<P>
<img src="786_d.gif"><P>
of degree-bound <I>n</I> at <img src="786_e.gif"> (that is, at the <I>n</I> complex <I>n</I>th roots of unity).<SUP>2</SUP> Without loss of generality, we assume that <I>n</I> is a power of 2, since a given degree-bound can always be raised--we can always add new high-order zero coefficients as necessary. We assume that <I>A</I> is given in coefficient form: <I>a </I>= (<I>a</I><SUB>0</SUB>, <I>a</I><SUB>1</SUB>, . . . , <I>a<SUB>n-</I>1</SUB>). Let us define the results <I>y<SUB>k</I></SUB>, for <I>k</I> = 0, 1, . . . , <I>n</I> - 1, by<P>
<img src="786_f.gif"><P>
<h4><a name="0993_1b18">(32.8)<a name="0993_1b18"></sub></sup></h4><P>
<a name="0993_1b16"><a name="0993_1b17">The vector <I>y</I> = (<I>y</I><SUB>0</SUB>, <I>y</I><SUB>1</SUB>, . . . , <I>y<SUB>n</I>-1</SUB>) is the <I><B>Discrete Fourier Transform (DFT) </I></B>of the coefficient vector <I>a</I> = (<I>a</I><SUB>0</SUB>, <I>a</I><SUB>1</SUB>, . . . , <I>a<SUB>n</I>-1</SUB>). We also write <I>y</I> = DFT<I><SUB>n</I></SUB>(<I>a</I>).<P>
<SUP>2</SUP>The length <I>n</I> is actually what we referred to as 2<I>n </I>in Section 32.1, since we double the degree-bound of the given polynomials prior to evaluation. In the context of polynomial multiplication, therefore, we are actually working with complex (2<I>n</I>)th roots of unity.+<P>
<P>







<h2>The FFT</h2><P>
<a name="0994_1b18"><a name="0994_1b19">By using a method known as the <I><B>Fast Fourier Transform (FFT),</I></B> which takes advantage of the special properties of the complex roots of unity, we can compute<I> </I>DFT<I><SUB>n</I></SUB>(<I>a</I>)<SUB> </SUB>in time <img src="../images/bound.gif">(<I>n</I> lg <I>n</I>), as opposed to the <img src="../images/bound.gif">(<I>n</I><SUP>2</SUP>) time of the straightforward method.<P>
The FFT method employs a divide-and-conquer strategy, using the even-index and odd-index coefficients of <I>A</I>(<I>x</I>) separately to define the two new degree-bound <I>n</I>/2 polynomials <I>A</I><SUP>[0]</SUP>(<I>x</I>) and <I>A</I><SUP>[1]</SUP>(<I>x</I>):<P>
<pre><I>A</I><SUP>[0]</SUP>(<I>x</I>)  =  <I>a</I><SUB>0</SUB> + <I>a</I><SUB>2</SUB><I>x</I> + <I>a</I><SUB>4</SUB><I>x</I><SUP>2</SUP> + <img src="../images/dot10.gif"> <img src="../images/dot10.gif"> <img src="../images/dot10.gif"> + <I>a<SUB>n</I>-2</SUB><I>x<SUP>n</I>/2-1</SUP>,</sub></sup></pre><P>
<pre><I>A</I><SUP>[1]</SUP>(<I>x</I>)  =  <I>a</I><SUB>1</SUB> + <I>a</I><SUB>3</SUB><I>x</I> + <I>a</I><SUB>5</SUB><I>x</I><SUP>2</SUP> + <img src="../images/dot10.gif"> <img src="../images/dot10.gif"> <img src="../images/dot10.gif"> + <I>a<SUB>n</I>-1</SUB><I>x<SUP>n</I>/2-1</SUP>.</sub></sup></pre><P>
Note that <I>A</I><SUP>[0]</SUP> contains all the even-index coefficients of <I>A</I> (the binary representation of the index ends in 0) and <I>A</I><SUP>[1]</SUP> contains all the odd-index coefficients (the binary representation of the index ends in 1). It follows that<P>
<pre><I>A</I>(<I>x</I>) = <I>A</I><SUP>[0]</SUP>(<I>x</I><SUP>2</SUP>) + <I>xA</I><SUP>[1]</SUP>(<I>x</I><SUP>2</SUP>),</sub></sup></pre><P>
<h4><a name="0994_1b1b">(32.9)<a name="0994_1b1b"></sub></sup></h4><P>
so that the problem of evaluating <I>A</I>(<I>x</I>) at <img src="787_a.gif"> reduces to<P>
1.     evaluating the degree-bound <I>n</I>/2 polynomials <I>A</I><SUP>[0]</SUP>(<I>x</I>) and <I>A</I><SUP>[1]</SUP>(<I>x</I>) at the points<P>
<img src="787_b.gif"><P>
<h4><a name="0994_1b1c">(32.10)<a name="0994_1b1c"></sub></sup></h4><P>
and then<P>
2.     combining the results according to equation (32.9).<P>
<a name="0994_1b1a">By the halving lemma, the list of values (32.10) consists not of <I>n</I> distinct values but only of the <I>n</I>/2 complex (<I>n</I>/2)th roots of unity, with each root occurring exactly twice. Therefore, the polynomials <I>A</I><SUP>[0]</SUP> and <I>A</I><SUP>[1]</SUP> of degree-bound <I>n</I>/2 are recursively evaluated at the <I>n</I>/2 complex (<I>n</I>/2)th roots of unity. These subproblems have exactly the same form as the original problem, but are half the size. We have now successfully divided an <I>n</I>-element DFT<I><SUB>n</I></SUB> computation into two <I>n</I>/2-element DFT<I><SUB>n</I>/2</SUB> computations. This decomposition is the basis for the following recursive FFT algorithm, which computes the DFT of an <I>n</I>-element vector <I>a</I> = (<I>a</I><SUB>0</SUB><I>, a</I><SUB>1</SUB>, . . . , <I>a<SUB>n - </I>1</SUB>), where <I>n</I> is a power of 2.<P>
<img src="788_a.gif"><P>
The <FONT FACE="Courier New" SIZE=2>RECURSIVE</FONT>-FFT procedure works as follows. Lines 2-3 represent the basis of the recursion; the DFT of one element is the element itself, since in this case<P>
<img src="788_b.gif"><P>
Lines 6-7 define the coefficient vectors for the polynomials <I>A</I><SUP>[0]</SUP> and <I>A</I><SUP>[1]</SUP>. Lines 4, 5, and 13 guarantee that <I>w</I> is updated properly so that whenever lines 11-12 are executed, <img src="788_c.gif">. (Keeping a running value of <I>w</I> from iteration to iteration saves time over computing <img src="788_d.gif"> from scratch each time through the <B>for</B> loop.) Lines 8-9 perform the recursive DFT<I><SUB>n</I>/2</SUB> computations, setting, for <I>k</I> = 0, 1, . . . , <I>n</I>/2 - 1,<P>
<img src="788_e.gif"><P>
or, since <img src="788_f.gif"> by the cancellation lemma,<P>
<img src="788_g.gif"><P>
Lines 11-12 combine the results of the recursive DFT<I><SUB>n</I>/2</SUB> calculations. For <I>y</I><SUB>0</SUB>, <I>y</I><SUB>1</SUB>, . . . , <I>y<SUB>n</I>/2 - 1</SUB>, line 11 yields<P>
<img src="788_h.gif"><P>
where the last line of this argument follows from equation (32.9). For <I>y<SUB>n</I>/2</SUB>,<I>y<SUB>n</I>/2+1</SUB>, . . . , <I>y<SUB>n </I>- 1</SUB>, letting <I>k</I> = 0, 1, . . . , <I>n</I>/2 - 1, line 12 yields<P>
<img src="789_a.gif"><P>
The second line follows from the first since <img src="789_b.gif">. The fourth line follows from the third because <img src="789_c.gif"> implies <img src="789_d.gif">. The last line follows from equation (32.9). Thus, the vector <I>y</I> returned by <FONT FACE="Courier New" SIZE=2>RECURSIVE</FONT>-FFT is indeed the DFT of the input vector <I>a</I>.<P>
To determine the running time of procedure <FONT FACE="Courier New" SIZE=2>RECURSIVE</FONT>-FFT, we note that exclusive of the recursive calls, each invocation takes time <img src="../images/bound.gif">(<I>n</I>), where <I>n</I> is the length of the input vector. The recurrence for the running time is therefore<P>
<pre><I>T</I>(<I>n</I>) = 2<I>T</I>(<I>n</I>/2) + <img src="../images/bound.gif">(<I>n</I>)</sub></sup></pre><P>
<pre>= <img src="../images/bound.gif">(<I>n</I> lg <I>n</I>) .</sub></sup></pre><P>
Thus, we can evaluate a polynomial of degree-bound <I>n</I> at the complex <I>n</I>th roots of unity in time <img src="../images/bound.gif"><I></I>(<I>n</I> lg <I>n</I>) using the Fast Fourier Transform.<P>
<P>







<h2>Interpolation at the complex roots of unity</h2><P>
<a name="0995_1b1b"><a name="0995_1b1c">We now complete the polynomial multiplication scheme by showing how to interpolate the complex roots of unity by a polynomial, which enables us to convert from point-value form back to coefficient form. We interpolate by writing the DFT as a matrix equation and then looking at the form of the matrix inverse.<P>
From equation (32.4), we can write the DFT as the matrix product <I>y</I> = <I>V<SUB>n</SUB>a</I>, where <I>V<SUB>n</I></SUB> is a Vandermonde matrix containing the appropriate powers of <I>w<SUB>n</I></SUB>:<P>
<img src="789_e.gif"><P>
The (<I>k, j</I>) entry of <img src="789_f.gif">, for <I>j,k</I> = 0, 1, . . . , <I>n</I> - 1, and the exponents of the entries of <I>V<SUB>n</I></SUB> form a multiplication table.<P>
For the inverse operation, which we write as <img src="789_g.gif">, we proceed by multiplying <I>y</I> by the matrix <img src="789_h.gif">, the inverse of <I>V<SUB>n</I></SUB>.<P>
<a name="0995_1b1e">Theorem 32.7<a name="0995_1b1e"><P>