今天一天尽享多项式!Í dag er「多项式」æfa! 今日は一日「多项式」三昧!

题目是在玩梗。甩你个冰岛语词典自己查吧。这个句子大概这么发音: ɪː tɑɣ ɪrʱ 多项式  ævɑ(肯定不太准)

拉格朗日插值法

插值是什么?

差不多就是以一个函数 f(x) 已知的若干点值代入反解,做出一个特定的函数求对于一个给定的 xf(x) 的近似值的过程。

如果这个特定函数是一个多项式,就称为插值多项式

其中一种方法是 O(n^3) 类的高消。

拉格郎日插值法在干什么?

构造一个对于所有情况都通用的插值基函数,也叫做拉格朗日基本多项式 , 然后可以一步步由它求出插值多项式。

为什么叫做基函数?因为最后的插值多项式是由这个函数线性组合得到的。

这个基函数要干什么?

我们通过 n+1 个点,构造若干个多项式(插值基函数)做线性组合,其中对于第 i 个多项式在 x_i 处的取值为 y_i ,在其他点处取值为 0

把这个多项式叫做 l_i(x) ,那么利用这个性质,最后的插值多项式就可以表示成

l(x) = \sum_{i=0}^n{y_il_i\left( x \right)}

尝试把给定的 n+1个点值代入,发现插值多项式对于所有的点值都成立。于是就偷税地插值完了。

构造

这里利用了一个性质: 0 乘任何数都为 01 乘任何数都不变。

所以按照上文, l_i(x) 就要满足 l_i(x_i)=1 \,,\, \forall j \in [0,n],j \ne i \,,\, l_i(x_j)=0

那么这个构造非常的暴力:

l_i\left( x \right) =\prod_{j=\text{0 }j\ne i}^n{\frac{x-x_j}{x_i-x_j}}

代入各个 x_i 算一下就可以发现满足该性质。

于是我们用 O(n^2) 的时间复杂度求出了拉格日插值多项式。就是这个东西:

l\left( x \right) =\sum_{i=0}^n{\left( y_i\prod_{j=\text{0,}j\ne i}^n{\frac{x-x_j}{x_i-x_j}} \right)}

如果给定的点值是 [0,n] 内的话,那么有方法可以优化到 O(n) 的复杂度。

拉格日插值的线代理解

显然插值基函数是线性无关的,所有的点都可以被线性表出。

所以我们只是构造了一组较好计算的基,使得可以代入 x 来确定坐标 (x,y) 而已。

为什么这样是对的?

我们回到朴素的插值法。这玩意其实就是解线性方程组:

\left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array} \right) =\left( \begin{matrix} 1& x_1& \cdots& x_{1}^{n-1} \\ 1& x_2& \cdots& x_{2}^{n-1} \\ \vdots& \vdots& \ddots& \vdots \\ 1& x_n& \cdots& x_{n}^{n-1} \end{matrix} \right) \left( \begin{array}{c} a_0 \\ a_1 \\ \vdots \\ a_n \end{array} \right)

简写一下就是

\vec{y}=A\vec{a}

如果 |A| \ne 0 那方程组就有唯一解。 观察到 |A| 其实是范德蒙德行列式。

|A|=\prod_{1\leqslant i<j\leqslant n}{\left( x_j-x_i \right)}

如果 x_i 互不相等,那么 |A| \ne 0

所以拉格朗日插值得到的多项式曲线是唯一的。

代码怎么写?

慢速傅里叶变换 (Foolish Fourier Transformation)

FFT 其实就是用分治法加速 DFT 和 IDFT 的过程。

我对大部分博文都不太满意的一点是,思路偏重点不满意。在大多数情况下 (除非实在找不到办法) ,我更偏向的是探究型的讲解而非灌输型的讲解。我觉得有机会去搞懂这里的方法、思路是如何发现的,为什么要使用这样的思路,这一点是最重要的。

不过话说回来, FFT 这类优化思路过于偏门,在 OI 中搞懂了 FFT 是怎么来的并没有助于你解题。

