PolyEval.muse

#title 一元多项式求值和插值

<contents>

; * 求幂算法

; 我们在因子分解一章中将要用到求幂算法,这里先稍加介绍.这里要提到的求幂算法,实际上是一种二进算法,我们可以将求幂的过程描述为下面的算法.见<cite>taocp2</cite>4.6.3节.

; <algorithm name="二进求幂算法1">
; 设$n$为一正整数,我们要求$x^n$,可将$n$表示为其二进制形式,并将二进制序列中的1用SX代替,0用S代替,并去掉第一个SX,然后自左至右,遇到S进行平方(square),遇到X则乘以x,最后得的结果为$x^n$.
; </algorithm>

; 例如$n=13$时,其二进制表示为$1101$,对应的序列为$SXSSX$,于是结果为$(((x^2)x)^2)^2x=x^{13}$.

; 该算法需要从高位到低位依次扫描$n$的位,对于计算机来说,更方便的应该是从低位到高位,因此还有下面的算法:

; <algorithm  name="二进求幂算法2">

;  1. 赋初值$y=1$,$z=x$,$m=n$,

;  2. 取$m$的最低位,若其为1,则$y=yz$,$z=z^2$,否则$z=z^2$,$y$不变,

;  3. $m$右移一位,

;  4. 若$m=0$,则输出$y$,否则转2步.
; </algorithm>

; <proof name="算法有效性">

; 我们只需要归纳地证明每一步之后均有$x^n=yz^m$.首先在初始条件下该式成立,每次执行第2,3步之后,若$m$最低位为1,设新的$y,z,m$用$y',z',m'$表示,则$y'=yz,z'=z^2,m'=(m-1)/2$,此时仍有$y'z'^{m'}=yz(z^2)^{(m-1)/2}=yz^m$.

; 同样地,若$m$最低位为0,则$y'=y,z'=z^2,m'=m/2$,$y'z'^{m'}=y(z^2)^{m/2}=yz^m$.
; </proof>

; 类似地,我们还有$m$进方法,主要利用$n$的$m$进制表示来减少求幂过程中的乘法运算.

; 在求幂过程中,通过对$n$的因子分解,也能减少乘法运算.例如对于$n=15$的情况, 如果用二进方法,需要做6次乘法运算,而用因子分解法,我们可以先利用二进方法,用2次乘法算出$y=x^3$,再用3次乘法算出$y^5$,总共只需5次乘法.

; 下面再介绍一种比二进算法和因子分解算法稍微优越的幂树算法.幂树的构造方法是:首先在第一层置1,当置完$k$层后,从左到右依次取第$k$层每个节点$n$,在其下附加节点$n+1,n+a_1,\ldots,n+a_{k-1}=2n$,其中$1,a_1,a_2,\ldots,a_{k-1}$是树根到$n$的到$n$的通路上的点,但已出现的节点不予添加.下图显示了前4层幂树:

; [[images/powertree.png][前4层幂树]]

; 有了幂树之后,当我们要计算某个$n$次幂时,只需要在该树上找到$n$,依次计算从根节点到$n$的路径上的幂次即可.


本章主要介绍一元多项式求值和插值的快速算法<cite>mca</cite>，它们本身都不是很复杂，因此不打算用太多的篇幅介绍。关于这一方面,有兴趣的读者也可以参阅<cite>cmf72</cite>,<cite>eh72</cite>.

由于我们从本章开始将进入多项式的相关算法，首先对一些通用的多项式记号作一简单说明。设有多项式
$$f=\sum_{i=0}^na_nx^n=a_0+a_1x+\cdots+a_nx^n,$$
则称$n$为$f$的次数，记为$\deg f$,称$a_n$为$f$的<index>领项系数</index>(Leading Coefficient)，记为$\lc (f)$，称$x^n$为$f$的<index>领项单项式</index>(Leading Monomial),记为$\lm{f}$,称$a_nx^n$为$f$的<index>领项</index>(Leading Term)，记为$\lt{f}$,后文中也会提到首项,即领项.定义$f$的首一化多项式为
$$\mathrm{monic}(f)=\frac{f}{\plc{f}}.$$

* 求值算法

