快速傅里叶变换 FFT

前置内容

卷积

两个长度分别为 $n+1$ 和 $m+1$ 的序列 $a,b$ 的卷积定义为一个长度为 $n+m+1$ 的序列 $c$ ，其中 $c$ 的每一项满足：

$$c_i=\sum_{j=0}^i a_jb_{i-j}$$

多项式

两个次数分别为 $n$ 和 $m$ 的多项式 $A(x)=\displaystyle\sum_{i=0}^n a_ix^i$ 和 $B(x)=\displaystyle\sum_{i=0}^m b_ix^i$ 的乘法结果定义为一个次数为 $n+m$ 的多项式 $C(x)=\displaystyle\sum_{i=0}^{n+m} c_ix^i$ ，其中 $c$ 的每一项满足：

$$c_i=\sum_{j=0}^i a_jb_{i-j}$$

不难发现，多项式乘法也可以看作是两个多项式系数序列的卷积。

引理 1 ：

称一个次数为 $n-1$ 的多项式为 $n$ 项多项式。那么一个 $n$ 项多项式 $F(x)=\displaystyle\sum_{i=0}^{n-1}f_ix^i$ 可以由恰好 $n$ 个点对 $(x_i,y_i)$ 确定。其中对于每个 $y_i$ 满足 $y_i=\displaystyle\sum_{j=0}^{n-1}f_jx_i$ 。

引理 2 ：

对于 $n$ 个点对 $(x_i,y_i)$ ，这些点对可以恰好确定一个项数不超过 $n$ 的多项式 $F(x)=\displaystyle\sum_{i=0}^{n-1}f_ix^i$ ，其中对于每个 $y_i$ 满足 $y_i=\displaystyle\sum_{j=0}^{n-1}f_jx_i$ 。

由上面的两个引理，可以轻松得到多项式有两个表达：一个是系数表达，一个是点值表达。

DFT 和 FFT

DFT（ Discrete Fourier Transform ，离散傅里叶变换）是一种将信号从时域转移到频域的方法。放在多项式中就是把多项式的系数表达转化成多项式的点值表达。

IDFT （Inverse Discrete Fourier Transform ，离散傅里叶逆变换）是一种将信号从频域转移到时域的方法。放在多项式中就是把多项式的点值表达转化成多项式的系数表达。

FFT （ Fast Fourier Transform ，快速傅里叶变换）和 IFFT （ Inverse Fast Fourier Transeform ，快速傅里叶逆变换）分别为 DFT 和 IDFT 的快速算法。

可以发现， DFT 和 IDFT 结合起来，可以轻松做到多项式乘法。原因是多项式乘法中，设有 $n+m+1$ 个点对 $(x_i,a_i)$ 表示 $n+1$ 项多项式 $A(x)$ ，有 $n+m+1$ 个点对 $(x_i,b_i)$ 表示 $m+1$ 项多项式 $B(x)$ ，有 $n+m+1$ 项多项式 $C(x)=A(x)B(x)$ ，有 $n+m+1$ 个点对 $(x_i,c_i)$ 表示 $C(x)$ 。那么对于任意的 $i$ ，都有 $c_i=a_ib_i$ ，这是显然的。那么只需要对 $A(x),B(x)$ 做 DFT ，对 $c$ 做 IDFT ，就可以完成多项式乘法。

而对于 DFT 和 IDFT 的计算，我们是有 $O(n^2)$ 的算法的。

DFT 的 $O(n^2)$ 算法：

显然对于 $n$ 项多项式 $A(x)$ ，随便取 $n$ 个 $x_i$ ，作为点值表达的横坐标。接下来就是多项式多点求值。这个是显然可以做到 $O(n^2)$ 的。于是 DFT 可以在 $O(n^2)$ 的时间内完成。

IDFT 的 $O(n^2)$ 算法：

引理 3 （拉格朗日插值）：

对于 $n$ 个点对 $(x_i,y_i)$ ，这些点可以唯一确定一个不超过 $n$ 项的多项式 $F(x)$ ：

$$F(x)=\sum_{i=0}^{n-1} y_i\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}$$

引理 3 证明：

显然多项式项数不超过 $n$ 且多项式唯一。

不难发现，对于任意的 $i$ ， $F(x_i)=y_i\frac{(x_i-x_0)(x_i-x_1)\cdots}{(x_i-x_0)(x_i-x_1)\cdots}+\displaystyle\sum_{k=0,k\ne i}^{n-1} y_k\frac{(x_i-x_0)(x_i-x_1)\cdots(x_i-x_i)\cdots}{(x_k-x_0)(x_k-x_1)\cdots(x_k-x_i)\cdots}=y_i+0=y_i$ 。

于是引理 3 得证。

那么 IDFT 可以用拉格朗日插值得到多项式系数，从而 $O(n^2)$ 计算。

复数

记实数集合为 $\mathbb R$ ，整数集合为 \mathbb N$ ，有限域为 $\mathbb F_p$ 。

由中学知识可得，对于任意 $x\in \mathbb R$ ，有 $x^2\ge 0$ 。

