多项式乘法
记多项式 F(x),G(x) 进行乘法,得到一个新的多项式 H(x) ,记 CPoly(i) 表示 Poly 这个多项式 xi 这一项的系数是多少。
那么
CH(i)=∑ij=0CF(j)CG(i−j)
我们也称这种形式的乘法为卷积。
朴素地计算
H(x)
的时间复杂度显然是
O(n2)
的。
但是我们可以通过点值和插值运算,得到新多项式的值的信息,再得到新的多项式。
点值与插值
点值运算
一个
n
次多项式
因为
n+1
个点值对确定的
n
次多项式是唯一的。
证明
记系数列向量AT=(c0,c1,⋯,cn) , 行向量 B=(f(x0),f(x1),⋯,f(xn))
那么 B=VnA ,其中 Vn 是 n 阶范德蒙德矩阵。
而范德蒙德矩阵可逆的充要条件是xi 互不相同。那么这里就可以求出唯一的逆矩阵 V−1n ,以及唯一的系数列向量 A=BV−1n插值运算
插值运算是点值运算的逆运算。
用矩阵形式来表达则是 A=BV−1n应用
假如我们可以通过某种手段获得两个多项式的点值表达,则他们的积的多项式的点值表达则可以 O(n) 的求出来。假如我们可以通过某种手段实现插值运算,则多项式问题就迎刃而解了。
注意这里假如两个多项式的次数为 n ,那么我们需要求
2(n+1) 个点值对。因为乘出来的多项式次数是 2n 的。
N次复数根
定义
主 n 次单位复数根是满足
wn=1 的复数 w ,记为wn
n 次复数根恰好有n 个,从图像上看它们平分了复平面上的单位圆。
数值上这些根为
wkn=e2π⋅kn
其中 k=0,1,⋯,n−1
特别的, w0n=1欧拉幅角公式
eiφ=cosφ+isinφ
证明
ex=1+x+x22!+x33!⋯
sinx=x−x33!+x55!⋯
cosx=1−x22!+x44!⋯eiφ=1+iφ−φ22!−iφ33!+φ44!+⋯
cosφ=1−φ22!+φ44!⋯
isinφ=iφ−iφ33!+iφ55!⋯
eiφ=cosφ+isinφ性质
以下的性质中可以将 n 视作
2 的幂,因为实际操作时我们也是这么做的。复数的乘法对应旋转,又由欧拉幅角公式可以看出这些复数的模长都为 1 ,于是就有
wkn=wk−1nwn=(wn)k 因为 wnn=w0n=1 ,故
- wpn=wpmodnn
- wanwbn=w(a+b)modnn
消去引理: wdkdn=wkn
证明
wdkdn=e2πdkdn=e2πkn=wnk推论:
- wn2n=w2=−1
- wn2+kn=−wkn
- (wkn)2=w2kn=wkn2
折半引理
n 次复数根的平方组成的集合恰好是
n2 次复数根组成的集合证明
wkn2=(wkn)2=(wk+n2n)2求和引理
∑n−1j=0(wkn)j=0
证明
∑n−1j=0(wkn)j=0=1−(wkn)n1−wkn=1−wknn1−wkn=1−1k1−wkn=0注意这里 k 不能是
n 的倍数,不然分母就为 0 了。
计算点值表达
如何快速地计算点值表达呢?
朴素地计算点值显然是
O(n2) 的。
但是我们可以巧妙地选取复数根作自变量的值,加速运算。核心思想
我们选取 n 次复数根作为自变量的值来求点值表达。
那么计算A(wkn)=c0+c1wkn+⋯+cn−1(wkn)n−1
令 A0(x)=c0+c2x2+⋯+c2kxk
A1(x)=c1+c3x+⋯+c2k+1xk
那么 A(wkn)=A0((wkn)2)+wknA1((wkn)2)
又根据折半引理,要求的自变量数减少了一半,但是总的多项式项数不变。直到要求的自变量数只有一个,也就是 w01=1 的时候,每个多项式的项数也只有一项常数项了,于是我们可以直接得到。根据上面的式子来递归计算就可以了。
具体实现
由于在实际实现程序的时候不能对整个数组进行复制,我们要考虑一下如何使用一些技巧来避免进行一整个数组的复制。
首先假如是递归地实现,大概是以下的流程。
- 如果多项式只剩下一项了,那么返回常数项。否则
- 将多项式分成奇数项组成的多项式和偶数项组成的多项式,以当前所有自变量的平方作为自变量递归计算两个多项式的值
- 对所有的自变量合并两个分多项式所得的结果
注意到在每一层的自变量都是 wkn 的形式,其中 k 取遍
0 到 n−1 。那么我们就没必要再开多余的空间来记录有哪些自变量。我们只需要简单地知道多项式的项数 n 就可以得到所有的自变量了。于是递归版本大概长这个样子
- 如果多项式只剩下一项了,那么返回常数项。否则
- 将多项式分成奇数项组成的多项式和偶数项组成的多项式,以当前单位复数根的平方传递递归计算两个多项式的值
- 对复数根的幂对应合并两个分多项式所得的结果
但是这里还是需要多开空间来存储递归计算的结果,这一部分的处理也是这个算法比较巧妙的地方。
假如我们可以按照某种顺序来储存,并用某种顺序计算,满足
- 寻址快速(O(1))
- 线性空间
- 无后效性那么一个可行的算法就出来了。
一个简单而又有效的方式,是按照它编号的二进制翻转后的顺序进行储存系数和答案。(有趣的是,翻转再翻转就是本身。也就是说我们进行两次这样的运算以后我们得到的是原顺序,之后会提到这个性质)
当 n=8 的时候,大概是这样的一个顺序:0,4,2,6,1,5,3,7
而此时每一层存储的答案数组的意义是不同的。假如当前正在处理第 k 层
[l,r=l+n2k−1) 这个区间。那么左半边的答案数组 fi(l≤i<mid) 代表的实际上是 [l,mid=l+n2k) 区间这个多项式,以 wimod2kn2k 为自变量的值。首先来看第一次递归。显然奇数项都在右边,偶数项都在左边。
第二次递归也是(编号整除二以后),每一次递归是将区间劈成两半,左边恰好是偶数项右边恰好是奇数项。这个结果是显然的,因为二进制翻转后最低位为 0 的编号总是排在左边,最低位为1 的编号总是排在右边,相当于把当前这一层的奇数项和偶数项分开成为各自连续的一段。
于是分开多项式的问题解决了。接下来就是转移的问题。因为 (wk2n)2=(wk+n2n)2 ,因此它们来自于上一层的转移是相同的。我们将所有的求值两两分成一组来考虑,它们是相互独立的。
A(wk2n)=A0(wkn)+wk2nA1(wkn)
A(wk+n2n)=A0(wkn)+wk+n2nA1(wkn)=A0(wkn)−wknA1(wkn)至此,层间的转移就基本上解决了。
大致的程序实现应该是这样的
- 将所有的编号翻转并建立起对应关系
- 递归求值,记当前层求值的区间为 [l,r)
- 如果项数为 1 ,将该位的值记为常数项并返回
- 否则递归两侧,得到一半的多项式的对应复数根为自变量的值后,按照上面的式子得到这一层的值。具体来说
- 枚举当前转移的是哪一个分治区间。( 若当前每个区间的长度为
2k(k>0) ,那么总共有 n2k 组,每一组内有 2k−1 对转移 )- 枚举当前分治区间是哪一个转移,转移之。
高效率的实现方法
注意不要搞混枚举的顺序,假如我们先枚举了是第几个转移,然后再枚举区间,就会使得内存的读取十分的零散,使得效率大大降低。
至此,求点值的部分已经实现完毕了。
计算插值运算
理论部分
在实际操作的时候,点值运算和差值运算实际上相差不多。它相当于是通过将点值矩阵乘上范德蒙德矩阵的逆矩阵,然后得到系数矩阵。
而重点就在于得到范德蒙德矩阵的逆矩阵,不妨如下构造一个。令 V 为
n×n 的范德蒙德矩阵,我们令 V−1i,j=1nVi,j 即可得到 V 的逆矩阵V−1 证明
主要是要证明 V−1V=E
- 当 i≠j 时, (V−1V)i,j=∑n−1k=0V−1i,kVk,j
又 Vk,j=(wjn)k,V−1i,k=1nVi,k=1n(wkn)i=1n(win)k
因此原式 =∑n−1k=0(wjn)k⋅1n(win)k=1n∑n−1k=0(wi−jn)k
因为 wi−jn≠0 ,根据求和引理,原式 =0- 当 i=j 时, 原式=1n∑n−1k=0(wi−jn)k=1n∑n−1k=01k=1
故 V−1V=E ,证毕。
实际操作
程序实现的时候,我们并不需要重新实现一段求插值的代码。因为求点值和求插值十分的相似,我们只需要修改一小部分的参数就可以了。具体来说
- 令 n 次单位复数根为
w−n1 - 重新使用上述递归过程
- 将结果的每一位除 n
就完成了。
快速傅里叶变换
上述的点值表达和插值运算综合在一起就称快速傅里叶变换(fast Fourier transfer)。
理论部分到此就结束了。
talk is cheap, show me your code.
延伸
- 模意义下的FFT
- 模数为
k2n+1 ,且是质数- 模数是质数
- 模数是一般数
- 快速威尔士变换(fast Walsh transform),即FFT的布尔运算形式