FFT

相信大家FFT已经掌握的很熟练了。

参考这一篇博客:浅谈算法——从多项式乘法到FFT

常数优化1——IDFT

若已知多项式A(x)A(x)的点值表示
<(ωn0,A(ωn0)),(ωn1,A(ωn1)), ,(ωnn1,A(ωnn1))> <(\omega_n^0,A(\omega_n^0)),(\omega_n^1,A(\omega_n^1)),\cdots,(\omega_n^{n-1},A(\omega_n^{n-1}))>
则多项式A(x)A(x)满足
[xrev(k)]A(x)=1ni=0n1ωnkiA(ωni) [x^{rev(k)}]A(x)=\frac{1}{n}\sum_{i=0}^{n-1}\omega_n^{ki}A(\omega_n^i)
其中
rev(k)={0k=0nk1k̸=0 rev(k)=\begin{cases} 0 & k=0\\ n-k-1 & k\not= 0 \end{cases}
就是说:IDFT只需要调用DFT函数然后std::reverse(a+1,a+n+1),最后同时除以nn即可。

证明:假设
A(x)=i=0n1aixi A(x)=\sum_{i=0}^{n-1}a_ix^i
那么
[xrev(k)]A(x)=1ni=0n1ωnkij=0n1ajωnij [x^{rev(k)}]A(x)=\frac{1}{n}\sum_{i=0}^{n-1}\omega_n^{ki}\sum_{j=0}^{n-1}a_j\omega_n^{ij}
交换两个求和符号
[xrev(k)]A(x)=1ni=0n1aij=0n1ωnijωnkj [x^{rev(k)}]A(x)=\frac{1}{n}\sum_{i=0}^{n-1}a_i\sum_{j=0}^{n-1}\omega_n^{ij}\omega_n^{kj}
显然
[xrev(k)]A(x)=1ni=0n1aij=0n1ωn(i+k)j [x^{rev(k)}]A(x)=\frac{1}{n}\sum_{i=0}^{n-1}a_i\sum_{j=0}^{n-1}\omega_n^{(i+k)j}
此时,若i=rev(k)i=rev(k),则i+k=0( mod n)i+k=0(\bmod n),此时
ωn(i+k)j=1 \omega_n^{(i+k)j}=1

j=0n1ωn(i+k)j=n \sum_{j=0}^{n-1}\omega_n^{(i+k)j}=n
否则由等比数列求和公式
j=0n1ωn(i+k)j=1ωn(i+k)n1ωni+k=0 \sum_{j=0}^{n-1}\omega_n^{(i+k)j}=\frac{1-\omega_n^{(i+k)n}}{1-\omega_n^{i+k}}=0
因此
[xrev(k)]A(x)=1narev(k)n=arev(k) [x^{rev(k)}]A(x)=\frac{1}{n}a_{rev(k)}n=a_{rev(k)}

常数优化2

假设有两个实多项式A(x),B(x)A(x),B(x),现在要求出他们两个的FFT后的结果。若
P(x)=A(x)+iB(x) P(x)=A(x)+iB(x)

A(ωnk)=P(ωnk)+conj(P(ωnrev(k)))2B(ωnk)=iP(ωnk)conj(P(ωnrev(k)))2 A(\omega_n^k)=\frac{P(\omega_n^k)+\mathrm{conj}(P(\omega_n^{rev(k)}))}{2}\\ B(\omega_n^k)=i\frac{P(\omega_n^k)-\mathrm{conj}(P(\omega_n^{rev(k)}))}{2}
其中i=1i=\sqrt{-1}conj(x)\mathrm{conj}(x)xx的共轭虚数。

证明:容易发现,我们只需要证明
conj(P(ωnrev(k)))=A(ωnk)iB(ωnk) \mathrm{conj}(P(\omega_n^{rev(k)}))=A(\omega_n^k)-iB(\omega_n^k)
即可得证。

假设
A(x)=i=0n1aixiB(x)=i=0n1bixi A(x)=\sum_{i=0}^{n-1}a_ix^i\\ B(x)=\sum_{i=0}^{n-1}b_ix^i
那么
A(ωnk)iB(ωnk)=j=0n1ajωnkjij=0n1bjωnkj A(\omega_n^k)-iB(\omega_n^k)=\sum_{j=0}^{n-1}a_j\omega_n^{kj}-i\sum_{j=0}^{n-1}b_j\omega_n^{kj}
由于
ωnk=cos2πkn+isin2πkn \omega_n^k=\cos\frac{2\pi k}{n}+i\sin\frac{2\pi k}{n}
带入得
A(ωnk)iB(ωnk)=j=0n1aj(cos2πkjn+isin2πkjn)+j=0n1bi(icos2πkjn+sin2πkjn) A(\omega_n^k)-iB(\omega_n^k)=\sum_{j=0}^{n-1}a_j(\cos \frac{2\pi kj}{n}+i\sin\frac{2\pi kj}{n})+\sum_{j=0}^{n-1}b_i(-i\cos \frac{2\pi kj}{n}+\sin \frac{2\pi kj}{n})
而由于ωnrev(k)=ωnk\omega_n^{rev(k)}=\omega_n^{-k}
conj(P(ωnrev(k)))=conj(A(ωnk)+iB(ωnk)) \mathrm{conj}(P(\omega_n^{rev(k)}))=\mathrm{conj}(A(\omega_n^{-k})+iB(\omega_n^{-k}))

conj(P(ωnrev(k)))=conj(j=0n1aj(cos(2πkjn)+isin(2πkjn))+j=0n1bj(icos(2πkjn)sin(2πkjn))) \mathrm{conj}(P(\omega_n^{rev(k)}))=\mathrm{conj}(\sum_{j=0}^{n-1}a_j(\cos(-\frac{2\pi kj}{n})+i\sin(-\frac{2\pi kj}{n}))+\sum_{j=0}^{n-1}b_j(i\cos(-\frac{2\pi kj}{n})-\sin(-\frac{2\pi kj}{n})))
由于sin(x)=sinx,cos(x)=cosx\sin (-x)=-\sin x,\cos(-x)=\cos x
conj(P(ωnrev(k)))=conj(j=0n1aj(cos2πkjnisin2πkjn)+j=0n1bj(icos2πkjn+sin2πkjn)) \mathrm{conj}(P(\omega_n^{rev(k)}))=\mathrm{conj}(\sum_{j=0}^{n-1}a_j(\cos \frac{2\pi kj}{n}-i\sin \frac{2\pi kj}{n})+\sum_{j=0}^{n-1}b_j(i\cos\frac{2\pi kj}{n}+\sin \frac{2\pi kj}{n}))
由于aj,bja_j,b_j是实数,因此
conj(P(ωnrev(k)))=A(ωnk)iB(ωnk) \mathrm{conj}(P(\omega_n^{rev(k)}))=A(\omega_n^k)-iB(\omega_n^k)
因此,我们将两个多项式的FFT变成了一个多项式的FFT。

MTT

对于两个多项式
A(x),B(x) A(x),B(x)
求多项式乘法, mod 109+7\bmod 10^9+7之类的质数。

不能NTT,FFT的精度误差太大,怎么办?

假设这个质数为pp,令p=pp'=\lfloor\sqrt{p}\rfloor
A(x)=pQA(x)+RA(x)B(x)=pQB(x)+RB(x) A(x)=p'Q_A(x)+R_A(x)\\ B(x)=p'Q_B(x)+R_B(x)
其中RA(x)R_A(x)RB(x)R_B(x)的每一项的系数都<p< p',那么
A(x)B(x)=p2QA(x)QB(x)+p(QA(x)RB(x)+QB(x)RA(x))+RA(x)RB(x) A(x)B(x)=p'^2Q_A(x)Q_B(x)+p'(Q_A(x)R_B(x)+Q_B(x)R_A(x))+R_A(x)R_B(x)
这样能大大减小精度误差。

常数优化

使用上面FFT部分中提到的常数优化2即可。

多项式求逆

A(x)F(x)=1mod  xn A(x)F(x)=1\mod{x^n}

其中nn是多项式A(x)A(x)的度数。

假设我们已经求出了
A(x)F0(x)=1mod  xn/2 A(x)F_0(x)=1\mod x^{\lceil n/2\rceil}
上面两式相减
F(x)F0(x)=0mod  xn/2 F(x)-F_0(x)=0\mod x^{\lceil n/2\rceil}
两边平方
F2(x)2F(x)F0(x)+F02(x)=0mod  xn F^2(x)-2F(x)F_0(x)+F_0^2(x)=0\mod x^n
左右同时乘以A(x)A(x)
F(x)2F0(x)+A(x)F02(x)=0mod  xn F(x)-2F_0(x)+A(x)F_0^2(x)=0\mod x^n

F(x)=2F0(x)A(x)F02(x)mod  xn F(x)=2F_0(x)-A(x)F_0^2(x)\mod x^n

常数优化

对于MTT在多项式求逆上的常数优化,参见这篇文章:一个多项式求逆的卡常技巧

多项式取对数

F(x)=lnA(x)mod  xn F(x)=\ln A(x)\mod x^{n}

考虑对两边同时求导
F(x)=A1(x)A(x)mod  xn F'(x)=A^{-1}(x)A'(x)\mod x^n
因此
F(x)=A1(x)A(x)dx F(x)=\int A^{-1}(x)A'(x)\mathrm dx

牛顿迭代

对于一个多项式F(x)F(x),满足
G(F(x))=0mod  xn G(F(x))=0\mod x^n
G(x)G(x)已知,现在要求F(x)F(x)

假设已经求出了
G(F0(x))=0mod  xn/2 G(F_0(x))=0\mod x^{\lceil n/2\rceil}
考虑对G(F(x))G(F(x))F0(x)F_0(x)处展开
G(F(x))=G(F0(x))+G(F0(x))1!(F(x)F0(x))+G(F0(x))2!(F(x)F0(x))2+ G(F(x))=G(F_0(x))+\frac{G'(F_0(x))}{1!}(F(x)-F_0(x))+\frac{G''(F_0(x))}{2!}(F(x)-F_0(x))^2+\cdots
注意这里的G(F(x))G'(F(x))G(F(x))G(F(x))F(x)F(x)取导数的结果。

由于F(x)F(x)F0(x)F_0(x)的最后n2\lceil \frac{n}{2}\rceil项相同,因此
(F(x)F0(x))k=0mod  xn (F(x)-F_0(x))^k=0\mod x^n
k2k\geq 2时都成立。

因此
G(F(x))=G(F0(x))+G(F0(x))(F(x)F0(x))mod  xn G(F(x))=G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))\mod x^n
由于
G(F0(x))=0mod  xn G(F_0(x))=0\mod x^n
因此
F(x)=F0(x)G(F0(x))G(F0(x))mod  xn F(x)=F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}\mod x^n