于是，数学家们定义， $i^2=-1$ ，对于任意实数 $a,b$ ，记 $a+bi$ 为复数，复数集合为 $\mathbb C$ 。

定义复数 $z=a+bi$ 的共轭为 $conj(z)=a-bi$ 。

复数四则运算：

1
2
3
4

(a+bi)+(c+di)=(a+c)+(b+d)i 
(a+bi)-(c+di)=(a-c)+(b-d)i 
(a+bi)(c+di)=(ac-bd)+(ad+bc)i 
(a+bi)/(c+di)=(a+bi)(c-di)/(c^2+d^2)

单位根

引理 4 （代数基本定理）：

对于任意 $n$ 项整系数多项式 $F(x)$ ， $F(x)=0$ 的复数解恰有 $n-1$ 个（重根也算）。

引理 5 （泰勒展开）：

本引理省去部分内容并省去使用条件及收敛条件。

对于任意函数 $F(x)$ ，记 $F'(x)$ 为 $F(x)$ 的导函数， $F^{(n)}(x)$ 为 $F(x)$ 的 $n$ 阶导函数，则有：

$$F(x)=\sum_{i\ge 0}\frac{F^{(i)}(k)}{i!}(x-k)^i$$

引理 6 （欧拉定理，泰勒展开推论）：

$$e^{ix}=\cos x+i\sin x$$

$$(\cos x+i\sin x)^n=e^{ixn}$$

基于此，定义 $\omega_n$ 为 $\cos\frac{2\pi}{n}+i\sin\frac{2\pi}{n}$ 。称 $\omega_n$ 和 $\omega_n^k$ 为 $n$ 次单位根。

可以得到， $\omega_n^n=(\cos\frac{2\pi}{n}+i\sin\frac{2\pi}{n})^n=e^{in\times\frac{2\pi}{n}}=e^{2i\pi}=1$ ，且 $\omega_n^k$ 在 $0 \le k <n$ 时各不相同。

快速傅里叶变换 FFT

FFT 核心是利用了单位根的特点，选用单位根作为 DFT 的点值横坐标，从而实现快速变换。

设 FFT 序列长度为 $n$ ，且存在一个非负整数 $k$ 使 $2^k=n$ 。

考虑分治计算点值。

首先我们要面临一个问题：怎么分治？

原问题是一个 $n$ 项多项式和 $n$ 个 $n$ 次单位根点值。我们想要分治必然要转化为 $\frac{n}{2}$ 项多项式和 $\frac{n}{2}$ 个点值，特别的，原来的点值是 $n$ 个 $n$ 次单位根，转化后的是 $\frac{n}{2}$ 个 $\frac{n}{2}$ 次单位根。

引理 7 （折半引理）：

$$(\omega_{2n}^k)^2=\omega_{2n}^{2k}=\omega_{n}^k$$

$$\omega_{n}^k=-\omega_n^{k+\frac{n}{2}}$$

证明显然。

我们设 $n$ 项多项式 $A(x)=\displaystyle\sum_{i=0}^{n-1}a_ix^i$ 为 DFT 的式子， $A_1(x)=\displaystyle\sum_{i=0}^{\frac{n}{2}-1}a_{2i}x^i$ ， $A_2(x)=\displaystyle\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}x^i$ 。那么有：

$$A(x)=A_1(x^2)+xA_2(x^2)$$

设 $k<\frac{n}{2}$ ，则有：

$$A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})=A_1(\omega_\frac{n}{2}^k)+\omega_n^kA_2(\omega_\frac{n}{2}^k)$$

$$A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k+n})=A_1(\omega_\frac{n}{2}^{k})-\omega_n^kA_2(\omega_\frac{n}{2}^{k})$$

我们成功的将原序列分治成立两个一样的子问题。

可以发现这样做的复杂度是 $O(n\log n)$ 的。

做完 FFT ，我们还要考虑 IFFT 。

设要变换序列为 $y$ ， IFFT 变换后的结果为 $a$ 。

使用 FTT 对 $y$ 变换出在 $\omega_n^{-k}$ 的点值，设为 $c$ 。这里和上面是有着相同的性质的，可以直接套用。

那么就有：

$$\begin{aligned}c_k & =\sum_{i=0}^{n-1} y_i(\omega_n^k)^i\\ & =\sum_{i=0}^{n-1} \sum_{j=0}^{n-1} a_j(\omega_n^i)^j(\omega_n^k)^i\\ & =\sum_{j=0}^{n-1} \sum_{i=0}^{n-1} a_j(\omega_n^{j-k})^i \end{aligned}$$

当 $j-k=0$ 时，有 $\displaystyle\sum_{i=0}^{n-1} a_j(\omega_n^{j-k})^i=n$ ，

当 $j-k\ne 0$ 时，有 $\displaystyle\sum_{i=0}^{n-1} a_j(\omega_n^{j-k})^i=\frac{1-1}{1-\omega_n^{j-k}}=0$ 。

所以有：

$$c_k=na_k+\displaystyle\sum 0=na_k$$

于是我们就成功的在 $o(n\log n)$ 的时间内做到了 FFT 和 IFFT 。

题目： 多项式乘法 。

