FFT&NTT

FFT快速傅里叶变换<NTT

FFT和NTT是 $O(nlogn)$ 处理两个多项式相乘的算法（FFT<NTT）

前置知识

复数

一个复数可以表示为

$$z=a+iba,b\in R$$

我们把他看做平面上的一个点，横轴代表实数部分，纵轴代表虚数部分

这个点就是 $(a,b)$

我们把它放在极坐标上

[没事，不会极坐标这里有](极坐标系 - 知乎 (zhihu.com))

令$\theta =\arctan \frac{a}{b},x=\sqrt{a^2+b^2}$

那么这个点就是 $(x,\theta )$

则有

$$a=x\cos \theta ,b=x\sin \theta$$

虚数变为了

$$z=x\cos \theta +ix\sin \theta =x(\cos \theta +i\sin \theta )$$

由欧拉公式得(欧拉公式)

$$z=x{e}^{i\theta }$$

这样任意一个虚数可以表示成这样

以上是复数前置知识

正式FFT

单位根

下文中，默认$n$ 为 $2$ 的正整数次幂

在复平面上，以原点为圆心，$1$ 为半径作圆，所得的圆叫单位圆。以圆点为起点，圆的 $n$ 等分点为终点，做 $n$ 个向量，设幅角为正且最小的向量对应的复数为 ${\omega }_n^0$​ ，称为 $n$ 次单位根，n代表长度。

根据复数乘法的运算法则，其余 $n-1$ 个复数为 ${\omega }_n^1,{\omega }_n^2,\dots ,{\omega }_n^{n-1}$ ​

上面是我在别的地方看到的内容，学习FFT的时候我一直不太懂上面的部分，我今天用一个更加简单的方式让大家理解这部分内容

我们本来的多项式是这样的：

$$A(x)=c_0+c_{1}x+c_2{x}^2+\dots +c_n{x}^n$$

因为复数$\omega$有一些特殊性质，我们需要用，并且我们把 $x$ 直接替换 $\omega$ 是没有关系的

所以就变成了：

$$A(\omega )=c_0+c_1\omega +c_2{\omega }^2+\dots +c_n{\omega }^n$$

${\omega }_1,{\omega }_2,{\omega }_3,...,{\omega }_n$不是一样的数，下面的n代表当前数列长度，算法里面不同时候用不同的，这个和他的性质有关系。

说了这么久 $\omega$ 的性质，到底是什么性质呢？

通过欧拉定理

$${\omega }_n^k=\cos k\frac{2\pi }{n}+i\sin k\frac{2\pi }{n}$$

证明一些单位根性质，下面需要用