提起一元多项式的求值算法,我们都会想起最著名的<index>Horner规则</index>,即对于$f(x)=\sum_{0\le i<n}f_ix^i$,利用下式来求其值:
$$f(x)=(\cdots(f_{n-1}x+f_{n-2})x+\cdots+f_1)x+f_0.$$
在该算法中,需要计算总共$n-1$次乘法和$n-1$次加法.如果要同时计算多项式在$n$个不同点处的值,则需要$O(n^2)$的计算量.

对于一般多项式的单点求值问题来说,Horner规则已经是最优的了.Pan<cite>pan66</cite>于1966年证明了Horner规则使用乘法的次数是最少的,即$n$次多项式的求值至少要$n$次乘法.

但是对于多点求值来说,Horner规则就不一定是最优的了.下面给出一种快速多点求值算法,其计算复杂度为$O(M(n)\log n)$,其中$M(n)$表示$n$次多项式乘法计算的复杂度(参考<cite>mca</cite>封三说明).为了说明算法,我们取$n=2^k$是2的一整数次幂,需要求值的点为$u_0,u_1,\ldots,u_{n-1}$,并且令$m_j=x-u_j(j=0,\ldots,n-1)$.下面构造一棵完全二叉树,以$M_{i,j}$表示从叶($i=0$)往上数第$i$层,从该层左往右数第$j$个结点,结点值为
$$M_{i,j}=\prod_{j\cdot 2^i\le l<(j+1)\cdot 2^i}m_l,$$
图<literal>\ref{fig:interpolationtree}</literal>表示其结构.

<latex>
\begin{figure}[h]
\centering\includegraphics[scale=0.75]{images/interpolationtree}
\caption{$M_{i, j}$~二叉树示意图}
\label{fig:interpolationtree}
\end{figure}
</latex>

其中每个叶结点都是一次式$M_{0,j}=m_j$.算法##al:bitree 给出了上面树的构造算法.

<algorithm name="$M_{i,j}$二叉树构造算法" label="al:bitree">

 1. 对$0\le j<n$,令$M_{0,j}=m_j$,

 2. 对$i$从$1$循环到$k$,做下面第3步,

 3. 对$j$从$0$循环到$2^{k-i}$,$M_{i,j}=M_{i-1,2j}M_{i-1,2j+1}$.
</algorithm>

利用上面的构造的二叉树,我们可以实现快速求值算法.

<algorithm  name="快速求值算法">

输入:多项式$f$和$n$(满足$\deg f<n$)个点$u_0,\ldots,u_{n-1}$,

输出:$f(u_0),\ldots,f(u_{n-1})$.

 1. 若$n=1$则输出$f$(此时$f$为常数), 

 2. $r_0=f\bmod M_{k-1,0}$,$r_1=f\bmod M_{k-1,1}$,

 3. 递归调用本算法求$r_0(u_0),\ldots,r_0(u_{n/2-1})$,

 4. 递归调用本算法求$r_1(u_{n/2}),\ldots,r_1(u_{n-1})$,

 5. 输出上面两步求出的$n$个值.
</algorithm>

我们举例来说明上述算法.

<problem label="ex:eval">
求多项式$f(x)=x^3+2x^2+3x+4$在$\{u_0,u_1,u_2,u_3\}=\{1,2,3,4\}$四个点处的值.
</problem>

<solution>
首先由$f\bmod M_{1,0}=f \bmod (x-1)(x-2)=16x-6$和$f\bmod M_{1,1}=f \bmod (x-3)(x-4)=54x-104$可以将问题化为求$16x-6$在$1,2$处的值和$54x-104$在$3,4$处的值.再做一次模运算（模一次多项式实际上相当于求值）可以求出四个点处的值分别为$10,26,58,112$.
</solution>

<remark>
以上讨论的都是$n$为2的整数次幂的情况,而一般情况下这是不能满足的,这时,我们可以采取添一些项,例如在多点求值时添一些不同的求值点以凑成2的整数次幂情形,或者不必得到完全二叉树,在递归建立$M_{i,j}$的树或求值时每次将待处理的点分成数目大致相同的两组进行处理.下一小节的插值算法当$n$不是2的整数次幂时,也可同样处理.
</remark>

* 插值算法

