chap32.htm

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

Thus, <I>n</I>-point evaluation and interpolation are well-defined inverse operations that transform between the coefficient representation of a polynomial and a point-value representation.<SUP>1</SUP> The algorithms described above for these problems take time <img src="../images/bound.gif">(<I>n</I><SUP>2</SUP>).<P>
<SUP>1</SUP>Interpolation is a notoriously tricky problem from the point of view of numerical stability. Although the approaches described here are mathematically correct, small differences in the inputs or round-off errors during computation can cause large differences in the result.<P>
The point-value representation is quite convenient for many operations on polynomials. For addition, if <I>C</I>(<I>x</I>) = <I>A</I>(<I>x</I>) + <I>B</I>(<I>x</I>), then <I>C</I>(<I>x<SUB>k</I></SUB>) = <I>A</I>(<I>x<SUB>k</I></SUB>) + <I>B</I>(<I>x<SUB>k</I></SUB>) for any point <I>x<SUB>k</I></SUB>. More precisely, if we have a point-value representation for <I>A</I>,<P>
<pre>{(<I>x</I><SUB>0</SUB>, <I>y</I><SUB>0</SUB>), (<I>x</I><SUB>1</SUB>, <I>y</I><SUB>1</SUB>), . . ., (<I>x<SUB>n</I>-1</SUB>, <I>y<SUB>n</I>-1</SUB>)} ,</sub></sup></pre><P>
and for <I>B</I>,<P>
<pre>{(x<SUB>0</SUB>, y'<SUB>0</SUB>), (x<SUB>1</SUB>, y'<SUB>1</SUB>), . . ., (x<I><SUB>n</I>-1</SUB>, y'<I><SUB>n</I>-1</SUB>)}</sub></sup></pre><P>
(note that <I>A</I> and <I>B</I> are evaluated at the <I>same n</I> points), then a point-value representation for <I>C</I> is<P>
<pre>{(x<SUB>0</SUB>, y<SUB>0 </SUB>+ y'<SUB>0</SUB>), (x<SUB>1</SUB>, y<SUB>1</SUB>+y'<SUB>1</SUB>), . . ., (x<SUB>n-1</SUB>, y<SUB>n-1</SUB>+y'<SUB>n-1</SUB>)}.</sub></sup></pre><P>
The time to add two polynomials of degree-bound <I>n</I> in point-value form is thus <img src="../images/bound.gif">(<I>n</I>).<P>
Similarly, the point-value representation is convenient for multiplying polynomials. If <I>C</I>(<I>x</I>)<I> = A</I>(<I>x</I>)<I>B</I>(<I>x</I>)<I>, </I>then<I> C</I>(<I>x<SUB>k</I></SUB>)<I> = A</I>(<I>x<SUB>k</I></SUB>)<I>B</I>(<I>x<SUB>k</I></SUB>) for any<SUB> </SUB>point <I>x<SUB>k</I></SUB>, and we can pointwise multiply a point-value representation for <I>A </I>by a point-value representation for <I>B</I> to obtain a point-value representation for <I>C</I>. We must face the problem, however, that the degree-bound of <I>C</I> is the sum of the degree-bounds for <I>A</I> and <I>B</I>. A standard point-value representation for <I>A</I> and <I>B</I> consists of <I>n</I> point-value pairs for each polynomial. Multiplying these together gives us <I>n</I> point-value pairs for <I>C</I>, but since the degree-bound of <I>C</I> is 2<I>n</I>, Theorem 32.1 implies that we need 2<I>n</I> point-value pairs for a point-value representation of <I>C</I>. We must therefore begin with "extended" point-value representations for <I>A</I> and for <I>B</I> consisting of 2<I>n</I> point-value pairs each. Given an extended point-value representation for <I>A</I>,<P>
<pre>{(<I>x</I><SUB>0</SUB>,<I>y</I><SUB>0</SUB>),(<I>x</I><SUB>1</SUB>, <I>y</I><SUB>1</SUB>),...,(<I>x</I><SUB>2<I>n</I>-1</SUB>,<I>y</I><SUB>2<I>n</I>-1</SUB>)},</sub></sup></pre><P>
and a corresponding extended point-value representation for <I>B</I>,<P>
<pre>{(<I>x</I><SUB>0</SUB>,<I>y</I>'<SUB>0</SUB>),(<I>x</I><SUB>1</SUB>,<I>y</I>'<SUB>1</SUB>),...,(<I>x</I><SUB>2<I>n</I>-1</SUB>,<I>y</I>'<SUB>2<I>n</I>-1</SUB>)},</sub></sup></pre><P>
then a point-value representation for C is<P>
<pre>{(<I>x</I><SUB>0</SUB>,<I>y</I><SUB>0</SUB><I>y</I>'<SUB>0</SUB>),(<I>x</I><SUB>1</SUB>,<I>y</I><SUB>1</SUB><I>y</I>'<SUB>1</SUB>),...,(<I>x</I><SUB>2<I>n</I>-1</SUB>,<I>y</I><SUB>2<I>n</I>-1</SUB><I>y</I>'<SUB>2<I>n</I>-1</SUB>)}.</sub></sup></pre><P>
Given two input polynomials in extended point-value form, we see that the time to multiply them to obtain the point-value form of the result is <img src="../images/bound.gif">(<I>n</I>), much less than the time required to multiply polynomials in coefficient form.<P>
Finally, we consider how to evaluate a polynomial given in point-value form at a new point. For this problem, there is apparently no approach that is simpler than converting the polynomial to coefficient form first, and then evaluating it at the new point.<P>
<P>







<h2>Fast multiplication of polynomials in coefficient form</h2><P>
<a name="098f_1b0a"><a name="098f_1b0b">Can we use the linear-time multiplication method for polynomials in point-value form to expedite polynomial multiplication in coefficient form? The answer hinges on our ability to convert a polynomial quickly from coefficient form to point-value form (evaluate) and vice-versa (interpolate).<P>
We can use any points we want as evaluation points, but by choosing the evaluation points carefully, we can convert between representations in only <img src="../images/bound.gif">(<I>n </I>lg <I>n</I>) time. As we shall see in Section 32.2, if we choose <FONT FACE="CG Times (W1)" SIZE=2>"</FONT>complex roots of unity" as the evaluation points, we can produce a point-value representation by taking the Discrete Fourier Transform (or DFT) of a coefficient vector. The inverse operation, interpolation, can be performed by taking the "inverse DFT<FONT FACE="CG Times (W1)" SIZE=2>"</FONT> of point-value pairs, yielding a coefficient vector. Section 32.2 will show how the FFT performs the DFT and inverse DFT operations in <img src="../images/bound.gif">(<I>n </I>lg <I>n</I>) time.<P>
Figure 32.1 shows this strategy graphically. One minor detail concerns degree-bounds. The product of two polynomials of degree-bound <I>n </I>is a polynomial of degree-bound 2<I>n</I>. Before evaluating the input polynomials <I>A</I> and <I>B</I>, therefore, we first double their degree-bounds to 2<I>n</I> by adding <I>n </I>high-order coefficients of 0. Because the vectors have 2<I>n</I> elements, we use &quot;complex (2<I>n</I>)th roots of unity,&quot; which are denoted by the <I>w</I><SUB>2<I>n</I></SUB> terms in Figure 32.1.<P>
Given the FFT, we have the following <img src="../images/bound.gif">(<I>n </I>lg <I>n</I>)-time procedure for multiplying two polynomials <I>A</I>(<I>x</I>) and <I>B</I>(<I>x</I>) of degree-bound <I>n</I>, where the input and output representations are in coefficient form. We assume that <I>n</I> is a power of 2; this requirement can always be met by adding high-order zero coefficients.<P>
<img src="782_a.gif"><P>
<h4><a name="098f_1b0c">Figure 32.1 A graphical outline of an efficient polynomial-multiplication process. Representations on the top are in coefficient form, while those on the bottom are in point-value form. The arrows from left to right correspond to the multiplication operation. The w<SUB>2n</SUB> terms are complex (2n)th roots of unity.<a name="098f_1b0c"></sub></sup></h4><P>
1.     <I>Double degree-bound:</I> Create coefficient representations of <I>A</I>(<I>x</I>) and <I>B</I>(<I>x</I>) as degree-bound 2<I>n</I> polynomials by adding <I>n</I> high-order zero coefficients to each.<P>
2.     <I>Evaluate:</I> Compute point-value representations of <I>A</I>(<I>x</I>) and <I>B</I>(<I>x</I>) of length 2<I>n</I> through two applications of the FFT of order 2<I>n</I>. These representations contain the values of the two polynomials at the (2<I>n</I>)th roots of unity.<P>
3.     <I>Pointwise multiply:</I> Compute a point-value representation for the polynomial <I>C</I>(<I>x</I>) = <I>A</I>(<I>x)B</I>(<I>x</I>) by multiplying these values together pointwise. This representation contains the value of <I>C</I>(<I>x</I>) at each (2<I>n</I>)th root of unity.<P>
4.     <I>Interpolate:</I> Create the coefficient representation of the polynomial <I>C</I>(<I>x</I>) through a single application of an FFT on 2<I>n</I> point-value pairs to compute the inverse DFT.<P>
Steps (1) and (3) take time <img src="../images/bound.gif">(<I>n</I>), and steps (2) and (4) take time <img src="../images/bound.gif">(<I>n </I>lg <I>n</I>). Thus, once we show how to use the FFT, we will have proven the following.<P>
<a name="098f_1b0d">Theorem 32.2<a name="098f_1b0d"><P>
The product of two polynomials of degree-bound <I>n</I> can be computed in time <img src="../images/bound.gif">(<I>n </I>lg <I>n</I>), with both the input and output representations in coefficient form.      <P>
<P>







<h2><a name="0990_1b10">Exercises<a name="0990_1b10"></h2><P>
<a name="0990_1b11">32.1-1<a name="0990_1b11"><P>
Multiply the polynomials <I>A</I>(<I>x</I>) = 7<I>x</I><SUP>3</SUP> - <I>x</I><SUP>2</SUP> + <I>x</I> - 10 and <I>B</I>(<I>x</I>) = 8<I>x</I><SUP>3</SUP> - 6<I>x</I> + 3 using equations (32.1) and (32.2).<P>
<a name="0990_1b12">32.1-2<a name="0990_1b12"><P>
<a name="0990_1b0c"><a name="0990_1b0d">Evaluating a polynomial <I>A</I>(<I>x</I>) of degree-bound <I>n</I> at a given point <I>x</I><SUB>0</SUB> can also be done by dividing <I>A</I>(<I>x</I>) by the polynomial (<I>x</I> - <I>x</I><SUB>0</SUB><I>) </I>to obtain a quotient polynomial <I>q</I>(<I>x</I>) of degree-bound <I>n</I> - 1 and a remainder <I>r</I>, such that<P>
<pre><I>A</I>(<I>x</I>) = <I>q</I>(<I>x</I>)(<I>x</I> - <I>x</I><SUB>0</SUB>) + <I>r </I>.</sub></sup></pre><P>
Clearly, <I>A</I>(<I>x</I><SUB>0</SUB>) = <I>r</I>. Show how to compute the remainder <I>r</I> and the coefficients of <I>q</I>(<I>x</I>) in time <img src="../images/bound.gif">(<I>n</I>) from <I>x</I><SUB>0</SUB> and the coefficients of <I>A</I>.<P>
<a name="0990_1b13">32.1-3<a name="0990_1b13"><P>
Derive a point-value representation for <img src="783_a.gif"> from a point-value representation for <img src="783_b.gif">, assuming that none of the points is 0.<P>
<a name="0990_1b14">32.1-4<a name="0990_1b14"><P>
Show how to use equation (32.5) to interpolate in time <img src="../images/bound.gif">(<I>n</I><SUP>2</SUP>). (<I>Hint</I>: First compute <img src="../images/piuc.gif"><I><SUB>j</I></SUB>(<I>x</I> - <I>x<SUB>k</I></SUB>) and <img src="../images/piuc.gif"><I><SUB>j</I></SUB>(<I>x<SUB>j</I></SUB> - <I>x<SUB>k</I></SUB>) and then divide by (<I>x</I> - <I>x<SUB>k</I></SUB>) and (<I>x<SUB>j</I></SUB> - <I>x<SUB>k</I></SUB>) as necessary for each term. See Exercise 32.1-2.)<P>
<a name="0990_1b15">32.1-5<a name="0990_1b15"><P>
Explain what is wrong with the "obvious" approach to polynomial division using a point-value representation. Discuss separately the case in which the division comes out exactly and the case in which it doesn't.<P>
<a name="0990_1b16">32.1-6<a name="0990_1b16"><P>
<a name="0990_1b0e"><a name="0990_1b0f">Consider two sets <I>A</I> and <I>B</I>, each having <I>n</I> integers in the range from 0 to 10<I>n</I>. We wish to compute the <I><B>Cartesian sum</I></B> of <I>A</I> and <I>B</I>, defined by<P>
<pre><I>C</I> = {<I>x</I> + <I>y</I>: <I>x</I> <img src="../images/memof12.gif"> <I>A</I> and <I>y</I> <img src="../images/memof12.gif"> <I>B</I>}.</sub></sup></pre><P>
Note that the integers in <I>C</I> are in the range from 0 to 20<I>n</I>. We want to find the elements of <I>C</I> and the number of times each element of <I>C</I> is realized as a sum of elements in <I>A</I> and <I>B</I>. Show that the problem can be solved in <I>O</I>(<I>n </I>lg <I>n</I>) time. (<I>Hint:</I> Represent <I>A</I> and <I>B</I> as polynomials of degree 10<I>n</I>.)<P>
<P>


<P>







<h1><a name="0991_0001">32.2 The DFT and FFT<a name="0991_0001"></h1><P>
In Section 32.1, we claimed that if we use complex roots of unity, we can evaluate and interpolate in <img src="../images/bound.gif">(<I>n </I>lg <I>n</I>) time. In this section, we define complex roots of unity and study their properties, define the DFT, and then show how the FFT computes the DFT and its inverse in just <img src="../images/bound.gif">(<I>n </I>lg <I>n</I>) time.<P>
<img src="784_a.gif"><P>
<h4><a name="0991_0002">Figure 32.2 The values of <img src="784_b.gif"> in the complex plane, where w<SUB>8</SUB><FONT FACE="Times New Roman" SIZE=2> = e<SUP>2</SUP><img src="../images/piuc.gif">i/8 <SUP><FONT FACE="Times New Roman" SIZE=2>is the principal 8th root of unity.<a name="0991_0002"></FONT></FONT></sub></sup></h4><P>





<h2>Complex roots of unity</h2><P>
<a name="0992_1b10"><a name="0992_1b11">A <I><B>complex nth root of unity</I></B> is a complex number <I>w</I> such that<P>
<pre><I>w<SUP>n</I></SUP> = 1 .</sub></sup></pre><P>
There are exactly <I>n</I> complex <I>n</I>th roots of unity; these are <I>e</I><SUP>2</SUP><img src="../images/piuc.gif"><I>ik/n</I> for <I>k</I> = 0, 1, . . . , <I>n</I> - 1. To interpret this formula, we use the definition of the exponential of a complex number:<P>
<pre><I>e<SUP>iu</I></SUP> = cos(<I>u</I>) + <I>i</I> sin(<I>u</I>).</sub></sup></pre><P>
<a name="0992_1b12">Figure 32.2 shows that the <I>n </I>complex roots of unity are equally spaced around the circle of unit radius centered at the origin of the complex plane. The value<P>
<pre><I>w<SUB>n</I></SUB> = <I>e</I><SUP>2</SUP><img src="../images/piuc.gif"><I>i</I>/<I>n</I></sub></sup></pre><P>
<h4><a name="0992_1b16">(32.6)<a name="0992_1b16"></sub></sup></h4><P>
is called <I><B>the principal nth root of unity</I></B>; all of the other complex <I>n</I>th roots of unity are powers of <I>w<SUB>n</SUB>.</I><P>
The <I>n</I> complex <I>n</I>th roots of unity,<P>