快速傅里叶变换|快速数论变换

FFT——快速傅里叶变换

什么是FFT

$FFT$全称($Fast \ Fourier \ Transformation$)，中文名：快速傅里叶离散变换。

这个名字听起来很高级，实际上也很高级，它是$DFT$和$IDFT$的加速版本，用于加速多项式乘法。

• 接下来我们把某个多项式记为“大写字母+元素”，如多项式：

$A(x)=a_nx^n+a_{n-1}x^{n-1}+…+a_2x^2+a_1x^1+a_0$

$B(x)=b_mx^m+b_{m-1}x^{m-1}+…+b_2x^2+b_1x^1+b_0$

你也可以理解为$A/B$是关于$x$的函数

现在，我们要求$A(x)*B(x)$，即两个多项式的乘积(卷积)，那么我们很容易想到的就是$O(N^2)$的做法：把每一项依次分别相乘再相加。

$$\sum^{n+m-1}_{i=0}(\sum^{i}_{j=0}a_i*b_{i-j})*x^i$$

这样的时间复杂度明显不尽人意，那么$FFT$就是把这种算式计算优化到$NlogN$的算法。

---

前置知识

• 复数

我们定义虚数 $i$， 且 $i^2 = -1$，即 $i = \sqrt {-1}$，那么复数就形如 $a+bi$，$i$ 是一个虚数单位，$a$ 叫做复数的实部， $b$ 叫做复数的虚部，那么 $a+bi$ 在复平面（高中会讲）可以类似于平面直角坐标系，表示为点 $(a, b)$。

复数的运算法则

1. 加法：实部和虚部分别相加。

2. 减法：实部和虚部分别相减。

3. 乘法：像一次多项式一样相乘。

4. 除法：分子分母同乘分母的共轭。

共轭：两个复数实部相同，虚部相反。

上面除法分子分母同乘分母的共轭，把分母有理化成实数。

    struct complex { //复数 
        long double a, b;
        complex() { a = b = 0;}
        //重载复数四则运算 
        complex operator + (const complex x) const {
            complex temp;
            temp.a = a + x.a;
            temp.b = b + x.b;
            return temp; 
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.a = a - x.a;
            temp.b = b - x.b;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.a = a * x.a - b * x.b;
            temp.b = a * x.b + b * x.a;
            return temp; 
        }

        complex operator / (const complex x) const {
            double div = x.a * x.a + x.b * x.b;
            complex temp;
            temp.a = (a * x.a + b * x.b) / div;
            temp.b = (b * x.a - a * x.b) / div;
            return temp; 
        }
    };

这边大家直接记一个重要结论：

复数相乘，辐角相加，模长相乘。

---

单位根

对于复数$w$,若整数 $n$ 满足 $w^n=1$，则称 $w$ 是 $n$ 次单位根。

扔个单位圆：

为什么扔个单位圆呢？因为单位圆上任意一点到原点的距离，即模长都为$1$。

那么，对于一个复数，只有模长为$1$的复数才有可能成为$n$次单位根!

现在，我们在复平面上把复数锁定在了这个单位圆上。

设复数$w$模长为$q$，其幅角为$mπ$，$m$为$[0,2]$间实数，那么，其$n$次方幅角为$nmπ$。

则有$nmπ=kπ,k∈N$，即$nm=k$。

1. 若$k=0$,因为$n>0$，所以方程有且只有$m=0$一解。

2. 若$k>0$，则$k=1,2,3,…$，易知$m$有$n-1$个不重解

综上，复数$w$的$n$次单位根有$n$个。

并且，复数的$n$次单位根$n$等分单位圆。

单位根的书写：

单位根用符号$ω$表示，从1开始，沿单位圆逆时针把单位根从$0$开始标号。

然后单位根还有个神奇的东西$ω^{k}_{n}=ω_{n}^{k\mod n}$

单位根的性质

1. 对于任意$n$，$ω_n^0=1$。

2. $ω_n^k=(ω_n^1)^k$。

3. $(ω_n^m)^k=(ω_n^k)^m$

4. $ω_n^i*ω_n^j=ω_{n}^{i+j}$。

5. $ω_{2n}^{2m}=ω_n^m$。

6. 若 $2 \mid n$，$ω_n^{m+\frac{n}{2}}=-ω_n^m$

---

多项式的表达方法

1. 系数表达法：这个如同我们上面写的多项式即为系数表达法。

2. 点值表达法

定义多项式$F(x)$：

我们知道，$1$个n元一次方程最少需要$n$个式子求解，而我们把系数表达$F(x)$抽象为一个函数$F(x)$，那么这个函数至少需要$n+1$个点确定下来，那么我们任取$n+1$个不同的数构成集合$U={k_1,k_2,…,k_{n+1}}$，对$F(x)$分别求值得到一个新的集合，那么将两个集合合并成点集，有：

$$T={(k_1,F(k_1)),k_2,F(k_2),…,(k_{n+1},F(n+1))}$$

这便是$F(x)$的点值表达式，并且，点值表达法下的乘法运算更加简单：