主要代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

const double pi=acos(-1);
int n,m,len,r[N],l;
struct C{
    double x,y;
    C(){x=0;y=0;}
    C(double x,double y):x(x),y(y){}
}a[N],b[N],c[N],fa[N],fb[N],fc[N];//复数类
C operator+(C x,C y){return C(x.x+y.x,x.y+y.y);}
C operator-(C x,C y){return C(x.x-y.x,x.y-y.y);}
C operator*(C x,C y){return C(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);}
void FFT(C *a,C *f,int t){
    for(int i=0;i<len;i++)if(i<r[i])swap(a[i],a[r[i]]);
    for(int i=0;i<len;i++)f[i]=a[i];
    for(int i=2;i<=len;i<<=1){
        C G(cos(pi/(i>>1)),t*sin(pi/(i>>1)));
        for(int j=0;j<len;j+=i){
            C g(1,0);
            for(int k=j;k<(j+(i>>1));k++,g=g*G){
                C x=f[k],y=g*f[k+(i>>1)];
                f[k]=x+y;
                f[k+(i>>1)]=x-y;
            }
        }
    }
}
int main(){
    cin>>n>>m;
    for(int i=0;i<=n;i++)cin>>a[i].x;
    for(int i=0;i<=m;i++)cin>>b[i].x;
    l=0;len=1;
    while(len<=n+m)len<<=1,l++;
    for(int i=0;i<len;i++)r[i]=(r[i>>1]>>1)+((i&1)<<(l-1));//奇数偶数分类
    FFT(a,fa,1);FFT(b,fb,1);
    for(int i=0;i<len;i++)fc[i]=fa[i]*fb[i];
    FFT(fc,c,-1);
    for(int i=0;i<=n+m;i++)cout<<(int)(c[i].x/len+0.5)<<' ';//最后要将答案/len
    return 0;
}

容易发现，这里算一次需要做三次 FFT ，再加上大量浮点运算和三角函数，常数与精度都不特别良好。

事实上，无论是算多项式乘法，还是对多个序列进行 FFT 和 IFFT ，只要变换的序列是实数，都可以在两次解决。

多项式乘法主要依托于如下结论： $(A+Bi)^2=(A^2-B^2)+2ABi$ 。

可以看到，只需要对 $(A+Bi)$ 进行 FFT ，自乘，然后 IFFT ，最后提取虚部就可以了。一共需要两次。

而对于两个实数序列的 FFT 和 IFFT ，主要依托于如下结论：

对于两个实数多项式 $A(x),B(x)$ ，

设 $P(x)=A(x)+iB(x),Q(x)=A(x)-iB(x)$ ， $P(x),Q(x)$ 的 FFT 变换结果序列分别为 $p,q$ 。

那么就有：

$$\begin{aligned} p_j & =\sum_{k=0}^{n-1} (a_k+ib_k) \omega_n^{jk} \\ & =\sum_{k=0}^{n-1} (a_k+ib_k) (\cos \frac{2jk\pi}{n}+i\sin \frac{2jk\pi}{n}) \\ q_j & =\sum_{k=0}^{n-1} (a_k-ib_k) (\cos \frac{2jk\pi}{n}+i\sin \frac{2jk\pi}{n}) \\ & =\sum_{k=0}^{n-1} (a_k\cos \frac{2jk\pi}{n}+b_k\sin\frac{2jk\pi}{n})+i(a_k\sin\frac{2jk\pi}{n}-b_k\cos \frac{2jk\pi}{n}) \\ & =conj(p_{n-j}) \end{aligned}$$

于是我们就可以在一次 FFT 中计算出 $A,B$ 的 FFT 变换序列。 IFFT 同理。

可以发现第二种算法比第一种广得多，本文将不在探讨第一种。

第一种 FFT 优化主要代码：

1
2
3
4

for(int i=0;i<=n;i++)cin>>a[i].x;
for(int i=0;i<=m;i++)cin>>a[i].y;
FFT(a,1);for(int i=0;i<len;i++)a[i]=a[i]*a[i];FFT(a,-1);
for(int i=0;i<=n+m;i++)cout<<(int)(a[i].y/(2*len)+0.5)<<' ';

第二种 FFT 优化主要代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

struct C{
    double x,y;
    C(){x=0;y=0;}
    C(double x,double y):x(x),y(y){}
}a[N],b[N];
C operator+(C x,C y){return C(x.x+y.x,x.y+y.y);}
C operator-(C x,C y){return C(x.x-y.x,x.y-y.y);}
C operator*(C x,C y){return C(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);}
C operator/(C x,int y){return C(x.x/y,x.y/y);}
C operator~(C x){return C(x.x,-x.y);}
int main(){
    cin>>n>>m;
    for(int i=0;i<=n;i++)cin>>a[i].x;
    for(int i=0;i<=m;i++)cin>>a[i].y;
    l=0;len=1;
    while(len<=n+m)len<<=1,l++;
    for(int i=0;i<len;i++)r[i]=(r[i>>1]>>1)+((i&1)<<(l-1));
    FFT(a,1);a[len]=a[0];
    for(int i=0;i<len;i++)b[i]=(a[i]+(~a[len-i]))*((~a[len-i])-a[i])/4*(C(0,1));
    FFT(b,-1);
    for(int i=0;i<=n+m;i++)cout<<(int)(b[i].x/len+0.5)<<' ';
    return 0;
}