————第 1 节————

先来考虑我们的任务:O(n \log n) 的多项式乘法。

首先第一点是我们发现系数表示法的复杂度下界是 O(n^2) ,因为对于第一个多项式的每一个系数,都要和第二个多项式的 n 个系数做一次乘法。

所谓系数表示法就是我们惯常的形式:

A\left( x \right) =a_0+a_1x+\cdots +a_nx^n=\left( \begin{matrix} a_0& a_1& \cdots& a_n \end{matrix} \right) \left( \begin{array}{c} 1 \\ x \\ \cdots \\ x^n \end{array} \right)

那么我们有一种优化是将多项式由系数表示法转为点值表示法。对于一个 n 次多项式,依据上文,可以用 n+1 个互不相等的点唯一确定:

f\left( x \right) =a_0+a_1x+a_2x^2\,\,\Leftrightarrow \,\,f\left( x \right) =\left\{ \left( x_1,y_1 \right) ,\left( x_2,y_2 \right) ,\left( x_3,y_3 \right) \right\}

然后就可以 O(n) 搞出乘积,如下:

f\left( x \right) =\left\{ \left( x_i,A_i \right) \right\}

g\left( x \right) =\left\{ \left( x_i,B_i \right) \right\}

\left( f\cdot g \right) \left( x \right) =\left\{ \left( x_i,A_iB_i \right) \right\}

————第 2 节————

但是问题是对于这个多项式,求 n 个点值的复杂度还是 O(n^2) 的,非常棘手。

这里我们就开始玩 tricky 的变换了。首先我们假设这个多项式是 2^m-1 次的(这样我们就需要求 2^m 个点值)

n 次多项式为(这里我比较喜欢用矩阵的写法,系数和变量分开,看着清楚)

A\left( x \right) =a_0+a_1x+\cdots +a_nx^n=\left( \begin{matrix} a_0& a_1& \cdots& a_n \end{matrix} \right) \left( \begin{array}{c} 1 \\ x \\ \cdots \\ x^n \end{array} \right)

我们可以按奇偶性分组:

A\left( x \right) =\left( a_0+a_2x^2+a_4x^4 \right) +\left( a_1x+a_3x^3+a_5x^5 \right)

=\left( a_0+a_2x^2+a_4x^4 \right) +x\left( a_1+a_3x^2+a_5x^4 \right)

然后如果我们设

A_0\left( x \right) =\left( \begin{matrix} a_0& a_2& a_4& \cdots \end{matrix} \right) \left( \begin{array}{c} 1 \\ x \\ x^2 \\ \cdots \end{array} \right)

A_1\left( x \right) =\left( \begin{matrix} a_1& a_3& a_5& \cdots \end{matrix} \right) \left( \begin{array}{c} 1 \\ x \\ x^2 \\ \cdots \end{array} \right)

那么我们就有

A\left( x \right) =A_0\left( x^2 \right) +xA_1 \left( x^2 \right)

这挺好的。我们可以递归式地计算每个点值。不过并不会优化什么运算。为什么?

注意到

x_1^2 , x_2^2 , \cdots , x_n^2

仍极有可能组成 n 个不同的数字。对于每个 A(x_i) ,我们都要分别递归计算 A_0(x_i)A_1(x_i) ,这样的复杂度是高于 O(n \log n)  的 O(n^2) ,唉。

那么我们的任务只剩下找一组特定的 \{x\} 使得它不组成 n 个不同的数字,而是 \frac{n}{2} 个,那么递归计算 A(x) 时规模即可减半。

————第 3 节————

比如说我让这组 \{x\}

\left\{ x_0,\cdots ,x_{\frac{n}{2}-1},-x_0,\cdots ,-x_{\frac{n}{2}-1} \right\}

那么平方之后不就减半了么?

但是这种关系不会传递下去。平方后得到的数字全是不小于 0 的数字。再次递归又回到了原来的形式。很尴尬。