多项式的插值是指已知多项式在$n$个不同点处的值,求出此多项式,当然是在模$x^n$的意义下,即次数不超过$n$时，此多项式是唯一的.首先我们想到了著名的<index>Lagrange插值</index>算法,设已知$n$个点$u_0,\ldots,u_{n-1}$和$n$次多项式$f$在这些点上的值$v_0,\ldots,v_{n-1}$,记$m_j=x-u_j$,$m=\prod_{0\le j<n}m_j$,$s_i=\prod_{j\neq i,0\le j<n}\displaystyle\frac{1}{u_i-u_j}$,则Lagrange插值的结果可以表示为:
$$f=\sum_{0\le i<n}\frac{v_is_im}{m_i}.$$
这是一个复杂度$O(n^2)$的算法(见<cite>mca</cite>定理5.1).一元插值的方法还有Newton插值法,其复杂度也是$O(n^2)$的<cite>rz93</cite>,见多元多项式插值部分的算法##al:newtoninterpolation .

下面提出的快速插值算法,其复杂度与求值算法一样,也为$O(M(n)\log n)$.为求插值多项式,首先我们要求出$s_i$,因为
$$m'(u_i)=\sum_{0\le j<n}\left.\frac{m}{m_j}\right|_{x=u_i}=\left.\frac{m}{m_i}\right|_{x=u_i}=\frac{1}{s_i},$$
故$s_i=\displaystyle\frac{1}{m'(u_i)}$.现在令$c_i=v_is_i$,并设$n=2^k$是2的整数次幂.

<algorithm  name="快速插值算法">

输入:$u_0,\ldots,u_{n-1}$,$c_0,\ldots,c_{n-1}$,

输出:$\displaystyle\sum_{0\le i<n}\frac{c_im}{m_i}$.

 1. 若$n=1$则输出$c_0$,

 2. 递归调用本算法,输入前$n/2$个点,求出$r_0$,

 3. 递归调用本算法,输入后$n/2$个点,求出$r_1$,

 4. 输出$M_{k-1,1}r_0+M_{k-1,0}r_1$.
</algorithm>

<proof  name="算法有效性">
我们只需证明算法能返回结果$\sum_{0\le i<n}c_i\prod_{0\le j<n,j\neq i}m_j$.采用递归论证的方法,我们假设在每一次递归调用执行下一级算法时,均会得到正确的结果,即
$$r_0=\sum_{0\le i<n/2}c_i\prod_{0\le j<n/2,j\neq i}m_j,\quad r_1=\sum_{n/2\le i<n}c_i\prod_{n/2\le j<n,j\neq i}m_j.$$
那么该步算法将返回结果:
$$M_{k-1,1}r_0+M_{k-1,0}r_1=r_0\prod_{n/2\le j<n}m_j+r_1\prod_{0\le j<n/2}m_j=\sum_{0\le i<n}c_i\prod_{0\le j<n,j\neq i}m_j,$$
即若下一级递归结果正确,则该步算法也将返回正确结果.另外,在算法递归终止处,即$n=1$情况下,算法返回$c_0$,此结果显然是正确的. 算法有效性得证.
</proof>

同样地,我们给出一个具体的例子来说明本算法.
<problem>
考虑$\{u_0,u_1,u_2,u_3\}=\{1,2,3,4\}$,$\{v_0,v_1,v_2,v_3\}=\{10,26,58,112\}$,求次数小于4的多项式$f$使得$f(u_i)=v_i(1\le i\le 4)$.
</problem>

<solution>
首先$m=\prod_{0\le i<4}(x-u_i)=x^4-10x^3+35x^2-50x+24$,则$m'=4x^3-30x^2+70x-50$.于是
$$s_0^{-1}=m'(1)=-6,s_1^{-1}=m'(2)=2,s_2^{-1}=m'(3)=-2,s_3^{-1}=m'(4)=6.$$
由$c_i=v_is_i$可得
$$c_0=-\frac{5}{3},c_1=13,c_2=-29,c_3=\frac{56}{3}.$$
由于$n=4$,情况比较简单,只需递归两次,很容易得到第一级算法中的
$$r_0=-\frac{5}{3}(x-2)+13(x-1)=\frac{34}{3}x-\frac{29}{3},r_1=-29(x-4)+\frac{56}{3}(x-3)=-\frac{31}{3}x+60,$$
则最后插值结果为:
$$(x-3)(x-4)r_0+(x-1)(x-2)r_1=x^3+2x^2+3x+4.$$
可以看到,此结果与例##ex:eval 是相符合的.
</solution>