经过优化的 FFT 的常数有了明显的减小。

前置内容

原根

原根最标准的定义很复杂，这里只介绍性质。

引理 8 ：

设一个数 $P$ 的原根为 $g$ ，则 $g^{\phi(P)}\equiv 1 \pmod P$ ，且 $g^k\equiv K \pmod P$ ，其中 $1\le k <\phi(P)$ 且对于每个 $k$ ， $K$ 互不相同。

引理 9 ：

只有 $1,2,4,p^a$ 有原根。其中 $p$ 为奇素数， $a\ge 1$ 。

引理 10 ：

$g^{(P-1)/n} \equiv \omega_n \pmod P$ ，其中 $g$ 为 $P$ 的原根。

证明较为简单，因为原根 $g$ 与 $\omega_n$ 同样具有循环的性质。

快速数论变换 NTT

由上面的引理 10 ，可以轻松将 FFT 转化为使用原根代替单位根的 NTT （ Number Theoretic Transforms ，快速数论变换）。相比 FFT ，常数小很多很多（连优化后的 FFT 都无法比拟）。唯三的缺点是 NTT 算的是模意义下的，如果数据范围较大可能会有问题；以及 NTT 无法进行像 FFT 一样的优化；还有模数必须为 $a2^b+1$ 的形式且 $b\ge n$ 。

主要代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

const int N=4000010,P=998244353;
int n,m,a[N],b[N],c[N],fa[N],fb[N],fc[N],len,lg[N],r[N],l;
int kpow(int a,int b){
    int t=1;
    while(b){
        if(b&1)t=1ll*t*a%P;
        a=1ll*a*a%P;
        b>>=1;
    }
    return t;
}
void NTT(int *a,int *f,int t){
    for(int i=0;i<len;i++)if(i<r[i])swap(a[i],a[r[i]]);
    for(int i=0;i<len;i++)f[i]=a[i];
    for(int i=2;i<=len;i<<=1){
        int G=kpow(t==1?3:kpow(3,P-2),(P-1)/i);//使用原根代替单位根。 998244353 的原根是 3 。
        for(int j=0;j<len;j+=i){
            int g=1;
            for(int k=j;k<(j+(i>>1));k++){
                int x=f[k],y=1ll*g*f[k+(i>>1)]%P;
                f[k]=(x+y)%P;
                f[k+(i>>1)]=(x-y+P)%P;
                g=1ll*G*g%P;
            }
        }
    }
}
int main(){
    cin>>n>>m;
    for(int i=0;i<=n;i++)cin>>a[i];
    for(int i=0;i<=m;i++)cin>>b[i];
    l=0;len=1;
    while(len<=n+m)len<<=1,l++;
    for(int i=0;i<len;i++)r[i]=(r[i>>1]>>1)+((i&1)<<(l-1));
    NTT(a,fa,1);NTT(b,fb,1);
    for(int i=0;i<len;i++)fc[i]=1ll*fa[i]*fb[i]%P;
    NTT(fc,c,-1);
    int in=kpow(len,P-2);
    for(int i=0;i<=n+m;i++)cout<<1ll*c[i]*in%P<<" ";
    return 0;
}

前置内容

牛顿迭代法

已知一个复合多项式 $G(F(x))\equiv 0 \pmod {x^n}$ ， $F_0(x)\equiv F(x) \pmod {x^\frac{n}{2}}$ 。

引理 11 ：

若 $A(x) \equiv B(x) \pmod {x^\frac{n}{2}}$ ，那么 $(A(x)-B(x))^k \equiv 0 \pmod {x^n}$ 在 $k\ge 2$ 时成立。

证明显然。

那么对 $G(F(x))$ 在 $F_0(x)$ 出进行展开：

$$G(F(x))=\sum_i \frac{G^{(i)}(F_0(x))}{i!}(F(x)-F_0(x))^i$$

由引理 11 ，在 $i \ge 2$ 时后面的项全为 $0$ ，所以有：

$$G(F(x)) \equiv G(F_0(x))+\frac{G'(F_0(x))}{F(x)-F_0(x)} \pmod {x^n}$$

又因为 $G(F(x)) \equiv 0 \pmod {x^n}$ ，那么整理式子可得：

$$F(x) \equiv F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))} \pmod {x^n}$$

这就是牛顿迭代了。可以发现，每进行一次迭代，多项式的项数就会翻倍。

多项式乘法逆

已知 $F(x)$ ，求一个 $G(x)$ 满足 $F(x)G(x) \equiv 1\pmod{x^n}$ 。

整理式子可以得到 $F(x)-\frac{1}{G(x)} \equiv 0\pmod {x^n}$ 。

定义 $H(t)=F(x)-\frac{1}{t}$ ，那么就有 $H(G(x)) \equiv 0 \pmod{x^n}$ 。

直接使用牛顿迭代即可得到：