不过这启迪我们。取相反数然后平方,可以把问题规模减半。但再一次递归就失效了。是否存在不失效的取值呢?

如果对于一个有偶数个元素的数集平方后得到的新集合,大小能够减半,且新集合也满足性质,这个性质就是折半性质。如果我们可以快速的找到一个满足折半性质的自变量 x 的取值集合,分治就是可行的。

是可以找到的。这个东西就是 n 次单位复根。

n 次单位复根是满足 \omega^n = 1 的复数 \omega , 其中由复数的运算法则(辐角相乘,模相加)可以很简单地得出 n次单位根有 n 个这个结论。

而又因为复数 \omega^n 在复数平面上的模都是 1 ,所以相乘之后还会是 1 ,那么所有的 w_i (1 \leqslant i <n) 就会均匀分布在单位圆上,类似当 n = 16 时它是这样的:

这里首先我列出几条性质:

\left( \cos \theta +i\sin \theta \right) ^k=\cos \left( k\theta \right) +i\sin \left( k\theta \right)

这个性质可以归纳法证出来。

还有很重要的两个性质

考虑欧拉公式

e^{2\pi i}=1=\omega ^n

易得

\omega =e^{\frac{2\pi i}{n}}

这时的单位根称之为主次单位根,记作

\omega_n =e^{\frac{2\pi i}{n}}

那么对于其他的单位根,记作

\omega _{n}^{k}=e^{\frac{2\pi ik}{n}},0\leqslant k<n

显然都是主次单位根整数幂

消去引理

\omega_{dn}^{dk}=\omega_{n}^{k}\,\,\left( d>0n\geqslant \text{0,}k\geqslant 0 \right)

证明:代入欧拉公式即可。

折半引理

\left( \omega _{n}^{k} \right) ^2=\omega_{\tfrac{n}{2}}^{k}

证明:代入消去引理即可。

于是我们证明 n 次单位根集合的折半性质

\left( \omega _{n}^{k+\frac{n}{2}} \right) ^2=\omega _{n}^{2k+n}=\omega _{n}^{2k}\times \omega _{n}^{n}=\omega _{n}^{2k}=\left( \omega _{n}^{k} \right) ^2

那么大功告成 !  \left\{ \left( \omega _{n}^{i} \right) ^2 \right\} 就等效于 \left\{ \omega _{\frac{n}{2}}^{0},\omega _{\frac{n}{2}}^{1},\cdots ,\omega _{\frac{n}{2}}^{\frac{n}{2}-1},\omega _{\frac{n}{2}}^{0},\cdots ,\omega _{\frac{n}{2}}^{\frac{n}{2}-1} \right\}

 x 的取值集合取单位复数根 , 不但满足折半的性质,而且还有一定的规律性,与原集合保持一致。这意味着我们只需要计算取 \left\{ \omega _{\frac{n}{2}}^{0},\omega _{\frac{n}{2}}^{1},\cdots ,\omega _{\frac{n}{2}}^{\frac{n}{2}-1} \right\} A_0A_1 即可。问题规模成功减半,复杂度降为 O(n \log n)

具体计算怎么计算呢?

A_{k}^{0}=A_0\left( \omega _{\frac{n}{2}}^{k} \right) =A_0\left( \left( \omega _{n}^{k} \right) ^2 \right)

A_{k}^{1}=A_1\left( \omega _{\frac{n}{2}}^{k} \right) =A_1\left( \left( \omega _{n}^{k} \right) ^2 \right)

那么我们要求的就是所有

A\left( \omega _{n}^{k} \right) =A_0\left( \omega _{n}^{2k} \right) +\omega _{n}^{k}A_1\left( \omega _{n}^{2k} \right) =A_0\left( \omega _{\frac{n}{2}}^{k} \right) +\omega _{n}^{k}A_1\left( \omega _{\frac{n}{2}}^{k} \right) =A_{k\,\text{mod }\frac{n}{2}}^{0}+\omega _{n}^{k}A_{k\,\,\text{mod }\frac{n}{2}}^{1}