多项式$P,Q$分别取点$(x,y_1)$和$(x,y_2)$，则$P*Q$可得到点$(x,y_1*y_2)$。令$S=P*Q$，则$C(x_i)=y_{1_i}*y_{2_i}$。

可见，点值表达法下，多项式乘法是$O(N)$的。

• （补充）乘法本质：卷积

根据我们上面暴力计算$A(x)*B(x)$，不妨设$C=A*B$，那么：

$C[k]=\sum^{k}_{i=0}A[i]B[k-i]$

$A[i],B[k-1]$分别为$A$中$x^i$的系数和$B$中$x^{k-i}$的系数，那么我们知道：$x^i*x^{k-i}=x^k$，故所有的$x^i*x^{k-i}$的系数对$x$的次项系数均有贡献，需要累加。

那么就可以写成：$C[k]=\sum\limits_{i+j=k}A[i]B[j]$

那么形如$C[k]=\sum\limits_{i⊕j=k}A[i]B[j]$($⊕$为某种运算)的式子称为卷积。

那么，多项式的乘法就是对于加法的卷积。

我们仍然可以对多项式相乘的模板打个暴力：

    #include<bits/stdc++.h>
    using namespace std;

    const int SIZE = 1e6 + 10;

    long long A[SIZE], B[SIZE], C[SIZE];
    int n, m;

    int main() {
        scanf("%d %d", &n, &m);
        n++, m++;
        for (int i = 0; i < n; i++) scanf("%lld", A + i);
        for (int i = 0; i < m; i++) scanf("%lld", B + i);
        for (int k = 0; k < n+m-1; k++)
            for (int i = 0; i <= k; i++)
                C[k] += A[i] * B[k - i];
        for (int i = 0; i< n + m - 1; i++) 
            printf("%lld ", C[i]);
        return 0;
    }

模板题亲测 44 分。

---

算法理论

上面的主要部分还是写了一个$O(N^2)$的暴力，那我们要寻找一个更快的办法：

就是它！$FFT$算法

前面我们已经说过：$FFT$算法是一个$O(NlogN)$的算法，但常数较大(相信你们会卡常)，并且，$FFT$是$DFT$和$IDFT$的快速算法：

---

思想：

首先你要会点值表达式，当然我们已经讲过。

然后我们来看$A(x)=x^2+3x-2$（红色）与$B(x)=-x^2+6x+4$(绿色)的取点得到点的点值表达式:

在$A(x)$上取三个点（紫色）$(-4,2)(-2,-4)(1,2)$；在$B(x)$上取三个点$(0,4)(2,12)(6,4)$，然后转为系数表达$O(N^2)$相乘然后再带入$x$得到$C$的点值表达式？这时绝对不行的，我们上面说过，点值表达式相乘可以在$O(n)$内解决，那怎么办呢？

我们可以取$3$个$x$相同的点：

这样就行了？那是不可能的。我们知道，两个二次项式相乘得到的是一个四次项式，而$C(x)=y_1*y_2$，只能得到三个点，这肯定是不行的。

这很明显我们需要$2n+1$个点，不够怎么办呢？再添两个不就好了？
下面添加了$I,J$(横坐标相同)，$G,H$(横坐标相同)，这样就行了。

那么，我们只需要进行$2n+1$次乘法即可，省略常数（都说了常数巨大了），复杂度$O(N)$

神仙问题：$tql$ $but$ 这有用吗

当然有用啊：

算法流程：

$$\text{系数表达式}\implies \text{点值表达式} \implies \text{系数表达式}$$

这不就有个用了吗？

但这个思路仅仅是$DFT+IDFT$，而$FFT$是用来加速$DFT$和$IDFT$的。

如果让我们求点值表达式，那我们肯定会带整数进去，因为这样好计算。但是傅里叶，（不知道是不是脑抽了），把复数带进了多项式中，然后神奇的事情就发生了……

我们把多项式如下拆开：

$F(x)=(a_0+a_2x^2+…+a_{n-2}x^{n-2})+(a_1x^1+a_3x^3+…+a_{n-1}x^{n-1})$

（保证$n=2^k,k∈N^*$)

然后设两个有$n/2$次项的多项式$FL(x)$和$FR(x)$

$$FL(x)=a_0+a_2x+…+a_{n-2}x^{\frac{n}{2}-1}$$

$$FR(x)=a_1+a_3x+…+a_{n-1}x^{\frac{n}{2}-1}$$

$$F(x)=FL(x^2)+xFR(x^2)$$

设$k<\frac{n}{2}$，把$ω_n^k$代入$F(x)$：$F(ω_n^k)=FL(ω_n^{2k})+ω_n^kFR(ω_n^{2k})$

接着，就有$F(ω_n^k)=FL(ω_{n/2}^{k})+ω_n^kFR(ω_{n/2}^{k})$

那么，如果我们知道$FL(x),FR(x)$分别在$ω^0_{n/2},ω^1_{n/2},…,ω^{n/2-1}_{n/2}$处的点值表示，那么我们就可以$O(N)$内求出$F(x)$在$ω^0_{n},ω^1_{n},…,ω^{n/2-1}_{n}$处的点值表示了。