$$G(x)\equiv 2G_0(x)-F(x)G_0^2(x) \pmod{x^n}$$

复杂度为 $T(n)=T(\frac{n}{2})+O(n\log n)=O(n\log n)$ 。

题目： 多项式乘法逆 。

主要代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

int main(){
    cin>>n;n--;
    for(int i=0;i<=n;i++)cin>>a[i];
    tmpb[0]=b[0]=kpow(a[0],P-2);tmpa[0]=a[0];tmpa[1]=a[1];
    len=2,l=1;
    for(int tt=2;((tt>>1)<=n);tt<<=1){
        len<<=1,l++;
        for(int i=0;i<len;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
        NTT(tmpa,1);NTT(tmpb,1);
        for(int i=0;i<len;i++)tmpa[i]=1ll*tmpa[i]*tmpb[i]%P*tmpb[i]%P;
        NTT(tmpa,-1);
        int in=kpow(len,P-2);
        for(int i=0;i<tt;i++)b[i]=(2ll*b[i]%P-1ll*tmpa[i]*in%P+P)%P;
        for(int i=tt;i<len;i++)b[i]=0;
        for(int i=0;i<len;i++)tmpa[i]=a[i],tmpb[i]=b[i];
    }
    for(int i=0;i<=n;i++)cout<<b[i]<<" ";
    return 0;
}

多项式 ln

已知多项式 $F(x)$ ，求多项式 $G(x)$ 满足 $G(x)\equiv \ln F(x)\pmod{x^n}$ 。

多项式 ln 存在的充要条件为 F(0)=1 。

求解这个得到过程相对简单。

$$G(x) \equiv \ln F(x)\pmod{x^n}$$

$$G'(x)\equiv \frac{F'(x)}{F(x)} \pmod{x^n}$$