多项式求指数

F(x)=eA(x)mod  xn F(x)=e^{A(x)}\mod x^n

两边同时取对数得到
lnF(x)=A(x)mod  xn \ln F(x)=A(x)\mod x^n

G(F(x))=lnF(x)A(x)=0mod  xn G(F(x))=\ln F(x)-A(x)=0\mod x^n
容易发现
G(F(x))=F1(x)mod  xn G'(F(x))=F^{-1}(x)\mod x^n
套牛顿迭代公式,假设我们已经求出了
lnF0(x)A(x)=0mod  xn/2 \ln F_0(x)-A(x)=0\mod x^{\lceil n/2\rceil}

F(x)=F0(x)F0(x)(lnF0(x)A(x))mod  xn F(x)=F_0(x)-F_0(x)(\ln F_0(x)-A(x))\mod x^n
更优雅的式子是
F(x)=F0(x)(1lnF0(x)+A(x))mod  xn F(x)=F_0(x)(1-\ln F_0(x)+A(x))\mod x^n

多项式开方

F2(x)=A(x)mod  xn F^2(x)=A(x)\mod x^n

考虑
G(F(x))=F2(x)A(x)=0mod  xn G(F(x))=F^2(x)-A(x)=0\mod x^n
容易发现
G(F(x))=2F(x)mod  xn G'(F(x))=2F(x)\mod x^n
套牛顿迭代公式,假设我们已经求出了
F02(x)A(x)=0mod  xn/2 F_0^2(x)-A(x)=0\mod x^{\lceil n/2\rceil}
即有
F(x)=F0(x)21F01(x)(F02(x)A(x))mod  xn F(x)=F_0(x)-2^{-1}F_0^{-1}(x)(F_0^2(x)-A(x))\mod x^n
更优雅的式子是
F(x)=21(F0(x)+A(x)F01(x))mod  xn F(x)=2^{-1}(F_0(x)+A(x)F_0^{-1}(x))\mod x^n

拉格朗日反演


G(F(x))=x G(F(x))=x
已知G(x)G(x)要求F(x)F(x)。其中,G(x)G(x)无常数项。

经过一些奥妙重重的变换可以得到
[xn]F(x)=1n[xn1]1Gn(x) [x^n]F(x)=\frac{1}{n}[x^{n-1}]\frac{1}{G'^n(x)}
其中
G(x)=G(x)x G'(x)=\frac{G(x)}{x}
为啥?我也不知道……

UPD:证明