接下来我们的目标就是求$F(x)$在$ω^{n/2}_{n},ω^{n/2+1}_{n},…,ω^{n-1}_{n}$处的点值表示。

然后……

继续设$k<n/2$，然后把$ω^{k+n/2}_{n}$代入$F(x)$中。

运用一系列单位根的性质，读者可以自行证明出：

$$F(ω^{k+n/2}_{n})=FL(ω^{k}_{n/2})-ω^{k}_{n}FR(ω^{k}_{n/2})$$

然后我们就可以求出$F(x)$在$ω^0_{n},ω^1_{n},…,ω^{n-1}_{n}$处的点值表示，接着，我们就实现了$DFT$（听着好累的样子……）

但是我们还不知道$FL(x),FR(x)$在$ω^0_{n},ω^1_{n},…,ω^{n/2-1}_{n}$处的点值表示。难道要暴力代入吗？这样反而会使时间复杂度上升。

但是我们想想：$FL(x)$和$FR(x)$是由原来问题$F(x)$转化来的，怎么转化的？看到这里，你应该想到了分治。我们知道，分治的时间复杂度是$logN$，好像完成了一半了？

然后，因为$FL(x),FR(x)$是$F(x)$的子问题，那么我们可以对它们继续进行分解成一个个小的子问题，然后再合并，类似于归并排序。

$FFT$分治的合并过程根据：

1. $F(ω^{k}_{n})=FL(ω^{k}_{n/2})+ω^{k}_{n}FR(ω^{k}_{n/2})$

2. $F(ω^{k+n/2}_{n})=FL(ω^{k}_{n/2})-ω^{k}_{n}FR(ω^{k}_{n/2})$

从理论到实现

上文的所有证明都基于$2=2^k,k∈N^*$，但是我们做题的时候并没有这一点，所以我们可以补一些系数为$0$的项，这对计算结果没有影响。

---

1. 预处理单位根

我们知道一个性质$ω_n^k=(ω_n^1)^k$

那么，我们只要知道$ω_n^1$就可以求出$ω_n^0,ω_n^1,ω_n^2,…,ω_n^{n-1}$。

怎么求$ω_n^1$？

解三角形不就好了？已知$|OA|=1$,幅角是$\frac{1}{n}$圆周，那么：

(由于$c++$下三角函数用弧度制表示，以下全部使用弧度制。)

设幅角$α=\frac{2π}{n}$,，则$ω_n^1(cos\big(\frac{2π}{n} \big),sin\big(\frac{2π}{n} \big))$。

这边要牢记，$x$是实轴，$y$是虚轴：

    struct complex { //复数 
        long double a, b;
        complex() { a = b = 0;}
            //重载复数四则运算 
        complex operator + (const complex x) const {
            complex temp;
            temp.a = a + x.a;
            temp.b = b + x.b;
            return temp; 
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.a = a - x.a;
            temp.b = b - x.b;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.a = a * x.a - b * x.b;
            temp.b = a * x.b + b * x.a;
            return temp; 
        }

        complex operator *= (const complex &x) {
            *this = *this * x;
            return *this;
        }

        complex operator / (const complex x) const {
            double div = x.a * x.a + x.b * x.b;
            complex temp;
            temp.a = (a * x.a + b * x.b) / div;
            temp.b = (b * x.a - a * x.b) / div;
            return temp; 
        }
    }w[SIZE],temp, buf;

    void unit_roots(int n) {
        const long double pi = 3.14159265358979;
        //complex ;
        temp.a = cos(2*pi/n), temp.b = sin(2*pi/n);
        w[0].a = 1, w[0].b = 0;
        buf.a = 1, buf.b = 0;
        for(int i = 1; i < n ; i++) {
            buf *= temp;
            w[i] = buf; 
        }
        for(int i = 0; i < n; i++)
            printf("w_%d %.3Lf %.3Lf\n", i, w[i].a, w[i].b);	
    } 

然后我们就可以实现$DFT$了

因为还没讲$IDFT$，先写$DFT$:

递归$DFT$实现（复数板子上面有了下面就不写了）

    int n, m;
    Complex f[SIZE];
    void dft(Complex *f, int len) { //当前子问题长度
        if(len == 0) return;
        Complex fl[len + 10], fr[len + 10];
        for(int i = 0 ; i < len; i++)
            fl[i] = f[i << 1], fr[i] = f[i << 1 | 1];
        dft(fl, len >> 1);
        dft(fr, len >> 1);
        len *= 2;
        Complex temp, buf;
        temp.a = cos(2*pi/n), temp.b = sin(2*pi/n);
        buf.a = 1;
        for(int i = 0 ; i < len / 2; i++) {
            f[i] = fl[i] + buf * fr[i];
            f[i + len / 2] = fl[i] - buf * fr[i];
            buf *= temp;
        }
    }
    int main() {

        scanf("%d", &n);
        for(int i = 0; i < n; i++) scanf("%Lf", &f[i].a);
        for(m = 1; m < n; m <<= 1);//补到2的幂次方
        dft(f, m >> 1);
        for(int i = 0; i < m; i++) printf("%.3Lf ", f[i].a);
        return 0;
    }