$1.{\omega }_n^0={\omega }_n^n=1$
这显而易见。
$2.{\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$
证明：

$${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k⟹{\omega }_n^{\frac{n}{2}}=-1$$

$${\omega }_n^{\frac{n}{2}}=\cos \frac{n}{2}\times \frac{2\pi }{n}+i\sin \frac{n}{2}\times \frac{2\pi }{n}=\cos \pi +i\sin \pi =-1$$

$3.{\omega }_{2n}^{2k}={\omega }_n^k$

证明：

$${\omega }_{2n}^{2k}=\cos 2k\frac{2\pi }{2n}+i\sin 2k\frac{2\pi }{2n}=\cos k\frac{2\pi }{n}+i\sin k\frac{2\pi }{n}={\omega }_n^k$$

$4.\sum _{i=0}^{n-1}({\omega }_n^k{)}^i=0$

证明：
第一步是根据等比数列求和公式：

$$\sum _{i=0}^{n-1}({\omega }_n^k{)}^i=\frac{({\omega }_n^k{)}^n-1}{{\omega }_n^k-1}=\frac{({\omega }_n^n{)}^k-1}{{\omega }_n^k-1}=\frac{(1{)}^k-1}{{\omega }_n^k-1}=0$$

傅里叶变换(学名不重要)

他的原理就是把系数表示法变成点值表示法，点值乘起来时间是O(n)，再转为系数表示法

我们取的点是${\omega }_n^0,{\omega }_n^1,...,{\omega }_n^{n-1}$

所对应的函数值我们定义一个函数 $h$ 表示

$h({\omega }_n^k)=c_0+c_1({\omega }_n^k{)}^1+c_1({\omega }_n^k{)}^2+...+c_{n-1}({\omega }_n^k{)}^{n-1}$

所以我们的步骤就变成了这样

那我们应该如何在$O(nlogn)$的复杂度内算出 $h$ 函数呢？

求 $h()$

$$h(x)=c_0+c_{1}x+c_2{x}^2+...+c_{n-1}{x}^{n-1}$$

我们把h函数分成偶数项和奇数项两部分

$$h_0(x)=c_0+c_{2}x+c_4{x}^2+...+c_{n-2}{x}^{\frac{n}{2}-1}$$

$$h_1(x)=c_1+c_{3}x+c_5{x}^2+...+c_{n-1}{x}^{\frac{n}{2}-1}$$

$$h(x)=h_0(x^2)+x{h}_1(x^2)$$

$$h({\omega }_n^k)=h_0({\omega }_n^{2k})+{\omega }_n^k{h}_1({\omega }_n^{2k})$$

通过可以推出${\omega }_{2n}^{2k}={\omega }_n^k$

$h({\omega }_n^k)=h_0({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k{h}_1({\omega }_{\frac{n}{2}}^k)$

同理，将 ${\omega }_n^{k+\frac{n}{2}}$ 代入得

$h({\omega }_n^{k+\frac{n}{2}})=h_0({\omega }_n^{2k+n})+{\omega }_n^{k+\frac{n}{2}}{h}_1({\omega }_n^{2k+n})$

因为${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$

$h({\omega }_n^{k+\frac{n}{2}})=h_0({\omega }_n^{2k}{\omega }_n^n)-{\omega }_n^k{h}_1({\omega }_n^{2k}{\omega }_n^n)$

因为${\omega }_n^n=1$

$h({\omega }_n^{k+\frac{n}{2}})=h_0({\omega }_n^{2k})-{\omega }_n^k{h}_1({\omega }_n^{2k})$

$h({\omega }_n^{k+\frac{n}{2}})=h_0({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k{h}_1({\omega }_{\frac{n}{2}}^k)$

发现 $h({\omega }_n^k)$ 和 $h({\omega }_n^{k+\frac{n}{2}})$ 刚好是一加一减

我们在枚举第一个式子的时候也可以求出第二个式子的值啦

$n$ 代表当前数列长度，每次减半，所以是$log(n)$

你是不是还是脑子里依托浆糊，我们来搞一个例子推一下(一个大括号里的前一个式子是我们需要的式子，后一个是顺便求出的)

假设我们求$h({\omega }_8^1)$这个数列一共八位

$$\{\begin{array}{l}h({\omega }_8^1)=h_0({\omega }_4^1)+{\omega }_8^1{h}_1({\omega }_4^1)\\ h({\omega }_8^5)=h_0({\omega }_4^1)-{\omega }_8^5{h}_1({\omega }_4^1)\end{array}$$

$h_{00}$ 代表偶数中的偶数，也就是这个数二进制下的末尾两位是不是00

$$\{\begin{array}{l}{h}_0({\omega }_4^1)=h_{00}({\omega }_2^1)+{\omega }_4^1{h}_{01}({\omega }_2^1)\\ h_0({\omega }_4^3)=h_{00}({\omega }_2^1)-{\omega }_4^3{h}_{01}({\omega }_2^1)\end{array}$$

$$\{\begin{array}{l}{h}_1({\omega }_4^1)=h_{10}({\omega }_2^1)+{\omega }_4^1{h}_{11}({\omega }_2^1)\\ h_1({\omega }_4^3)=h_{10}({\omega }_2^1)-{\omega }_4^3{h}_{11}({\omega }_2^1)\end{array}$$

继续推下去

$$\{\begin{array}{l}{h}_{00}({\omega }_2^1)=h_{000}({\omega }_1^1)+{\omega }_2^1{h}_{001}({\omega }_1^1)\\ h_{00}({\omega }_2^2)=h_{000}({\omega }_1^1)-{\omega }_2^2{h}_{001}({\omega }_1^1)\end{array}$$

$$\{\begin{array}{l}{h}_{01}({\omega }_2^1)=h_{010}({\omega }_1^1)+{\omega }_2^1{h}_{011}({\omega }_1^1)\\ h_{01}({\omega }_2^2)=h_{010}({\omega }_1^1)-{\omega }_2^2{h}_{011}({\omega }_1^1)\end{array}$$

$$\{\begin{array}{l}{h}_{10}({\omega }_2^1)=h_{100}({\omega }_1^1)+{\omega }_2^1{h}_{101}({\omega }_1^1)\\ h_{10}({\omega }_2^2)=h_{100}({\omega }_1^1)-{\omega }_2^2{h}_{101}({\omega }_1^1)\end{array}$$

$$\{\begin{array}{l}{h}_{11}({\omega }_2^1)=h_{110}({\omega }_1^1)+{\omega }_2^1{h}_{111}({\omega }_1^1)\\ h_{11}({\omega }_2^2)=h_{110}({\omega }_1^1)-{\omega }_2^2{h}_{111}({\omega }_1^1)\end{array}$$

$$h_{000}=c_0,h_{001}=c_1,...,h_{111}=8$$

这是一种递归，我讨厌递归，他很慢，所以我们考虑能不能把递归换成递推

递推部分：

看了这个图你一定可以懂、

规则：

当$n=8$具体操作（配合代码食用更佳）

我们还有一步就是把 $h$ 函数转回去

我们考虑这样做：

我们把 $h$ 放入一个数列

$<h({\omega }_n^1),h({\omega }_n^2),h({\omega }_n^3),...,h({\omega }_n^n)>$

把这个数列在进行一次傅里叶变换，得出这个序列的离散傅里叶级数，但是是负的

$$g({\omega }_n^{-k})=h(w_n^0)+h(w_n^1){\omega }^{-k}+...+h(w_n^{n-1}{\omega }^{-(n-1)})$$

你成功学会了FFT，当然NTT比他好十倍甚至九倍，也是一样简单

NTT

NTT

代码

FFT代码

#include<bits/stdc++.h>
#define llf double 
using namespace std;
const llf PI=acos(-1);
const int N=5e6+10;
int n,m,x,rev[N],logO=0,mn;
struct Cmop{
    llf x,y;
}a[N],b[N],c[N];
Cmop operator + (Cmop a ,Cmop b){return {a.x+b.x,a.y+b.y};}
Cmop operator - (Cmop a ,Cmop b){return {a.x-b.x,a.y-b.y};}
Cmop operator * (Cmop a ,Cmop b){return {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
void FFT(Cmop *c,int len,int op){
    for(int i = 0;i < len; ++i)if(i < rev[i]) swap(c[i], c[rev[i]]);
    for(int k = 1;k < len;k <<= 1){//当前有多少行 
        Cmop omega = {cos(PI / k),sin(PI / k * op)};
        for(int j = 0; j < len;j += (k << 1)){//j列和j+1列操作 
            Cmop o = {1,0};
            for(int i = 0; i < k; ++i){//把i行进行操作 
                Cmop u = c[i + j], v = o * c[i + j + k];
                c[i + j] = u + v;
                c[i + j + k] = u - v;
                o = o * omega;                
            }
        }
    }
}
void input(){
    scanf("%d%d", &n, &m);
    for(int i = 0;i <= n; ++i){
        scanf("%d", &x);
        a[i] = {x*1.0, 0};
    }
    for(int i = 0;i <= m; ++i){
        scanf("%d", &x);
        b[i] = {x*1.0, 0};
    }
}
void op(){
    mn=1;
    while(mn <= n+m) mn <<= 1, ++logO;
    for(int i = 0;i < mn; ++i)rev[i] = (rev[i>>1] >> 1) | (( i & 1) << (logO-1));
    FFT(a, mn, 1);
    FFT(b, mn, 1);
    for(int i = 0; i < mn; ++i){
        c[i] = a[i] * b[i];
    }
    FFT(c, mn, -1);
    for(int i = 0; i < mn; ++i){
        c[i].x /= mn;
    }
    for(int i = 0; i <=m+n; ++i){
        printf("%d ", (int)(c[i].x + 0.1));
    }
}
int main(){
    input();
    op(); 
    return 0;
}

NTT代码

#include<bits/stdc++.h>
#define ll long long 
using namespace std;
const ll N = 4e6+10,MOD = 998244353,g = 3,revG = 332748118;
int n,m,rev[N],logO=0,mn;
ll a[N], b[N], c[N];
inline ll mpow(ll a,ll k){
    ll ans = 1;
    while(k){
        if(k & 1) ans=(ans * a) % MOD;
        a=(a * a) % MOD;
        k >>= 1;
    }
    return ans%MOD;
}
inline void FFT(ll *c,int len,int op){
    for(int i = 0;i < len; ++i)if(i < rev[i]) swap(c[i], c[rev[i]]);
    for(int k = 1;k < len;k <<= 1){//当前有多少行 
        ll g_k = mpow(op == 1 ? g : revG , (MOD-1) / (k << 1));
        for(int j = 0; j < len;j += (k << 1)){//j列和j+1列操作 
            ll o = 1;
            for(int i = 0; i < k; ++i){//把i行进行操作 
                ll u = c[i + j], v = o * c[i + j + k] % MOD;
                c[i + j] = (u + v) % MOD;
                c[i + j + k] =(u - v + MOD) % MOD;
                o = o * g_k % MOD;                
            }
        }
}
inline void input(){
    scanf("%d%d", &n, &m);
    for(int i = 0;i <= n; ++i){
        scanf("%lld", &a[i]);
    }
    for(int i = 0;i <= m; ++i){
        scanf("%lld", &b[i]);
    }
}
inline void op(){
    mn=1;
    while(mn <= n+m) mn <<= 1, ++logO;
    for(int i = 0;i < mn; ++i)rev[i] = (rev[i>>1] >> 1) | (( i & 1) << (logO-1));
    FFT(a, mn, 1);
    FFT(b, mn, 1);
    for(int i = 0; i < mn; ++i){
        c[i] = a[i] * b[i] % MOD;
    }
    FFT(c, mn, -1);
    ll inv = mpow(mn, MOD-2);
    for(int i = 0; i <=m + n; ++i){
        printf("%lld ", c[i]*inv%MOD);
    }
}
int main(){
    input();
    op(); 
    return 0;
}