$$G(x) \equiv \int\frac{F'(x)}{F(x)} \pmod{x^n}$$

复杂度 $O(n\log n)$ 。

题目： 多项式 ln 。

主要代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

void Inv(int *a){
    len=2,l=1;
    b[0]=tmpb[0]=kpow(a[0],P-2);tmpa[0]=a[0],tmpa[1]=a[1];
    for(int tt=2;(tt>>1)<=n;tt<<=1){
        len<<=1,l++;
        for(int i=0;i<len;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
        NTT(tmpa,1);NTT(tmpb,1);
        for(int i=0;i<len;i++)tmpa[i]=1ll*tmpa[i]*tmpb[i]%P*tmpb[i]%P;
        NTT(tmpa,-1);
        int in=kpow(len,P-2);
        for(int i=0;i<tt;i++)b[i]=(2ll*b[i]%P-1ll*tmpa[i]*in%P+P)%P;
        for(int i=0;i<len;i++)tmpa[i]=a[i],tmpb[i]=b[i];
    }
    for(int i=0;i<=n;i++)a[i]=b[i];
}
void _intx(int *a){
    for(int i=n-1;i>=0;i--)a[i+1]=1ll*kpow(i+1,P-2)*a[i]%P;
    a[0]=0;
}
void _dx(int *a){
    for(int i=1;i<=n;i++)a[i-1]=1ll*a[i]*i%P;
    a[n]=0;
}
void ln(int *a){
    memset(c,0,sizeof(c));
    for(int i=0;i<=n;i++)c[i]=a[i];
    Inv(a);_dx(c);len=1,l=0;
    while(len<=2*n)len<<=1,l++;
    for(int i=0;i<len;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    NTT(a,1);NTT(c,1);for(int i=0;i<len;i++)a[i]=1ll*a[i]*c[i]%P;NTT(a,-1);
    int in=kpow(len,P-2);
    for(int i=0;i<=n;i++)a[i]=1ll*a[i]*in%P;
    for(int i=n+1;i<len;i++)a[i]=0;
    _intx(a);
}

分治 FFT

给定 $G(x)=\sum_i g_ix^i$ ，求一个 $F(x)$ ，满足 $f_0=1,f_i=\displaystyle\sum_{j=1}^i f_{i-j}g_j$ 。

这个是可以用多项式乘法逆解决的。只需要将其看作 $F(x)=F(x)G(x)+1$ 即可。复杂度 $O(n\log n)$ 。

但是有一些时候是不能用分治 FFT 解决的，这时我们就需要用分治 FFT 。

顾名思义，分治 FFT 使用分治计算这个。

假设现在分治要算 $f_{l\cdots r}$ 。

那么我们可以先算 $f_{l\cdots mid}$ ，后算 $f_{mid_1\cdots r}$ 。

重点就是算 $f_{l\cdots mid}$ 对 $f_{mid_1\cdots r}$ 的贡献。

首先，对于右边，我要算的其实是 $f_i=\displaystyle\sum_{j=l}^if_jg_{i-j}$ 。

而对于右边，已经算的有 $f_i=\displaystyle\sum_{j=mid+1}^if_jg_{i-j}$ 。

那么两者的差恰好为 $\displaystyle\sum_{j=l}^mid f_jg_{i-j}$ ，这恰恰是 $f_{l\cdots mid}$ 和 $g_{0\cdots r-l}$ 的卷积，可以 $O(n\log n)$ 计算。

总复杂度 $T(n)=2T(\frac{n}{2})+O(n\log n)=O(n\log^2n)$ 。

题目： 分治 FFT 。

主要代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14

void merge(int l,int r){
    if(l==r)return ;
    int mid=(l+r)>>1;
    merge(l,mid);
    len=r-l+1;
    for(int i=0;i<len;i++)ra[i]=(ra[i>>1]>>1)|((i&1)<<(lg[len]-1));
    memcpy(ta,a,sizeof(int)*len);
    memcpy(tf,f+l,sizeof(int)*(len>>1));
    memset(tf+(len>>1),0,sizeof(int)*(len>>1));
    NTT(ta,1);NTT(tf,1);for(int i=0;i<len;i++)tf[i]=1ll*tf[i]*ta[i]%P;NTT(tf,-1);
    int in=kpow(len,P-2);
    for(int i=mid+1;i<=r;i++)f[i]=(f[i]+1ll*in*tf[i-l]%P)%P;
    merge(mid+1,r);
}

既然多项式乘法逆能解决部分分治 FFT 的问题，那么分治 FFT 也可以解决所有多项式乘法逆问题。

多项式 exp

给定 $F(x)$ ，求 $G(x)$ 满足 $G(x)\equiv e^{F(x)} \pmod{x^n}$ 。

G(x) 存在的充要条件是 F(0)=0 。

题目： 多项式 exp 。

牛顿迭代解法

$$G(x)\equiv e^{F(x)}\pmod{x^n}$$

$$\ln G(x)\equiv F(x)\pmod{x^n}$$

$$\ln G(x)-F(x)\equiv 0\pmod{x^n}$$

令 $H(t)=\ln t-F(x)$ ，那么有 $H(G(x))\equiv 0\pmod{x^n}$ 。

直接牛顿迭代可得：

$$G(x)\equiv G_0(x)(1-\ln G_0(x)+F(x))\pmod{x^n}$$

总复杂度 $O(n\log n)$

主要代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

void exp(int *a,int n){
    d[0]=tmpd[0]=1;
    len=2,l=1;
    for(int tt=2;(tt>>1)<=n;tt<<=1){
        len<<=1,l++;
        int tlen=len,tl=l;
        for(int i=0;i<tt;i++)tmp[i]=tmpd[i];ln(tmp,tt-1);
        e[0]=(1+a[0]-tmp[0]+P)%P;
        for(int i=1;i<tt;i++)e[i]=(a[i]-tmp[i]+P)%P;
        len=tlen,l=tl;
        for(int i=0;i<len;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
        for(int i=tt;i<len;i++)e[i]=tmpd[i]=0;
        NTT(tmpd,1);NTT(e,1);
        for(int i=0;i<len;i++)tmpd[i]=1ll*tmpd[i]*e[i]%P;
        NTT(tmpd,-1);
        int in=kpow(len,P-2);
        for(int i=0;i<tt;i++)d[i]=1ll*tmpd[i]*in%P;
        for(int i=tt;i<len;i++)d[i]=0;
        for(int i=0;i<len;i++)tmpd[i]=d[i];
    }
    for(int i=0;i<=n;i++)a[i]=d[i];
}

分治 FFT 解法

可以预见的是，分治 FFT 做法复杂度为 $O(n\log^2n)$ 。

虽然复杂度不及多项式 exp ，但是分治 FFT 解法常数很小，相比于码量很大，常数也大（本来牛顿迭代的 $O(n\log n)$ 常数就大，多项式 exp 的牛顿迭代更是二次牛顿迭代，多项式 ln 一个，主流程一个）的多项式 exp ，分治 FFT 优势还是很多的。

依旧推公式：

$$G(x)\equiv e^{F(x)}\pmod{x^n}$$

$$G'(x)\equiv e^{F(x)}F'(x)\pmod{x^n}$$

$$G'(x)\equiv G(x)F'(x)\pmod{x^n}$$

$$G(x)\equiv \int G(x)F'(x)\pmod{x^n}$$

这个形式就可以用分治 FFT 求解了。

设 $F'(x)=\displaystyle\sum_i f_ix^i,G(x)=\displaystyle\sum_i g_ix^i$ ，那么就有：

$$g_n=\frac{1}{n}[x^{n-1}]G(x)F'(x)=\frac{1}{n}\sum_{i=0}^{n-1}f_ig_{n-i-1}$$

这就是标准的分治 FFT 形式了。

主要代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14

void cdq(int l, int r)
{
    if (l + 1 == r) return F[l] = l ? (ll)F[l] * ksm(l, mod - 2) % mod : 1, void();
    int mid = (l + r) >> 1; cdq(l, mid); init(r - l);
    for (int i = 0; i < mid - l; ++i) A[i] = F[i + l];
    for (int i = 0; i < r - l - 1; ++i) B[i] = G[i];
    for (int i = mid - l; i < lim; ++i) A[i] = 0;
    for (int i = r - l; i < lim; ++i) B[i] = 0;
    NTT(A, lim, 1); NTT(B, lim, 1);
    for (int i = 0; i < lim; ++i) A[i] = (ll)A[i] * B[i] % mod;
    NTT(A, lim, mod - 2);
    for (int i = mid - l - 1; i < r - l - 1; ++i) (F[i + l + 1] += A[i]) %= mod;
    cdq(mid, r);
}

特殊多项式快速幂 & 多项式开根

给定 $F(x),k$ ，求 $G(x)$ 满足 $G(x)\equiv F(x)^k\pmod{x^n}$ 。保证 $F(0)=1$ 。

这个就相当简单了：

$$F(x)^k=e^{\ln F(x)^k}=e^{k\ln F(x)}$$

只需要多项式 ln 和 exp 即可。

总复杂度 $O(n\log n)$ 。

题目： 多项式快速幂 。

1
2
3

for(int i=0;i<=n;i++)cin>>a[i],aa[i]=a[i];
ln(aa,n);for(int i=0;i<=n;i++)aa[i]=1ll*aa[i]*k%P;exp(aa,n);
for(int i=0;i<=n;i++)cout<<aa[i]<<' ';

多项式开根也一样。注意要保证 F(0)=1 。

题目： 多项式开根 。

主要代码：

1
2
3

for(int i=0;i<=n;i++)cin>>a[i],aa[i]=a[i];
ln(aa,n);for(int i=0;i<=n;i++)aa[i]=1ll*aa[i]*k%P;exp(aa,n);
for(int i=0;i<=n;i++)cout<<aa[i]<<' ';

一般多项式快速幂 & 多项式开根

其实只是不保证 $F(0)=1$ 而已。

当 $F(0)=0$ 时，可以令 $F_1(x)=\frac{F(x)}{x}$ ，然后对 $F_1(x)$ 做快速幂。

当 $F(0)>1$ 时，可以令 $F_1(x)=\frac{F(x)}{F(0)}$ ，然后对 $F_1(x)$ 做快速幂。

题目： 一般多项式快速幂 。

主要代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

int main(){
    scanf("%d%s",&n,str);
    int ln = strlen(str);
    for(reg int i=0;i!=ln;++i){
        k1 = (10ll*k1+str[i]-'0')%p;
        k2 = (10ll*k2+str[i]-'0')%(p-1);
    }
    for(reg int i=0;i!=min(6,ln);++i) k3 = 10*k3+str[i]-'0';
    for(reg int i=0;i!=n;++i) scanf("%d",&F[i]);
    if(F[0]==0&&k3>=n){
        for(reg int i=0;i!=n;++i) printf("0 ");
        return 0;
    }
    init(n<<1|1);
    power(F,n-1,k1,k2);
    return 0;
}

多项式开根也一样。不过必须保证 $F(0)$ 在模意义下时二次剩余，不然不存在模意义下答案。

题目： 一般多项式开根 。

主要代码：

1
2
3
4
5
6
7
8
9

int main()
{
    read(n);
    limit=1;
    while (limit<n) limit<<=1;
    for (i=0;i<n;i++) read(f[i]);gg(limit<<1);
    polyksm1(f,g,limit,ksm(2,p-2));
    for (i=0;i<n;i++) printf("%d ",g[i]);
}

第二类斯特林数·列（分治 FFT ）

题目： 第二类斯特林数·列 。

第二类斯特林数递推公式：

$$S(n,k)=S(n-1,k-1)+kS(n-1,k)$$

打包成生成函数。令 $F_k=\displaystyle\sum_i S_{i,k}x^i$ ，则有：

$$F_k=xF_{k-1}+kxF_k$$

$$F_k=\frac{x}{1-kx}F_{k-1}$$

$$F_k=x^k\prod_{i=1}^k \frac{1}{1-ix}$$

这是 $k$ 个一次函数的乘积。考虑分治 FFT ，每次递归左右两边，然后合并答案。总复杂度 $O(n\log^2n)$ 。

主要代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85

void NTT(int *f,int t){
    for(int i=0;i<len;i++)if(i<R[i])swap(f[i],f[R[i]]);
    for(int i=2;i<=len;i<<=1){
        int G=kpow(t==1?3:kpow(3,P-2),(P-1)/i);
        for(int j=0;j<len;j+=i){
            int g=1;
            for(int k=j;k<(j+(i>>1));k++){
                int x=f[k],y=1ll*g*f[k+(i>>1)]%P;
                f[k]=(x+y)%P;
                f[k+(i>>1)]=(x-y+P)%P;
                g=1ll*G*g%P;
            }
        }
    }
    if(t==-1){
        int in=kpow(len,P-2);
        for(int i=0;i<len;i++)f[i]=1ll*f[i]*in%P;
    }
}
void NTT(vector<int>&f,int t){
    int len=f.size();
    for(int i=0;i<len;i++)if(i<R[i])swap(f[i],f[R[i]]);
    for(int i=2;i<=len;i<<=1){
        int G=kpow(t==1?3:kpow(3,P-2),(P-1)/i);
        for(int j=0;j<len;j+=i){
            int g=1;
            for(int k=j;k<(j+(i>>1));k++){
                int x=f[k],y=1ll*g*f[k+(i>>1)]%P;
                f[k]=(x+y)%P;
                f[k+(i>>1)]=(x-y+P)%P;
                g=1ll*G*g%P;
            }
        }
    }
    if(t==-1){
        int in=kpow(len,P-2);
        for(int i=0;i<f.size();i++)f[i]=1ll*f[i]*in%P;
    }
}
vector<int> merge(int l,int r){
    if(l==r){
        vector<int>tt(2);
        if(l<=m)
        tt[1]=P-l,tt[0]=1;
        else tt[1]=0,tt[0]=1;
        return tt;
    }
    int mid=(l+r)>>1;
    vector<int>a=merge(l,mid);
    vector<int>b=merge(mid+1,r);
    int L=lg[r-l+1]+1;
    for(int i=0;i<(1<<L);i++)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    a.resize(1<<L);b.resize(1<<L);
    NTT(a,1);NTT(b,1);
    for(int i=0;i<a.size();i++)a[i]=1ll*a[i]*b[i]%P;
    NTT(a,-1);
    return a;
}
signed main(){
    cin>>n>>m;
    for(int i=2;i<=n;i++)lg[i]=lg[i>>1]+1;
    vector<int>a=merge(1,m);
    a.resize(800010);
    tmpb[0]=b[0]=kpow(a[0],P-2);tmpa[0]=a[0];tmpa[1]=a[1];
    len=2,l=1;
    n-=m;
    for(int tt=2;((tt>>1)<=n);tt<<=1){
        len<<=1,l++;
        for(int i=0;i<len;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(l-1));
        NTT(tmpa,1);NTT(tmpb,1);
        for(int i=0;i<len;i++)tmpa[i]=1ll*tmpa[i]*tmpb[i]%P*tmpb[i]%P;
        NTT(tmpa,-1);
        int in=kpow(len,P-2);
        for(int i=0;i<tt;i++)b[i]=(2ll*b[i]%P-1ll*tmpa[i]%P+P)%P;
        for(int i=tt;i<len;i++)b[i]=0;
        for(int i=0;i<len;i++)tmpa[i]=a[i],tmpb[i]=b[i];
    }
    if(n<0){
        for(int i=0;i<=n+m;i++)cout<<"0 ";
        return 0;
    }
    for(int i=0;i<m;i++)cout<<"0 ";
    for(int i=0;i<=n;i++)cout<<b[i]<<" ";
    return 0;
}

前置内容

减法卷积

如同卷积， $a$ 和 $b$ 的减法卷积 $c$ 定义为 $c_i=\displaystyle\sum_{j-k=i}a_jb_k$，或者说：$c_i=\displaystyle\sum_j a_jb_{j+i}$ 。

解决这个的方法也很简单。将 b 序列翻转后 $c_i=\displaystyle\sum_{j=0}^{n-i} a_jb_{n-j-i}$ ，构成卷积，可以求解。

函数平移

设一个 $n$ 项多项式函数 $F(x)$ ，和 $G(x)=F(x+k)$ 。

设 $F(x)=\displaystyle\sum_i f_ix^i,G(x)=\displaystyle\sum_i g_ix^i$ 。

则有：

$$\begin{aligned} G(x) &=\sum_{i=0}^{n-1} f_i(x+k)^i \\ &=\sum_{i=0}^{n-1}\sum_{j=0}^i C_i^jf_ix^jk^j \\ &=\sum_{j=0}^{n-1}\sum_{i=j}^{n-1} \frac{i!}{j!(i-j)!}f_ix^jk^j \\ &=\sum_{j=0}^{n-1} x^j\frac{k^j}{j!}\sum_{i=j}^{n-1} f_i\frac{i!}{(i-j)!} \end{aligned}$$

容易发现后面构成减法卷积。令 $h_i=f_{n-1-i}(n-1-i)!$ ，则有：

$$\begin{aligned} g_j&=\frac{k^j}{j!}\sum_{i=j}^{n-1} f_i\frac{i!}{(i-j)!} \\ &=\frac{k^j}{j!}\sum_{i=0}^{n-1-j} f_{i+j}\frac{(i+j)!}{i!} \\ &=\frac{k^j}{j!}\sum_{i=0}^{n-1-j} h_{n-1-i-j}\frac{1}{i!} \end{aligned}$$

这是一个卷积形式，可以直接求解。

所以函数平移可以在 $O(n\log n)$ 时间内求解。

主要代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

void get(vector<int>&x){
    vector<int>y(x.begin(),x.end()),z(x.begin(),x.end());
    for(int i=0;i<y.size();i++){
        y[i]=1ll*y[i]*fac[i]%P*kpow(n>>1,i)%P;
    }
    reverse(y.begin(),y.end());
    for(int i=0;i<x.size();i++){
        z[i]=inv[i];
    }
    mul(z,y);
    z.resize((n>>1)+1);
    reverse(z.begin(),z.end());
    for(int i=0;i<z.size();i++)z[i]=1ll*z[i]*kpow(n>>1,P-1-i)%P*inv[i]%P;
    mul(x,z);
    x.resize(n+1);
}