事实上虚数可以用(小心TLE)：

            complex <double> f;

然后我们得到了一堆点值表达……

这并没有什么用

于是，我们现在来学习$IDFT$——$DFT$的逆运算。

刚刚的多项式为$F(x)$，系数是$A$数列，点值表示成$M$，$M[k]=\sum \limits_{i=0}^{n-1}(ω_n^{k})^i*F[i]$

然后就有$F[k]=\frac{1}{n}\sum \limits_{i=0}^{n-1}(ω_n^{-k})^i*M[i]$

结论记住就好了，可以自己尝试证明。

我们可以先求出和式，然后再除以$n$。

接下来我们求出$ω_n^0,ω_n^{-1},…ω_n^{-n+1}$

并且仍有：$ω_n^{-j}=(ω_n^{-1})^j$

所以我们还是求出$ω_0^{-1}$即可，并且，$ω_0^{-1}$与$ω_0^{1}$关于$x$轴对称。

    #include<iostream>
    #include<cstdio>
    #include<cmath>
    using namespace std;
    const long double pi = 3.141592653589793;
    const int SIZE = 4e6 + 10;
    struct complex { //复数 
        long double real, imag;
        complex() { real = imag = (double)0;}
            //重载复数四则运算 
        complex operator + (const complex x) const {
            complex temp;
            temp.real = real + x.real;
            temp.imag = imag + x.imag;
            return temp; 
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.real = real - x.real;
            temp.imag = imag - x.imag;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.real = real * x.real - imag * x.imag;
            temp.imag = real * x.imag + imag * x.real;
            return temp; 
        }

        complex operator *= (const complex &x) {
            *this = *this * x;
            return *this;
        }

        complex operator / (const complex x) const {
            double div = x.real * x.real + x.imag * x.imag;
            complex temp;
            temp.real = (real * x.real + imag * x.imag) / div;
            temp.imag = (imag * x.real - real * x.imag) / div;
            return temp; 
        } 
        void conj() { imag = -imag; } 
    }a[SIZE], b[SIZE];

    complex omega(int n, int k) {
        complex temp;
        temp.real = cos(2 * pi * k / n);
        temp.imag = sin(2 * pi * k / n);
        return temp;
    }	
    int n, m, p;		
    void fft(complex *f, int len, int inv) {
        if(len == 1) return;
        int mid = len >>1;
        complex fl[mid + 10], fr[mid + 10];
        for(int i = 0; i < mid; i++)
            fl[i] = f[i << 1], fr[i] = f[i << 1 | 1];
        fft(fl, mid, inv);
        fft(fr, mid, inv);
        complex temp, buf; 
        temp.real = cos(2 * pi / len);
         temp.imag = sin(2 * pi / len);
        buf.real = 1, buf.imag = 0;
        if(inv) temp.conj();
        for(int i = 0; i < mid; i++, buf *= temp) {
            complex s = buf * fr[i];
            f[i] = fl[i] + s;
            f[i + mid] = fl[i] - s;
        }
    }

    int main() {

        scanf("%d %d", &n, &m);
        for(p = 1; p < n + m + 1; p <<= 1);
        for(int i = 0; i <= n; i++) 
            scanf("%Lf", &a[i].real);
        for(int i = 0; i <= m; i++) 
            scanf("%Lf", &b[i].real);
        fft(a, p, 0);
        fft(b, p, 0);
        for(int i = 0; i < p; i++) 
            a[i] *= b[i];
        fft(a, p, 1);
        for(int i = 0; i <= n + m; i++)
            printf("%d ", (int)(a[i].real / p + 0.5));
	return 0;
}

接着TLE了……(丝毫没有卡常的欲望)

因为递归常数巨大，并且很容易爆空间。(但实际上是因为写了$long \ double$炸了)