另外我们要折半,这样才能 O(n \log n) 。假设上面式子的 k < \frac{n}{2} ,注意到

A\left( \omega _{n}^{k} \right) =A_{k}^{0}+\omega _{n}^{k}A_{k}^{1}

A\left( \omega _{n}^{k+\frac{n}{2}} \right) =A_0\left( \omega _{n}^{2k+n} \right) +\omega _{n}^{k+\frac{n}{2}}A_1\left( \omega _{n}^{2k+n} \right) =A_{k}^{0}-\omega _{n}^{k}A_{k}^{1}

后一个式子用到的几步变换无非就是把单位根的指数上拆了罢了。首先是指数是加了 \frac{n}{2} 相当于转了半圈,变为相反数。其次是 \omega_n^{2k+n} = \omega_n^n \times \omega_n^{2k}

然后就做完了 DFT 。再来理一下大致思路

想要转换点值表达式乘起来

→ 变换一下发现可以分治

→ 要算的 x 值太多,复杂度不对

→ 引入单位复根,折半算值

→ 求出点值表达式

————第 4 节————

我们来从线代角度考虑 DFT 究竟在做什么。

无非就是取了单位复根矩阵去算点值

\left( \begin{array}{c} A_0 \\ A_1 \\ \vdots \\ A_{n-1} \end{array} \right) =\left( \begin{matrix} \left( \omega _{n}^{0} \right) ^0& \left( \omega _{n}^{0} \right) ^1& \cdots& \left( \omega_{n}^{0} \right) ^{n-1} \\ \left( \omega _{n}^{1} \right) ^0& \left( \omega _{n}^{1} \right) ^1& \cdots& \left( \omega _{n}^{1} \right) ^{n-1} \\ \vdots& \vdots& \ddots& \vdots \\ \left( \omega _{n}^{n-1} \right) ^0& \left( \omega _{n}^{n-1} \right) ^1& \cdots& \left( \omega _{n}^{n-1} \right) ^{n-1} \end{matrix} \right) \left( \begin{array}{c} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{array} \right)

因为方程显然可解,那么这个矩阵存在逆矩阵。 所谓 IDFT 也即从点值表达还原为系数表达就是在点值表达矩阵上乘这个大矩阵的逆矩阵。

由于某些包括证明过长和知道也没啥用的缘由,这里不再给出证明,只给出性质:这个矩阵的逆矩阵为

E=\frac{1}{n}\left( \begin{matrix} \left( \omega _{n}^{-0} \right) ^0& \left( \omega _{n}^{-0} \right) ^1& \cdots& \left( \omega _{n}^{-0} \right) ^{n-1} \\ \left( \omega _{n}^{-1} \right) ^0& \left( \omega _{n}^{-1} \right) ^1& \cdots& \left( \omega _{n}^{-1} \right) ^{n-1} \\ \vdots& \vdots& \ddots& \vdots \\ \left( \omega _{n}^{-\left( n-1 \right)} \right) ^0& \left( \omega _{n}^{-\left( n-1 \right)} \right) ^1& \cdots& \left( \omega _{n}^{-\left( n-1 \right)} \right) ^{n-1} \end{matrix} \right)

然后这样 IDFT 实现上就是让 DFT 时的每个 \omega_n^k 变为 \omega_n^{-k} ,然后再乘 \frac{1}{n} 即可。

————第 5 节————

这里并不打算讲那个没啥思维含量还和 FFT 本身不太搭边的优化。不过 FFT 代码还是要搬上来的。完整版的话去翻别的文吧。

呆子数论变换 ( Nerd Theory Transformation )

 

 

 

 

 

 

 

 

 

 

 

 

 

“今天一天尽享多项式!Í dag er「多项式」æfa! 今日は一日「多项式」三昧!”的2个回复

GNAQ进行回复 取消回复

邮箱地址不会被公开。 必填项已用*标注

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据