FFT理性瞎扯
xtx1092515503
2021-06-29 08:55:35
> 大家好，这里是FFT理性瞎扯。
>
> 什么？你问[这](https://www.luogu.com.cn/blog/Troverld/fft-gan-xing-xia-che)是什么？
>
> 啊，那是FFT感性瞎扯。
>
> 也就是说，有一份感性瞎扯和一份理性瞎扯？
>
> 没错！上方链接中给出的，是很久以前我写的，近似于黑历史的粗浅的FFT感性瞎扯。而这份，是不那么像黑历史的不那么粗浅的FFT理性瞎扯。
>
> 实际上，前一份感性瞎扯注重的是FFT的器用，类似于黑匣子原理解说，而这份理性瞎扯更注重于FFT的本质。
>
> 也因此，更多奇奇怪怪的卷积科技，也会一起被塞在这篇被不知道哪个蠢货一拍大腿就题名为“FFT理性瞎扯”的学习笔记里吧。
>
> 什么？你问具体有哪些科技？
>
> ……目前是打算把循环卷积、单位根反演、Bluestein等东西全部加进去吧。
>
> 加完了吗？
>
> （……咕咕咕）
# I.DFT&IDFT的本质
DFT的本质是什么？
是变换：$f(x)=\sum\limits_{i=0}^{n-1}f_ix^i\rightarrow\hat f(x)=\sum\limits_{i=0}^{n-1}x^i\sum\limits_{j=0}^{n-1}f_j\omega_n^{ij}$。
也是变换：$\begin{bmatrix}f_0\\f_1\\\vdots\\f_{n-1}\end{bmatrix}\rightarrow\begin{bmatrix}f_0\\f_1\\\vdots\\f_{n-1}\end{bmatrix}\begin{bmatrix}1&1&\dots&1\\1&\omega_n^1&\dots&\omega_n^{n-1}\\\vdots&\vdots&\ddots&\vdots\\1&\omega_n^{n-1}&\dots&\omega_n^{(n-1)^2}\end{bmatrix}$
其中，$\omega_n$ 是 $n$ 次单位根，实质是虚数 $\cos\dfrac{2\pi}{n}+i\sin\dfrac{2\pi}{n}$。
这其中，我们令转移矩阵 $G=\begin{bmatrix}1&1&\dots&1\\1&\omega_n^1&\dots&\omega_n^{n-1}\\\vdots&\vdots&\ddots&\vdots\\1&\omega_n^{n-1}&\dots&\omega_n^{(n-1)^2}\end{bmatrix}$。同时，我们将多项式 $f$ 的向量形式记作 $\mathbf f$。则有变换本质是 $\mathbf f\rightarrow\mathbf fG=\hat{\mathbf f}$。
假如这里是感性瞎扯，我就会不厌其烦地分析它的复杂度、实现方法之类的吧。
可惜，这里是理性瞎扯，因此就没有愉快的OI环节了。
IDFT的本质是什么？
是变换 $\hat f(x)=\sum\limits_{i=0}^{n-1}x^i\sum\limits_{j=0}^{n-1}f_j\omega_n^{ij}\rightarrow f(x)=\sum\limits_{i=0}^{n-1}f_ix^i$。
也是 $\hat{\mathbf f}\rightarrow\hat{\mathbf f}G^{-1}=\mathbf f$。
这里我们引入**范德蒙德方阵**的概念：序列 $\{x_0,x_1,\dots,x_{n-1}\}$ 对应的范德蒙德方阵的定义如下：
$$V(\{x\})=\begin{bmatrix}1&1&\dots&1\\x_0&x_1^1&\dots&x_n^{n-1}\\\vdots&\vdots&\ddots&\vdots\\x_0^{n-1}&x_1^{n-1}&\dots&x_{n-1}^{n-1}\end{bmatrix}$$
则，$G=V(\{1,\omega_n,\omega_n^2,\dots,\omega_n^{n-1}\})$。
可以发现，范德蒙德方阵的实质就是多点求值，而其逆就是插值吧。
DFT的实现是什么？
那套分治理论，应该很好理解吧。
IDFT的实现是什么？
恐怕没几个人会能理解，其与DFT的唯一区别就在于单位根要取倒数并在搞完后除以 $n$ 的原因吧。
但是，借用范德蒙德方阵的理论，会发现假如DFT的作用是用 $G$ 乘上 $\mathbf f$ 的话，IDFT的作用就是 $G^{-1}$ 乘上 $\mathbf f$ 罢。
然后，我们IDFT的实操，明显就是使用了 $\dfrac{V(1,\omega_n^{-1},\omega_n^{-2},\dots,\omega_n^{-(n-1)})}n$ 这个矩阵吧。
那么，是否有$G^{-1}=\dfrac{V(1,\omega_n^{-1},\omega_n^{-2},\dots,\omega_n^{-(n-1)})}n$ 呢？
我们就来得出 $G^{-1}$ 到底是什么吧。
考虑拉格朗日插值。
$$f=\sum\limits_{i=0}^{n-1}\hat f_i\prod\limits_{j\neq i}\dfrac{x-\omega_n^j}{\omega_n^i-\omega_n^j}$$
稍微变换一下
$$f=\sum\limits_{i=0}^{n-1}\hat f_i\dfrac{\prod\limits_{j\neq i}x-\omega_n^j}{\prod\limits_{j\neq i}\omega_n^i-\omega_n^j}$$
这时候可以考虑快速插值时的分析方法。考虑分子的东西，首先令
$$G(x)=\prod\limits_{0\leq j<n}x-\omega_n^j$$
这时候，通过暴力展开/理性分析/猜，我们都可以得出 $G(x)=x^n-1$，因为方程 $G(x)=0$ 显然有 $n$ 个根 $\{1,\omega_n,\omega_n^2,\dots,\omega_n^{n-1}\}$，而这些根全都是方程 $x^n-1=0$ 的根，于是就得出 $G(x)\propto x^n-1$，而通过展开最高项得出正比例系数为 $1$，于是就有 $G(x)=x^n-1$。
然后关于分子上的东西就有
$$\begin{aligned}\prod\limits_{j\neq i}x-\omega_n^j&=\dfrac{G(x)}{x-\omega_n^i}\\&=\dfrac{x^n-1}{x-\omega_n^i}\\&=\dfrac{1}{\omega_n^i-x}-\dfrac{x^n}{\omega_n^i-x}\\&=\omega_n^{-i}\times\Big(\dfrac{1}{1-\omega_n^{-i}x}-x^n\times\dfrac{1}{1-\omega_n^{-i}x}\Big)\\&=\omega_n^{-i}\times\Big(\sum\limits_{j=0}\omega_n^{-ij}x^j-x^n\sum\limits_{j=0}\omega_n^{-ij}x^j\Big)&(\text{展开成数列形式})\\&=\omega_n^{-i}\times\Big(\sum\limits_{j=0}\omega_n^{-ij}x^j-\sum\limits_{j=n}\omega_n^{-i(j-n)}x^j\Big)\\&=\omega_n^{-i}\times\Big(\sum\limits_{j=0}\omega_n^{-ij}x^j-\sum\limits_{j=n}\omega_n^{-ij}x^j\Big)&(\omega_n^n=1)\\&=\omega_n^{-i}\times\sum\limits_{j=0}^{n-1}\omega_n^{-ij}x^j\end{aligned}$$
然后关于分母上的东西就有
$$\prod\limits_{j\neq i}\omega_n^i-\omega_n^j=\lim\limits_{x\rightarrow\omega_n^i}\dfrac{G(x)}{x-\omega_n^i}\xlongequal{\text{l'H(此符号表示洛必达法则)}}G'(\omega_n^i)=n\omega_n^{i(n-1)}=n\omega_n^{-i}$$
于是
$$f=\sum\limits_{i=0}^{n-1}\hat f_i\dfrac{\omega_n^{-i}\times\sum\limits_{j=0}^{n-1}\omega_n^{-ij}x^j}{n\omega_n^{-i}}=\sum\limits_{i=0}^{n-1}\hat f_i\dfrac{\sum\limits_{j=0}^{n-1}\omega_n^{-ij}x^j}n$$
我们发现，其就是 $\dfrac{V(1,\omega_n^{-1},\omega_n^{-2},\dots,\omega_n^{-(n-1)})}n$ 这个矩阵。
于是我们便理解了IDFT的本质。
# II.FFT的本质
FFT的本质是什么？
是卷积。
再说一遍？
是……是循环卷积。
没错，假如我们对两个长度为 $n$ 的多项式用FFT相乘，若DFT只展开到长度为 $n$，则IDFT时作用的就也是一个长度为 $n$ 的序列，但实际上的卷积长度可是 $2n-1$。一个 $n$ 个点的点值序列，怎么想都不可能插出一个长度为 $2n-1$ 的多项式，只有可能插出长度为 $n$ 的多项式。
那么，那些更高的位跑哪去了呢？
我们发现，因为DFT $n$ 位时用的是 $\omega_n$ 作为单位根，所以IDFT时，对于高于 $n$ 的位，有 $\omega_n^{i}=\omega_n^{i\bmod n}$。也就是说，这实际上是循环卷积。
于是，若 $h=f\times g$，但是这里的 $\times$ 是FFT展开 $n$ 位的乘法，则有 $h_i=\sum\limits_{j+k\equiv i\pmod n}f_j\times g_k$。
只是在普通卷积的时候，我们FFT的数组开得很长，所以不需要担心是否循环起来了而已。
# III.单位根反演
此章节特别鸣谢：[command-block神仙的博客](https://www.luogu.com.cn/blog/command-block/dan-wei-gen-fan-yan-xiao-ji)。
$$\boxed{[n|a]=\dfrac{\sum_{i=0}^{n-1}\omega_n^{ai}}{n}}$$
> 证：
>
> 考虑等比数列求和。当 $n\nmid a$ 时，$\sum\limits_{i=0}^{n-1}\omega_n^{ai}=\dfrac{1-\omega_n^{an}}{1-\omega_n^{a}}=\dfrac{1-1}{1-\omega_n^a}$，因在任何情形下分母均非零，故此式恒为 $0$。而当 $n|a$ 时，此式显然为 $1$。
推论：
$$\boxed{[a\equiv b\pmod n]=\dfrac{\sum_{i=0}^{n-1}\omega_n^{(a-b)i}}{n}}$$
## III.I.[[LOJ#6485]LJJ 学二项式定理](https://loj.ac/p/6485)
$$[i\equiv j\pmod4]=\dfrac{\sum\limits_{k=0}^3\omega_4^{(i-j)k}}4$$
于是
$$\begin{aligned}\sum\limits_{i=0}^n\dbinom ni\times s^i\times a_{i\bmod 4}&=1/4\times\sum\limits_{j=0}^3\sum\limits_{i=0}^n\dbinom nis^ia_j\sum\limits_{k=0}^3\omega_4^{(i-j)k}\\&=1/4\times\sum\limits_{j=0}^3\sum\limits_{k=0}^3a_j\omega_4^{-jk}\sum\limits_{i=0}^n\dbinom nis^i\omega_4^{ik}\\&=1/4\times\sum\limits_{j=0}^3\sum\limits_{k=0}^3a_j\omega_4^{-jk}(s\omega_4^k+1)^n\end{aligned}$$
$\omega_4$ 在 $\bmod998244353$ 意义下存在。故直接算即可。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
const int phi=mod-1;
int ksm(int x,int y){int z=1;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;return z;}
const int G=3;
const int w=ksm(G,(mod-1)/4);
typedef long long ll;
int T,s,res,a[4];
ll n;
int main(){
scanf("%d",&T);
while(T--){
scanf("%lld%d",&n,&s),res=0;
n%=phi;
for(int i=0;i<4;i++)scanf("%d",&a[i]);
for(int j=0;j<4;j++)for(int k=0;k<4;k++)(res+=1ll*a[j]*ksm(w,phi-j*k)%mod*ksm((1ll*s*ksm(w,k)+1)%mod,n)%mod)%=mod;
res=1ll*res*ksm(4,phi-1)%mod;
printf("%d\n",res);
}
return 0;
}
```
## III.II.[小猪佩奇学数学](https://www.luogu.com.cn/problem/P5591)
$$\begin{aligned}&\sum_{i=0}^n\dbinom nip^i\left\lfloor\dfrac{i}{k}\right\rfloor\\=&\sum_{i=0}^n\dbinom nip^i\times\dfrac{i-i\bmod k}{k}&(\text{整除的实际意义})\\=&\dfrac{\sum\limits_{i=0}^n\dbinom nip^ii-\sum\limits_{i=0}^n\dbinom nip^i(i\bmod k)}{k}\\=&\dfrac1k\Bigg(\sum_{i=1}^n\dbinom{n-1}{i-1}np^i-\sum_{i=0}^n\dbinom nip^i\sum_{j=0}^{k-1}j[i\bmod k=j]\Bigg)&(\text{左式使用了吸收恒等式}\dbinom nmm=\dbinom{n-1}{m-1}n)\\=&\dfrac1k\Bigg(np\sum_{i=0}^{n-1}\dbinom{n-1}{i}p^i-\sum_{i=0}^n\dbinom nip^i\sum_{j=0}^{k-1}j\times\dfrac{\sum\limits_{d=0}^{k-1}\omega_k^{(i-j)d}}k\Bigg)\\=&\dfrac1k\Bigg(np(p+1)^{n-1}-\dfrac1k\sum\limits_{d=0}^{k-1}\sum_{j=0}^{k-1}j\omega_k^{-jd}\sum_{i=0}^n\dbinom nip^i\omega_k^{id}\Bigg)\\=&\dfrac1k\Bigg(np(p+1)^{n-1}-\dfrac1k\sum\limits_{d=0}^{k-1}(1+p\omega_k^d)^n\sum_{j=0}^{k-1}j\omega_k^{-jd}\Bigg)\end{aligned}$$
这时候，考虑令 $A=\sum\limits_{j=0}^{k-1}j\omega_k^{-jd}$。
则 $\omega_k^{-d}A=\sum\limits_{j=0}^{k-1}j\omega_k^{-(j+1)d}$。
二者相减，得到 $(1-\omega_k^{-d})A=\sum\limits_{j=0}^{k-1}j\omega_{k}^{-jd}-\sum\limits_{j=1}^k(j-1)\omega^{-jd}=\sum\limits_{j=1}^{k-1}\omega_k^{-jd}-(k-1)\omega_k^{-kd}$。
显然，有 $\omega_k^{-kd}=1$，则有 $(1-\omega_k^{-d})A=\sum\limits_{j=1}^{k-1}\omega_k^{-jd}-k+1$。
显然，有 $\omega_k^{0}=1$。于是有 $(1-\omega_k^{-d})A=\sum\limits_{j=0}^{k-1}\omega_k^{-jd}-k$。
套上等比数列求和，得到 $(1-\omega_k^{-d})A=\dfrac{1-\omega_k^{-kd}}{1-\omega_k^{-d}}-k$。
显然，又出现了 $\omega_k^{-kd}=1$，则有 $(1-\omega_k^{-d})A=-k$。于是 $A=\dfrac{k}{\omega_k^{-d}-1}$。
带回原式可得
$$\begin{aligned}&\sum_{i=0}^n\dbinom nip^i\left\lfloor\dfrac{i}{k}\right\rfloor\\=&\dfrac1k\Bigg(np(p+1)^{n-1}-\dfrac1k\sum\limits_{d=0}^{k-1}(1+p\omega_k^d)^n\sum_{j=0}^{k-1}j\omega_k^{-jd}\Bigg)\\=&\dfrac1k\Bigg(np(p+1)^{n-1}-\dfrac1k\sum\limits_{d=0}^{k-1}(1+p\omega_k^d)^n\times\dfrac{k}{\omega_k^{-d}-1}\Bigg)\\=&\dfrac1k\Bigg(np(p+1)^{n-1}-\sum\limits_{d=0}^{k-1}\dfrac{(1+p\omega_k^d)^n}{\omega_k^{-d}-1}\Bigg)\end{aligned}$$
这时候，我们看向数据范围，就会发现我们可以枚举这个 $d$ 计算。这样，时间复杂度 $O(k\log n)$。
**但是，需要注意的是！**
当 $\omega_k^{-d}=1$ 时，我们上述关于 $A$ 的分析是不成立的。
这时候，我们重新看向 $A$，发现此时有 $A=\sum\limits_{j=0}^{k-1}j=\dfrac{(k-1)k}2$。
这时就需要特判一下。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
int ksm(int x,int y=mod-2){int z=1;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;return z;}
const int G=3;
int n,p,k,w,res;
int main(){
scanf("%d%d%d",&n,&p,&k),w=ksm(G,(mod-1)/k);
for(int d=0,W=1;d<k;d++,W=1ll*W*w%mod)(res+=1ll*ksm((1ll*W*p+1)%mod,n)*(W!=1?ksm((ksm(W)+mod-1)%mod):1ll*(k-1)*ksm(2)%mod)%mod)%=mod;
printf("%d\n",(1ll*n*p%mod*ksm(p+1,n-1)%mod+mod-res)%mod*ksm(k)%mod);
return 0;
}
```
## III.III.[[UOJ#450][集训队作业2018]复读机](https://uoj.ac/problem/450)
- $d=1$
非常显然，答案即为 $k^n$。
- 其余情况
这时就要用EGF了。
单个颜色的EGF是 $f(x)=\sum\limits_{d|i}\dfrac{x^i}{i!}$。
全部的EGF是 $f^k$。
则要求的就是 $n!f^k_n$。
但是 $f$ 的封闭形式难搞，连带着 $f^k$ 的封闭形式也难搞。
咋办呢？
单位根反演。
$$[d|i]=\dfrac{\sum\limits_{j=0}^{d-1}\omega_d^{ij}}{d}$$
带进去就得到
$$f(x)=\dfrac1d\sum\limits_{i=0}\dfrac{\sum\limits_{j=0}^{d-1}\omega_d^{ij}}{i!}x^i=\dfrac1d\sum\limits_{j=0}^{d-1}\sum\limits_{i=0}\dfrac{(x\omega_d^j)^i}{i!}=\dfrac1d\sum\limits_{j=0}^{d-1}e^{x\omega_d^j}$$
- $d=2$
$$f(x)=\dfrac{e^x+e^{-x}}{2}$$
则答案为 $n![x^n]\Big(\dfrac{e^x+e^{-x}}2\Big)^k$。
二项式展开得到 $\dfrac{n!}{2^k}[x^n]\sum\limits_{i=0}^k\dbinom kie^{ix}\times e^{-(k-i)x}=\dfrac{n!}{2^k}[x^n]\sum\limits_{i=0}^k\dbinom kie^{(2i-k)x}$。
提取系数得到 $\dfrac{\sum\limits_{i=0}^k\dbinom ki(2i-k)^n}{2^k}$。
$k\log n$ 直接算即可。
- $d=3$。
$$f(x)=\dfrac{e^x+e^{\omega_3x}+e^{\omega_3^2x}}{3}$$。
答案为 $n![x^n]\Big(\dfrac{e^x+e^{\omega_3x}+e^{\omega_3^2x}}{3}\Big)^k$。
……展开！
$$\dfrac{n![x^n]\sum\limits_{i=0}^k\sum\limits_{j=0}^{k-i}\dbinom ki\dbinom{k-i}je^{ix}\times e^{j\omega_3x}\times e^{(k-i-j)\omega_3^2x}}{3^k}$$
……提取！
$$\dfrac{\sum\limits_{i=0}^k\sum\limits_{j=0}^{k-i}\dbinom ki\dbinom{k-i}j\Big(i+j\omega_3+(k-i-j)\omega_3^2\Big)^n}{3^k}$$
$k^2\log n$ 直接算即可。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
const int mod=19491001;
const int w=663067;
int n,k,d,res,fac[500100],inv[500100];
int ksm(int x,int y=mod-2){int z=1;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;return z;}
int main(){
scanf("%d%d%d",&n,&k,&d);
if(d==1)return printf("%d\n",ksm(k,n)),0;
fac[0]=1;for(int i=1;i<=k;i++)fac[i]=1ll*fac[i-1]*i%mod;
inv[k]=ksm(fac[k]);for(int i=k;i;i--)inv[i-1]=1ll*inv[i]*i%mod;
if(d==2){
for(int i=0;i<=k;i++)(res+=1ll*fac[k]*inv[i]%mod*inv[k-i]%mod*ksm((2*i-k+mod)%mod,n)%mod)%=mod;
res=1ll*res*ksm(ksm(2,k))%mod;
}else{
for(int i=0;i<=k;i++)for(int j=0;i+j<=k;j++)(res+=1ll*fac[k]*inv[i]%mod*inv[j]%mod*inv[k-i-j]%mod*ksm((i+1ll*j*w+1ll*(k-i-j)*ksm(w,2))%mod,n)%mod)%=mod;
res=1ll*res*ksm(ksm(3,k))%mod;
}
printf("%d\n",res);
return 0;
}
```
# IV.循环卷积理论（`Bluestein`算法）
怎么循环卷积？
如果 $n$ 是 $2$ 的整数次幂的话，那么直接FFT就行了吧。
那如果 $n$ 不是 $2$ 的整数次幂呢？
那就正常开两倍数组，然后手动循环吧。
这样如果太慢了呢？
……
$$\begin{aligned}h_i&=\sum_{j+k\equiv i\pmod n}f_j\times g_k\\&=\dfrac1n\sum\limits_{j,k}f_jg_k\sum\limits_{t=0}^{n-1}\omega_n^{(j+k-i)t}\\&=\dfrac1n\sum\limits_{t=0}^{n-1}\omega_n^{-it}\sum\limits_{j,k}f_j\omega_n^{jt}g_k\omega_n^{kt}\end{aligned}$$
……然后发现这就从另一种形式上证出了DFT和IDFT的原理。
考虑DFT的式子：
$$\hat f(x)=\sum\limits_{i=0}^{n-1}x^i\sum\limits_{j=0}^{n-1}f_j\omega_n^{ij}$$
我们考虑一些东西来凑出 $ij$。
$$\dfrac{i^2+j^2-(i-j)^2}2=ij$$
代入得到
$$\hat f(x)=\sum\limits_{i=0}^{n-1}x^i\sum\limits_{j=0}^{n-1}f_j\omega_n^{\big(i^2+j^2-(i-j)^2\big)/2}$$
又有 $\omega_n^k=\omega_{2n}^{2k}$。于是有
$$\hat f(x)=\sum\limits_{i=0}^{n-1}\omega_{2n}^{i^2}x^i\sum\limits_{j=0}^{n-1}\omega_{2n}^{j^2}f_j\omega_{2n}^{-(i-j)^2}$$
用 $A(x)=\sum\limits_{i=0}^{n-1}\omega_{2n}^{i^2}f_i$ 和 $B(x)=\sum\limits_{i=0}^{n-1}\omega_{2n}^{-i^2}$ 卷一块即可做到DFT。注意这里的 $B$ 需要倍长，因为 $i-j$ 可以为正或负。
考虑IDFT，式子就变成
$$f(x)=\dfrac 1n\sum\limits_{i=0}^{n-1}x^i\sum\limits_{j=0}^{n-1}\hat f_j\omega_n^{-ij}$$
也可以类似处理。
这就是循环卷积的 `Bluestein` 算法。
由百度翻译，`stein` 有啤酒杯的意思，所以有时本人也会用“蓝酒杯”来称呼这个算法（但是这个称呼不怎么流行）。
## IV.I.[[POJ2821]TN's Kingdom III - Assassination](http://poj.org/problem?id=2821)
题意比较晦涩，反正就是已知多项式 $\alpha,\beta$ 循环卷积卷成了 $\gamma$。现在已知 $\beta,\gamma$，要你倒推出 $\alpha$。
因为 $n$ 次多项式与 $n$ 个点值在循环卷积意义下一一对应，所以循环卷积不需要烦恼甚么求逆之类，只需要对 $\beta,\gamma$ DFT，除一下得出 $\alpha$ 的点值，然后IDFT回来即可。
不知道SPJ出了甚么锅，输出 `double` 时写 `"%.4lf\n"` 死活过不去，改成 `"%.4f\n"` 就过了，不知道咋回事。
代码：
```cpp
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
const double pi=acos(-1);
const int N=(1<<19)+100;
int rev[N],lim=1,LG,n;
struct cp{
double x,y;
cp(){x=y=0;}
cp(double X,double Y){x=X,y=Y;}
friend cp operator+(const cp&u,const cp&v){return cp(u.x+v.x,u.y+v.y);}
friend cp operator-(const cp&u,const cp&v){return cp(u.x-v.x,u.y-v.y);}
friend cp operator*(const cp&u,const double&v){return cp(u.x*v,u.y*v);}
friend cp operator/(const cp&u,const double&v){return cp(u.x/v,u.y/v);}
friend cp operator*(const cp&u,const cp&v){return cp(u.x*v.x-u.y*v.y,u.x*v.y+u.y*v.x);}
friend cp operator/(const cp&u,const cp&v){return cp((u.x*v.x+u.y*v.y)/(v.x*v.x+v.y*v.y),(u.y*v.x-u.x*v.y)/(v.x*v.x+v.y*v.y));}
void print(){printf("(%lf,%lf)",x,y);}
}A[N],B[N],G[N],H[N];
void FFT(cp*a,int tp){
for(int i=0;i<lim;i++)if(i>rev[i])swap(a[i],a[rev[i]]);
for(int md=1;md<lim;md<<=1){
cp rt(cos(pi/md),tp*sin(pi/md));
for(int stp=md<<1,pos=0;pos<lim;pos+=stp){
cp w(1,0);
for(int i=0;i<md;i++,w=w*rt){
cp x=a[pos+i],y=w*a[pos+md+i];
a[pos+i]=x+y;
a[pos+md+i]=x-y;
}
}
}
if(tp==-1)for(int i=0;i<lim;i++)a[i]=a[i]/lim;
}
void Bluestein(cp*a,int tp){
for(int i=0;i<lim;i++)A[i]=B[i]=cp(0,0);
for(int i=0;i<n;i++)A[i]=cp(cos(pi*i*i/n),tp*sin(pi*i*i/n))*a[i];
for(int i=0;i<(n<<1);i++)B[i]=cp(cos(pi*(i-n)*(i-n)/n),-tp*sin(pi*(i-n)*(i-n)/n));
FFT(A,1),FFT(B,1);
for(int i=0;i<lim;i++)A[i]=A[i]*B[i];
FFT(A,-1);
for(int i=0;i<n;i++)a[i]=A[i+n]*cp(cos(pi*i*i/n),tp*sin(pi*i*i/n));
if(tp==-1)for(int i=0;i<n;i++)a[i]=a[i]/n;
}
int main(){
scanf("%d",&n);
while(lim<(n<<2))LG++,lim<<=1;
for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(LG-1));
for(int i=0;i<n;i++)scanf("%lf",&G[i].x);
for(int i=0;i<n;i++)scanf("%lf",&H[i].x);
// for(int i=0;i<n;i++)G[i].print();puts("");
// for(int i=0;i<n;i++)H[i].print();puts("");
Bluestein(G,1),Bluestein(H,1);
for(int i=0;i<n;i++)G[i]=H[i]/G[i];
Bluestein(G,-1);
for(int i=0;i<n;i++)printf("%.4f\n",G[i].x);
return 0;
}
```
## IV.II.[[HNOI2019]白兔之舞](https://www.luogu.com.cn/problem/P5293)
首先先列式子。
我们设乐曲长度是 $i$，则我们从 $L$ 层当中任选 $i$ 层停留即可。这个共有 $\dbinom Li$ 种方案。
然后再设转移矩阵是 $A$，我们便列出了式子为
$$[x,y]\sum\limits_{i=0}^Lx^{i\bmod k}\dbinom LiA^i$$
其中暂记 $[x,y]$ 表示矩阵 $[x,y]$ 位置的元素。老套路，单位根反演一下
$$[x,y]\dfrac1k\sum\limits_{i=0}^L\sum\limits_{t=0}^{k-1}x^t\sum\limits_{j=0}^{k-1}\omega_k^{(i-t)j}\dbinom LiA^i$$
调整顺序
$$[x,y]\dfrac1k\sum\limits_{t=0}^{k-1}x^t\sum\limits_{j=0}^{k-1}\omega_k^{-tj}\sum\limits_{i=0}^L\dbinom Li\omega_k^{ij}A^i$$
二项式展开
$$[x,y]\dfrac1k\sum\limits_{t=0}^{k-1}x^t\sum\limits_{j=0}^{k-1}\omega_k^{-tj}(\omega_k^jA+I)^L$$
这时我们发现全式中仅有 $(\omega_k^jA+I)^L$ 一项是矩阵，且该项结果仅与 $j$ 有关。于是我们设 $F_k=[x,y](\omega_k^jA+I)^L$。则上式变为
$$\dfrac1k\sum\limits_{t=0}^{k-1}x^t\sum\limits_{j=0}^{k-1}\omega_k^{-tj}F_j$$
我们翻出IDFT的式子
$$f(x)=\dfrac1n\sum\limits_{i=0}^{n-1}x^i\sum\limits_{j=0}^{n-1}\hat f_j\omega_n^{-ij}$$
发现两者长得惊人相似——实际上，对 $F_d$ 在长度为 $k$ 的意义下作循环IDFT，就能得到答案。
直接使用蓝酒杯即可。
需要注意的是，本题模数 $p$ 并不保证是NTT模数，故要使用MTT。
然后，$k$ 次单位根，就先求出 $p$ 的原根后，再求其 $\dfrac{p-1}k$ 次幂即可。
但是，需要注意的是！我们上文中提到的 `Bluestein` 的算式 $\dfrac{i^2+j^2-(i-j)^2}2=ij$ 中会在指数上出现 $\dfrac?2$，也即要对某个数开二次剩余，但是其可能不存在。
那咋办？
没关系，`Bluestein` 已经预见到了这种情形，并给了我们另一个算式：$\dbinom{i+j}2-\dbinom i2-\dbinom j2=ij$。
于是
$$\hat F(x)=\dfrac1k\sum\limits_{t=0}^{k-1}x^t\sum\limits_{j=0}^{k-1}\omega_k^{-tj}F_j=\dfrac1k\sum\limits_{t=0}^{k-1}x^t\sum\limits_{j=0}^{k-1}\omega_k^{\tbinom t2+\tbinom j2-\tbinom{t+j}2}F_j$$
然后就可以拆辣。
$$\hat F(x)=\dfrac1k\sum\limits_{t=0}^{k-1}x^t\omega_k^\tbinom t2\sum\limits_{j=0}^{k-1}\omega_k^{\tbinom j2}F_j\omega_k^{-\tbinom{t+j}2}$$
这样之后，设 $A_j=\omega_k^{\tbinom j2}F_j$，$B_j=\omega_k^{-\tbinom{2k-1-j}2}$。于是 $\hat F(x)=\dfrac1k\sum\limits_{t=0}^{k-1}x^t\omega_k^\tbinom t2\sum\limits_{j=0}^{k-1}A_jB_{2k-1-t-j}$。
然后就没有然后了。时间复杂度 $O(n\log n)$。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1<<18;
namespace Poly{
const int G=3;
int lim,rev[N];
struct Polynomial{
int mod,invlim,A[N],B[N];
int ksm(int x,int y){int z=1;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;return z;}
int NTT(int*a,int tp){
for(int i=0;i<lim;i++)if(i>rev[i])swap(a[i],a[rev[i]]);
for(int md=1;md<lim;md<<=1){
int rt=ksm(G,(mod-1)/(md<<1));
if(tp==-1)rt=ksm(rt,mod-2);
for(int stp=md<<1,pos=0;pos<lim;pos+=stp)for(int i=0,w=1;i<md;i++,w=1ll*w*rt%mod){
int x=a[pos+i],y=1ll*w*a[pos+md+i]%mod;
a[pos+i]=(x+y)%mod;
a[pos+md+i]=(x+mod-y)%mod;
}
}
if(tp==-1)for(int i=0;i<lim;i++)a[i]=1ll*a[i]*invlim%mod;
}
void mul(int*a,int*b){
invlim=ksm(lim,mod-2);
for(int i=0;i<lim;i++)A[i]=B[i]=0;
for(int i=0;i<(lim>>1);i++)A[i]=a[i],B[i]=b[i];
NTT(A,1),NTT(B,1);
for(int i=0;i<lim;i++)A[i]=1ll*A[i]*B[i]%mod;
NTT(A,-1);
}
}P[3];
int ksm(int x,int y,int mod){int z=1;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;return z;}
int mod0=P[0].mod=998244353,mod1=P[1].mod=469762049,mod2=P[2].mod=1004535809;
const int inv0_1=ksm(mod0,mod1-2,mod1);
const int inv01_2=ksm(1ll*mod0*mod1%mod2,mod2-2,mod2);
int mod,phi;
void mul(int*a,int*b,int*c,int LG){
lim=1<<LG;for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(LG-1));
// printf("%d:\n",lim);
// for(int i=0;i<lim;i++)printf("%d ",a[i]);puts("");
// for(int i=0;i<lim;i++)printf("%d ",b[i]);puts("");
for(int i=0;i<3;i++)P[i].mul(a,b);
// for(int k=0;k<3;k++){for(int i=0;i<lim;i++)printf("%d ",P[k].A[i]);puts("");}puts("");
for(int i=0;i<lim;i++){
ll tmp=(0ll+P[0].A[i]+1ll*(P[1].A[i]-P[0].A[i]%mod1+mod1)%mod1*inv0_1%mod1*mod0)%(1ll*mod0*mod1);
c[i]=(tmp%mod+1ll*(P[2].A[i]-tmp%mod2+mod2)*inv01_2%mod2*mod0%mod*mod1%mod)%mod;
}
}
}
using namespace Poly;
int K,omega,F[N],A[N],B[N];
namespace Initialization{
int n,L,x,y;
struct Matrix{
int a[3][3];
Matrix(){memset(a,0,sizeof(a));}
int*operator[](const int&x){return a[x];}
friend Matrix operator*(Matrix&u,Matrix&v){
Matrix w;
for(int i=0;i<n;i++)for(int j=0;j<n;j++)for(int k=0;k<n;k++)(w[i][j]+=1ll*u[i][k]*v[k][j]%mod)%=mod;
return w;
}
}W,M;
int KSM(int T){
Matrix R;
for(int i=0;i<n;i++)R[i][i]=1;
for(;T;T>>=1,M=M*M)if(T&1)R=R*M;
return R[x][y];
}
vector<int>v;
int G;
void init(){
scanf("%d%d%d%d%d%d",&n,&K,&L,&x,&y,&mod),x--,y--;
for(int i=0;i<n;i++)for(int j=0;j<n;j++)scanf("%d",&W[i][j]);
phi=mod-1;
int tmp=phi;
for(int i=2;i*i<=tmp;i++){
if(tmp%i)continue;
v.push_back(i);
while(!(tmp%i))tmp/=i;
}
if(tmp!=1)v.push_back(tmp);
for(G=1;;G++){
bool ok=true;
for(auto i:v)if(ksm(G,phi/i,mod)==1){ok=false;break;}
if(ok)break;
}
omega=ksm(G,phi/K,mod);
for(int i=0,w=1;i<K;i++,w=1ll*w*omega%mod){
for(int j=0;j<n;j++)for(int k=0;k<n;k++)M[j][k]=1ll*W[j][k]*w%mod;
for(int j=0;j<n;j++)(++M[j][j])%=mod;
F[i]=KSM(L);
}
}
}
int binom(int x){return (1ll*x*(x-1)/2)%K;}
int main(){
Initialization::init();
int LG=0;while((1<<LG)<(K<<2))LG++;
for(int i=0;i<K;i++)A[i]=1ll*F[i]*ksm(omega,binom(i),mod)%mod;
for(int i=0;i<2*K;i++)B[i]=ksm(omega,K-binom(2*K-1-i),mod);
// for(int i=0;i<K;i++)printf("%d ",A[i]);puts("");
// for(int i=0;i<2*K;i++)printf("%d ",B[i]);puts("");
mul(A,B,F,LG);
int invk=ksm(K,mod-2,mod);
for(int i=0;i<K;i++)printf("%d\n",1ll*invk*F[2*K-1-i]%mod*ksm(omega,binom(i),mod)%mod);
return 0;
}
```