下面是AC的递归代码，以后不要写long double 了

    #include<iostream>
    #include<cstdio>
    #include<cmath>
    using namespace std;
    const double pi = 3.141592653589793;
    const int SIZE = 4e6 + 10;
    struct complex { //复数 
        double real, imag;
        complex() { real = imag = (double)0;}
            //重载复数四则运算 
        complex operator + (const complex x) const {
            complex temp;
            temp.real = real + x.real;
            temp.imag = imag + x.imag;
            return temp; 
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.real = real - x.real;
            temp.imag = imag - x.imag;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.real = real * x.real - imag * x.imag;
            temp.imag = real * x.imag + imag * x.real;
            return temp; 
        }

        complex operator *= (const complex &x) {
            *this = *this * x;
            return *this;
        }

        complex operator / (const complex x) const {
            double div = x.real * x.real + x.imag * x.imag;
            complex temp;
            temp.real = (real * x.real + imag * x.imag) / div;
            temp.imag = (imag * x.real - real * x.imag) / div;
            return temp; 
        } 

        void conj() { imag = -imag;}
    }a[SIZE], b[SIZE];

    complex omega(int n, int k) {
        complex temp;
        temp.real = cos(2 * pi * k / n);
        temp.imag = sin(2 * pi * k / n);
        return temp;
    }	
    int n, m, p;		
    void fft(complex *f, int len, int inv) {
        if(len == 1) return;
        int mid = len >>1;
        complex fl[mid + 10], fr[mid + 10];
        for(int i = 0; i < mid; i++)
            fl[i] = f[i << 1], fr[i] = f[i << 1 | 1];
        fft(fl, mid, inv);
        fft(fr, mid, inv);
        complex temp, buf; 
        temp.real = cos(2 * pi / len);
         temp.imag = sin(2 * pi / len);
        buf.real = 1, buf.imag = 0;
        if(inv) temp.conj();
        for(int i = 0; i < mid; i++, buf *= temp) {
            complex s = buf * fr[i];
            f[i] = fl[i] + s;
            f[i + mid] = fl[i] - s;
        }
    }

    int main() {

        scanf("%d %d", &n, &m);
        for(p = 1; p < n + m + 1; p <<= 1);
        for(int i = 0; i <= n; i++) 
            scanf("%lf", &a[i].real);
        for(int i = 0; i <= m; i++) 
            scanf("%lf", &b[i].real);
        fft(a, p, 0);
        fft(b, p, 0);
        for(int i = 0; i < p; i++) 
            a[i] *= b[i];
        fft(a, p, 1);
        for(int i = 0; i <= n + m; i++)
            printf("%d ", (int)(a[i].real / p + 0.5));
        return 0;
    }

实际上跑的很快的递归评测记录。

那有没有更优化的办法呢？这是肯定的。

---

蝴蝶变换

先看一张图

这张图是什么？这不是我们分治的顺序吗？

我在一开始和结束都打上了二进制，你有没有什么发现？

是的，每一个数的二进制都反过来了。

然后我们的任务就转换成求二进制反序：

    rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << l)//l是二进制位数
    //rev[i]是i翻转得到的

结论牢记即可。

还记得$DP$吗？（貌似没什么关系）

我们用循环实现$FFT$，也就像$DP$一样定义状态、阶段决策：

• 第一层循环枚举变换区间大小（阶段），第二层循环枚举开头（决策），第三层处理区间。

    #include<iostream>
    #include<cstdio>
    #include<cmath>
    using namespace std;
    const double pi = 3.141592653589793;
    const int SIZE = 4e6 + 10;
    struct complex { //复数 
        double real, imag;
        complex() { real = imag = (double)0;}
            //重载复数四则运算 
        complex operator + (const complex x) const {
            complex temp;
            temp.real = real + x.real;
            temp.imag = imag + x.imag;
            return temp; 
        }

        complex operator += (const complex &x) {
            *this = *this + x;
            return *this;
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.real = real - x.real;
            temp.imag = imag - x.imag;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.real = real * x.real - imag * x.imag;
            temp.imag = real * x.imag + imag * x.real;
            return temp; 
        }

        complex operator *= (const complex &x) {
            *this = *this * x;
            return *this;
        }

        complex operator / (const complex x) const {
            double div = x.real * x.real + x.imag * x.imag;
            complex temp;
            temp.real = (real * x.real + imag * x.imag) / div;
            temp.imag = (imag * x.real - real * x.imag) / div;
            return temp; 
        } 
    }a[SIZE], b[SIZE];
    int n, m, p, L = -1;
    int rev[SIZE];

    void fft(complex *f, int inv) {
        for(int i = 0 ; i < p ; i++) 
            if(i < rev[i]) 
                swap(f[i], f[rev[i]]);
        for(int len = 2; len <= p; len <<= 1) {
            int length = len >> 1;
            complex temp;
            temp.real = cos(pi / length);
            temp.imag = sin(pi / length);
            if(inv) temp.imag = -temp.imag;
            for(int l = 0; l < p; l += len) {
                complex buf; buf.real = 1;
                for(int k = l; k < l + length; k++) {
                    complex t = buf * f[k + length];
                    f[k + length] = f[k] - t;
                    f[k] += t;
                    buf *= temp;
                }
            }
        }
    }

    int main() {
        scanf("%d %d", &n, &m);
        for(int i = 0; i <= n; i++)
            scanf("%lf", &a[i].real);
        for(int i = 0; i <= m; i++)
            scanf("%lf", &b[i].real);
        for(p = 1; p < n + m + 1; p <<= 1, L++);
        //l是位数，因为会多算一位，一开始初值直接给-1
        for(int i = 1; i < p; i++)
            rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L);
        fft(a, 0), fft(b, 0);
        for(int i = 0; i < p; i++)
            a[i] *= b[i];
        fft(a, 1);
        for(int i = 0; i <= n + m; i++)
            printf("%d ", (int)(a[i].real / p + 0.5));
        return 0;
    } 

评测记录

NTT——快速数论变换

