FFT
相信大家FFT已经掌握的很熟练了。
参考这一篇博客:浅谈算法——从多项式乘法到FFT
常数优化1——IDFT
若已知多项式A(x)的点值表示
<(ωn0,A(ωn0)),(ωn1,A(ωn1)),⋯,(ωnn−1,A(ωnn−1))>
则多项式A(x)满足
[xrev(k)]A(x)=n1i=0∑n−1ωnkiA(ωni)
其中
rev(k)={0n−k−1k=0k̸=0
就是说:IDFT只需要调用DFT函数然后std::reverse(a+1,a+n+1)
,最后同时除以n即可。
证明:假设
A(x)=i=0∑n−1aixi
那么
[xrev(k)]A(x)=n1i=0∑n−1ωnkij=0∑n−1ajωnij
交换两个求和符号
[xrev(k)]A(x)=n1i=0∑n−1aij=0∑n−1ωnijωnkj
显然
[xrev(k)]A(x)=n1i=0∑n−1aij=0∑n−1ωn(i+k)j
此时,若i=rev(k),则i+k=0(modn),此时
ωn(i+k)j=1
即
j=0∑n−1ωn(i+k)j=n
否则由等比数列求和公式
j=0∑n−1ωn(i+k)j=1−ωni+k1−ωn(i+k)n=0
因此
[xrev(k)]A(x)=n1arev(k)n=arev(k)
常数优化2
假设有两个实多项式A(x),B(x),现在要求出他们两个的FFT后的结果。若
P(x)=A(x)+iB(x)
则
A(ωnk)=2P(ωnk)+conj(P(ωnrev(k)))B(ωnk)=i2P(ωnk)−conj(P(ωnrev(k)))
其中i=−1,conj(x)为x的共轭虚数。
证明:容易发现,我们只需要证明
conj(P(ωnrev(k)))=A(ωnk)−iB(ωnk)
即可得证。
假设
A(x)=i=0∑n−1aixiB(x)=i=0∑n−1bixi
那么
A(ωnk)−iB(ωnk)=j=0∑n−1ajωnkj−ij=0∑n−1bjωnkj
由于
ωnk=cosn2πk+isinn2πk
带入得
A(ωnk)−iB(ωnk)=j=0∑n−1aj(cosn2πkj+isinn2πkj)+j=0∑n−1bi(−icosn2πkj+sinn2πkj)
而由于ωnrev(k)=ωn−k
conj(P(ωnrev(k)))=conj(A(ωn−k)+iB(ωn−k))
即
conj(P(ωnrev(k)))=conj(j=0∑n−1aj(cos(−n2πkj)+isin(−n2πkj))+j=0∑n−1bj(icos(−n2πkj)−sin(−n2πkj)))
由于sin(−x)=−sinx,cos(−x)=cosx
conj(P(ωnrev(k)))=conj(j=0∑n−1aj(cosn2πkj−isinn2πkj)+j=0∑n−1bj(icosn2πkj+sinn2πkj))
由于aj,bj是实数,因此
conj(P(ωnrev(k)))=A(ωnk)−iB(ωnk)
因此,我们将两个多项式的FFT变成了一个多项式的FFT。
MTT
对于两个多项式
A(x),B(x)
求多项式乘法,mod109+7之类的质数。
不能NTT,FFT的精度误差太大,怎么办?
假设这个质数为p,令p′=⌊p⌋
A(x)=p′QA(x)+RA(x)B(x)=p′QB(x)+RB(x)
其中RA(x)和RB(x)的每一项的系数都<p′,那么
A(x)B(x)=p′2QA(x)QB(x)+p′(QA(x)RB(x)+QB(x)RA(x))+RA(x)RB(x)
这样能大大减小精度误差。
常数优化
使用上面FFT部分中提到的常数优化2即可。
多项式求逆
A(x)F(x)=1modxn
其中n是多项式A(x)的度数。
假设我们已经求出了
A(x)F0(x)=1modx⌈n/2⌉
上面两式相减
F(x)−F0(x)=0modx⌈n/2⌉
两边平方
F2(x)−2F(x)F0(x)+F02(x)=0modxn
左右同时乘以A(x)
F(x)−2F0(x)+A(x)F02(x)=0modxn
即
F(x)=2F0(x)−A(x)F02(x)modxn
常数优化
对于MTT在多项式求逆上的常数优化,参见这篇文章:一个多项式求逆的卡常技巧
多项式取对数
F(x)=lnA(x)modxn
考虑对两边同时求导
F′(x)=A−1(x)A′(x)modxn
因此
F(x)=∫A−1(x)A′(x)dx
牛顿迭代
对于一个多项式F(x),满足
G(F(x))=0modxn
G(x)已知,现在要求F(x)。
假设已经求出了
G(F0(x))=0modx⌈n/2⌉
考虑对G(F(x))在F0(x)处展开
G(F(x))=G(F0(x))+1!G′(F0(x))(F(x)−F0(x))+2!G′′(F0(x))(F(x)−F0(x))2+⋯
注意这里的G′(F(x))是G(F(x))对F(x)取导数的结果。
由于F(x)和F0(x)的最后⌈2n⌉项相同,因此
(F(x)−F0(x))k=0modxn
在k≥2时都成立。
因此
G(F(x))=G(F0(x))+G′(F0(x))(F(x)−F0(x))modxn
由于
G(F0(x))=0modxn
因此
F(x)=F0(x)−G′(F0(x))G(F0(x))modxn
多项式求指数
F(x)=eA(x)modxn
两边同时取对数得到
lnF(x)=A(x)modxn
令
G(F(x))=lnF(x)−A(x)=0modxn
容易发现
G′(F(x))=F−1(x)modxn
套牛顿迭代公式,假设我们已经求出了
lnF0(x)−A(x)=0modx⌈n/2⌉
则
F(x)=F0(x)−F0(x)(lnF0(x)−A(x))modxn
更优雅的式子是
F(x)=F0(x)(1−lnF0(x)+A(x))modxn
多项式开方
F2(x)=A(x)modxn
考虑
G(F(x))=F2(x)−A(x)=0modxn
容易发现
G′(F(x))=2F(x)modxn
套牛顿迭代公式,假设我们已经求出了
F02(x)−A(x)=0modx⌈n/2⌉
即有
F(x)=F0(x)−2−1F0−1(x)(F02(x)−A(x))modxn
更优雅的式子是
F(x)=2−1(F0(x)+A(x)F0−1(x))modxn
拉格朗日反演
令
G(F(x))=x
已知G(x)要求F(x)。其中,G(x)无常数项。
经过一些奥妙重重的变换可以得到
[xn]F(x)=n1[xn−1]G′n(x)1
其中
G′(x)=xG(x)
为啥?我也不知道……
UPD:证明