• 首先，学习$NTT$算法你要学会$FFT$（理解也可）（刚刚我们才讲完你不会忘了吧……）（我认为我写的很详细了呢$QAQ$）。

• 接下来，你要确保你会原根这个东西（起码要知道），不过我在后面还会再提，具体可以先看我的这一篇博客：同余相关，在比较下面的部分（由于还待完善，所以只是写了一点）。

• 然后就是希望大家看了这一篇博客对$NTT$有所了解，可能我讲的不好，但是起码你要知道$NTT$是什么。

---

前置知识

（怎么又是前置知识）

$FFT$

上面已经说了先理解$FFT$，并且我也给了写的比较详细的博客，这里不再展开。

但其实$FFT$并不是重点，而是里面用于优化的很重要的一个东西：$\color{red}{\text{单位根}}$。

---

• 原根

但是如果 $FFT$ 要求取模呢？

$$\color{red}{What?\text{你不会？}}$$

好像因为 $FFT$ 用的单位根不太好取模（复数），所以我们可能需要找一个东西来替代下单位根，那么就要求这个东西和单位根具有类似的性质。然后科学家们就开始寻找了……（也不知道是谁找到的）接下来他们找到了原根。

---

首先我们要了解一下啥是原根。（如果你看了上面给出的博客你应该知道了）

原根的定义

其实在那一篇博客里的原根定义很容易让人摸不着头脑……为啥呢？

$$\color{green}{δ_m(a)=φ(m) \text{？} what \ is \ this?}$$

因为它很不好求……那我们不妨换一种说法：

设 $p$ 为质数，则根据 $\color{orange}{\text{费马小定理}}$，有：

$$\text{对于任意 }p\nmid a,a^{p-1}\equiv1 \pmod{p}$$

$$\color{red}{\text{那么称 $g$ 为模 $p$ 的原根，当且仅当}\{g^0,g^1,…,g^{p-2}\}\text{对} \ p \ \text{的模数互不相同。}}$$

这样可能会更好理解些？至于为什么是上面那个绿绿的，我也不知道啊……

• 补充：本原单位根

其实定义很简单，对于一些比较特殊的$n$次单位根，它的$0,1,…,n-1$次幂恰好生成了所有的$n$次单位根，我们把这样的单位根称为本原单位根。

在复数域 $C$ 上，$\omega_n=\exp\frac{2\pi i}{n}=\cos\frac{2\pi}{n}+i\sin\frac{2\pi}{n}$ 是一个本原单位根。

（以上不是重点，下面这句话才是重点）

在有限域（即模质数）上，本原单位根和数论中的原根有关！

这也是我补充本原单位根的原因了……因为后面我们还会提到它。

并且我们可以证明的是，模质数的原根总是存在的（不知道为啥你就问度娘吧）。

并且原根的性质和单位根非常相似哦~

$NTT$ 里原根的性质

换句话说，在模 $p$ 意义下，$g$ 可以被看做一个 $p-1$ 次本原单位根。

反正感性理解就好，结论不就是用来记的嘛……

设 $n\mid p-1$ ,我们不妨令 $\omega_n=g^{\frac{p-1}{n}}$（这个其实很像单位根的求法），然后你就可以得到一大串东西了：（好像也不多）

$$\omega_n^n=(g^{\frac{p-1}{n}})^n=g^{p-1}\equiv1\pmod p$$

并且因为单位根的定义，你会知道 $\omega_n^0,\omega_n^0,…,\omega_n^{n-1}$ 互不相同……这不就完了嘛 $QAQ$

因此此时的 $\omega_n$ 在模 $p$ 意义下具有 $n$ 次本原单位根的性质，所以我们利用这个东西来定义 $DFT$ 即可，这被称为 $\color{red}{NTT}$。

原根的求法

求模质数原根的办法： 将 $p-1$ 分解质因数，即 $p-1=p_1^{a_1}p_2^{a2}…p_{n}^{a_n}$ 为标准的$p-1$分解式，那么若 $g^{\frac{p-1}{p_i}}\ne 1\pmod p$ ，则$g$是模 $p$ 意义下的一个原根。

求模合数原根的办法： 其实你只需要把上面的 $p-1$ 换成 $φ(p)$ 即可。

证明：（假设法）
  先对于第一种模质数的情况：
  假设存在一个 $t<φ(p)=p-1$ 使得 $g^t\equiv1\pmod p$ 且 $\forall i∈[1,n],g^{\frac{p-1}{p_i}}\ne 1 \pmod p$。由裴蜀定理可知：一定存在一组 $(a,b)$ 使得：$a*t=(p-1)*b+gcd(t,p-1)$ 成立。由欧拉定理/费马小定理：$g^{p-1}\equiv 1 \pmod p$。所以 $g^{a*t}\equiv g^{(p-1)*b+gcd(t,p-1)}\equiv1\pmod p$。因为 $t<p-1$，所以又有 $gcd(t,p-1)<p-1$。又因为 $gcd(t,p-1)\mid p-1$，于是 $\exists i ∈[1,n],gcd(t,p-1)\mid{g^{\frac{p-1}{p_i}}}$。那么就可以得到 $g^{\frac{p-1}{p_i}}\equiv g^{gcd(t, p-1)}\equiv 1 \pmod p$。
 因此假设不成立。所以原命题成立。
$Q.E.D.$

那么上述结论自然可以推广到模合数而把 $p-1$ 改为 $φ(p)$ 了。

以上的证明是针对于这一句话的：$δ_m(a)=φ(m)$，也就是上述求法成立的条件。我们从另一个性质出发，即 $g^0,g^1,…,g^{p-1}$ 模$p$的余数各不相同，那么可以得知这些余数的集合为：$\{1,2,…p-1\}$（不分顺序）。我们知道取模是满足乘法性质的。也就是说，如果$\exists i∈[1,n]$使得$g^{\frac{p-1}{p^i}}=1$，我们知道唯有 $g^0 \bmod p=1 \bmod p = 1$ ，那么很明显 $\frac{p-1}{p^i}$ 是不会为 $0$ 的，那么一旦出现第一个数字的重复，也就是 $1$，那么很明显不满足原根的上述性质。

那么我就很好奇，为什么重复的$1$一定会出现在某个$g^{\frac{p-1}{p^i}}$当中呢？上面其实也已经证明，如果存在某个 $t<φ(p)$ 使得 $g^t\equiv1\pmod p$ 成立，那么必定会存在一个 $g^\frac{p-1}{p^i}\bmod p=1$ 。并且我认为这样可能更好理解。

反正证明这种东西感性理解一下就好了……

然后原根你也会求了，所以用原根 $g$ 去代替 $FFT$ 中的 $\omega_n$ (即单位根),就可以完成 $FFT$ 了。

NTT 模数

首先根据 $FFT$ 我们知道 $g^{\frac{p-1}{n}}$中的 $n$ 是 $2$ 的幂次，所以在做 $FFT$ 的时候我们最好构造形如 $p=a*2^k+1$的质数，这样就可以满足上面的需求了。这样的质数最好取大一点（精度），所以这样的质数可以取：

$p=a*2^k+1$	$\ $ $a$	$k$	$g$
$3$	$1$	$1$	$2$
$5$	$1$	$2$	$2$
$17$	$1$	$4$	$3$
$97$	$3$	$5$	$5$
$193$	$3$	$6$	$5$
$257$	$1$	$8$	$3$
$7681$	$15$	$9$	$17$
$12289$	$3$	$12$	$11$
$40961$	$5$	$13$	$3$
$65537$	$1$	$16$	$3$
$786433$	$3$	$18$	$10$
$5767169$	$11$	$19$	$3$
$7340033$	$7$	$20$	$3$
$23068673$	$11$	$21$	$3$
$104857601$	$25$	$22$	$3$
$167772161$	$5$	$25$	$3$
$469762049$	$7$	$26$	$3$
$\color{red}{998244353}$	$119$	$23$	$3$
$\color{red}{1004535809}$	$479$	$21$	$3$
$1998585857$	$953$	$21$	$3$
$2013265921$	$15$	$27$	$31$
$2281701377$	$17$	$27$	$3$
$3221225473$	$3$	$30$	$5$
$75161927681$	$35$	$31$	$3$
$77309411329$	$9$	$33$	$7$
$206158430209$	$3$	$36$	$22$
$2061584302081$	$15$	$37$	$7$
$2748779069441$	$5$	$39$	$3$
$6597069766657$	$3$	$41$	$5$
$39582418599937$	$9$	$42$	$5$
$79164837199873$	$9$	$43$	$5$
$263882790666241$	$15$	$44$	$7$
$1231453023109121$	$35$	$45$	$3$
$1337006139375617$	$19$	$46$	$3$
$3799912185593857$	$27$	$47$	$5$
$4222124650659841$	$15$	$48$	$10$
$7881299347898369$	$7$	$50$	$6$
$31525197391593473$	$7$	$52$	$3$
$180143985094819841$	$5$	$55$	$6$
$1945555039024054273$	$27$	$56$	$5$
$4179340454199820289$	$29$	$57$	$3$

理论到实现

第一件事，我们先来回忆一下 $FFT$，而 $FFT$ 分为两个步骤：$DFT$ 与 $IDFT$。

然后我们知道单位根的定义是对于复数域 $C$，有 $z^n=1$ 的复数就是一个 $n$次单位根，$n$ 次单位根有 $n$ 个，为：$e^{\frac{2\pi k}{n}i}, \ k = \{0,1,2…,n-1\}$（$i$是虚数单位）。（你用三角表示也是可以的，但是不好理解）

然后这里就用 $g^{\frac{p-1}{n}}$ 代替 上面那个就可以，因为你会发现他们具有相同的性质。

唯一具有卷积性质的变换就是 $DFT$，按照上面的式子，$DFT$ 的变换式就是：

$$X_k=\sum_{n=0}^{N-1}x_ne^{-\frac{2\pi}{N}nki} \ \ \ \ \ \ k=0,1,2…,N-1$$

$DFT$ 反演（也就是逆变换）公式有：

$$x_n=\frac{1}{N}\sum_{k=0}^{N-1}X_ke^{\frac{2\pi}{N}nki} \ \ \ \ \ \ n=0,1,2…N-1$$

那么由于 $g^{\frac{p-1}{n}}$ 和 $e^{\frac{2\pi k}{n}i}, \ k = \{0,1,2…,n-1\}$ 等价，所以：

$$e^{-\frac{2\pi}{N}nki}\equiv g^{\frac{p-1}{N}}\pmod p$$

接着我们就得到了快速数论变换的公式：

$$X_k=\sum_{n=0}^{N-1}x_ng^{nk}\pmod p \ \ \ \ \ \ k=0,1,2…,N-1$$

反演：

$$x_n=\frac{1}{N}\sum_{k=0}^{N-1}X_kg^{-nk}\pmod p\ \ \ \ \ \ n=0,1,2…N-1$$

这样我们就成功把复数域转到了整数域，然后剩下的东西在$\mod p$意义下考虑即可。上面也已经讲过素数的构造方法。

• 补充：原根为何可以代替单位根？

考虑单位根的性质（也就是我们为什么使用单位根）：

1. $n$ 个单位根互不相同。那么$g^0,g^1,…,g^{p-2}$ 在模 $p$ 意义下也不相同，成立。

2. $z^n=1$（$z∈C$）。根据费马小定理：$g^{p-1}\equiv 1\pmod p$，成立。

3. 单位根对称分布。其实这对于原根也成立。

证明：
  $∵(g^{\frac{p-1}{2}})^2\equiv g^{p-1}\equiv 1 \pmod p$
  又$∵g^i \ne g^j(i\ne j \text{且}0≤i,j≤p-1)$
  $∴g^{\frac{p-1}{2}}= \ne 1,\text{则}g^{\frac{p-1}{2}}=-1$
$Q.E.D.$

好了，我们真的可以使用 $NTT$ 来完成 $FFT$ 所能完成和所不能完成的事了。

可以尝试过一下模板题：$P3803$【模板】多项式乘法($FFT$)

虽然题目没有要求模数，但是仍然可以取一个比较大的模数从而完成。

这里取 $p=998244353,g=3$。

我们结合代码来理解 $NTT$ :

#include<bits/stdc++.h>
using namespace std;
const int N = (1 << 21) + 100, P = 998244353, G = 3;

inline int read() {
	static int x, f;
	static char c;
	while(!isdigit(c = getchar()) && c != '-');
	x = (f = c == '-') ? 0 : c - '0';
	while(isdigit(c = getchar())) 
		x = (x << 3) + (x << 1) + (c ^ 48);
	return f ? -x + 1 : x;
}

typedef long long ll;
inline ll power(ll a, ll k) {
	ll base = 1;
	for(; k; k >>= 1) {
		if(k & 1) base = (base * a) % P;
		a = (a * a) % P;
	}	
	return base;
}

int rev[1 << 21], limit = 1;
ll a[N], b[N];
inline void NTT(ll *f, int inv) {
	for(int i = 0; i < limit; i++)
		if(i < rev[i]) swap(f[i], f[rev[i]]);
	
	for(int len = 1; len < limit; len <<= 1) {
		ll temp = power(G, (P - 1) / (len << 1));
		if(inv) temp = power(temp, P - 2);
		//反NTT时，除法全部变为乘逆元 
		for(int l = 0; l < limit; l += len * 2) {
			ll w = 1;
			for(int k = 0; k < len; k++) {
				int x = f[l + k], y = w * f[l + k + len] % P;
				f[l + k] = (x + y) % P;
				f[l + k + len] = (x - y + P) % P;
				w = (w * temp) % P;
			}
		}
	} 
}

int main() {
	int n = read(), m = read(), L = -1;
	for(int i = 0; i <= n; i++) a[i] = (read() + P) % P;
	for(int i = 0; i <= m; i++) b[i] = (read() + P) % P;
	for( ; limit <= n + m; limit <<= 1, L++); 
	for(int i = 0; i < limit; i++)
		rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L);
	NTT(a, 0), NTT(b, 0);
	for(int i = 0; i < limit; i++)
		a[i] = (a[i] * b[i]) % P;
	NTT(a, 1);
	ll inv = power(limit, P - 2);
	for(int i = 0; i <= n + m; i++)
		printf("%lld ", (a[i] * inv) % P);
	
	return 0;
}

然后真的比$FFT$快好多呀……模板题评测记录：

$\text{递归}FFT$：点我

$\text{蝴蝶变换}FFT$：点我

$NTT$：点我

---

Upd20200918:文章当中还有一些错误，近期本人正在矫正。

__EOF__

本文作者：Ning-H
本文链接：https://www.cnblogs.com/Ning-H/p/12072626.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本文版权归作者和博客园共有，欢迎转载，但未经作者同意必须在文章页面给出原文连接，否则保留追究法律责任的权利。
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！