FFT 快速傅里叶变换 与 NTT 快速数论变换 学习笔记

这个东西是个玄乎的东西，反正最后还是背板子 理论知识也不用管太多吧

这里搬一些简单的证明

复数运算法则 设复数$A = a + bi, B = c + di$

$A + B = (a + c) + (b + d)i$
$A - B = (a - c) + (b - d)i$
$A × B = (ac - bd) + (ad + bc)i$
除法转化为乘共轭复数就好了

一般系统有自带的$std::complex<double>$可以直接用
但是手写速度要有保证一些

首先整篇文章讨论的多项式的次数都是$2^k - 1$次的，以下不作说明

$ω_n^k = cos(k\frac{2π}{n}) + i×sin(k\frac{2π}{n})$ 这个根据复数乘法可知

$ω_n^k = ω_{2n}^{2k}$ 这个根据上面那条东西可以得到

$ω_n^{k + \frac{n}{2}} = -ω_n^k$ 代入第一个式子得出$ω^\frac{n}{2} = -1$即可

那么我们假设原多项式是$A(x) = \sum_{i = 0}^{n - 1}a_ix^i$

奇偶分组后两个多项式分别为
$A_1(x) = a_0 + a_2 x + ... + a_{n - 2}x^{\frac{n - 1}{2}}$
$A_2(x) = a_1 + a_3 x + ... + a_{n - 1}x^{\frac{n - 1}{2}}$
则$A(x) = A_1(x^2) + x×A_2(x^2)$

我们试着推导以下已知$A_1, A_2$怎么还原出$A$
设$k < \frac{n}{2}$, 则$k和k+\frac{n}{2}$可以取遍$[0, n - 1]$的所有值
$A(ω_n^k) = A_1(ω_n^{2k}) + ω_n^k×A_2(ω_n^{2k})$
$\ \ \ \ \ \ \ \ \ \ =A_1(ω_\frac{n}{2}^k) + ω_n^k×A_2(ω_\frac{n}{2}^k)$
$A(ω_n^{k + \frac{n}{2}}) = A_1(ω_n^{2k+n}) + ω_n^{k+\frac{n}{2}}×A_2(ω_n^{2k+n})$
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ =A_1(ω_\frac{n}{2}^k) - ω_n^k×A_2(ω_\frac{n}{2}^k)$
具体推式子只需要之前推出的单位复根的性质即可
这样子就可以分治下去直到一个常数项为止，我们得到了一个在$O(nlogn)$时间把多项式从系数转化为点值的算法

一般点值表示法不是很常用，所以我们还要考虑怎么把一个多项式从点值转换为系数

设$c_k = \sum_{i = 0}^{n - 1}y_i(ω_n^{-k})^i$
也就是多项式 $C(x) = \sum_{i = 0}^{n - 1}y_ix^i$ 在$x$取${ω_n^{-k}}时的点值$
其中 $y_i$表示$x取{ω_n^{i}}时的A(x)的值$
下面开始推式子

$c_k = \sum_{i = 0}^{n - 1}y_i(ω_n^{-k})^i$
$c_k = \sum_{i = 0}^{n - 1} (ω_n^{-k})^i×(\sum_{j = 0}^{n - 1}a_j(ω_n^i)^j)$
$c_k = \sum_{i = 0}^{n - 1} \sum_{j = 0}^{n - 1}(ω_n^{j - k})^ia_j$
$c_k = \sum_{j = 0}^{n - 1} a_j\sum_{i = 0}^{n - 1}(ω_n^{j - k})^i$
我们发现$\sum_{i = 0}^{n - 1}(ω_n^{j - k})^i$是一个等比数列求和的形式
带入高中数学和蔼可亲的公式$S_n = a_1\frac{1 - q^n}{1 - q}$
我们发现公比$q = ω_n^{j - k}，a_1 = 1$
则$S_n = \frac{1 - (ω_n^{j - k})^n}{1 - q}$
当分母不为$0$时， $S_n = 0$， 否则$S_n = n$
则$c_k = a_kn$
即$a_k = \frac{c_k}{n}$
我们可以通过之前得到的$O(nlogn))$计算多项式点值的方法来实现
最后记得除以$n$就好了

总体就是快速将系数转为点值为了$O(n)$实现两个多项式相乘
最后再快速转为系数，通过单位复根优秀的性质来实现

贴个递归的板子

Codes

#include<bits/stdc++.h>
#define For(i, a, b) for(register int i = a; i <= b; ++ i)

using namespace std;

const int maxn = 4e6 + 10;
const double Pi = acos(-1);
int n, m;

struct Complex {
    double x, y;
    Complex (double xx = 0, double yy = 0) {
        x = xx, y = yy;
    }
    Complex operator + (const Complex &T) {
        return Complex{x + T.x, y + T.y};
    }
    Complex operator - (const Complex &T) {
        return Complex{x - T.x, y - T.y};
    }
    Complex operator * (const Complex &T) {
        return Complex{x * T.x - y * T.y, x * T.y + y * T.x};
    }
}a[maxn], b[maxn];

void FFT(int limit, Complex *a, int fh) {
    if(limit == 1) return;
    Complex a1[limit >> 1], a2[limit >> 1];
    for(int i = 0; i <= limit; i += 2)  
        a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];    
    FFT(limit >> 1, a1, fh), FFT(limit >> 1, a2, fh);
    Complex wn1 = Complex(cos(2.0 * Pi / limit), sin(2.0 * Pi / limit) * fh), w = Complex(1, 0);
    for(int i = 0; i < (limit >> 1); ++ i, w = w * wn1) 
        a[i] = a1[i] + w * a2[i], a[i + (limit >> 1)] = a1[i] - w * a2[i];
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("FFT.in", "r", stdin);
    freopen("FFT.out", "w", stdout);
#endif
    int x, limit = 1; scanf("%d%d", &n, &m);
    For(j, 0, n) scanf("%d", &x), a[j].x = x;
    For(j, 0, m) scanf("%d", &x), b[j].x = x;

    while(limit <= n + m) limit <<= 1;
    FFT(limit, a, 1), FFT(limit, b, 1);
    For(i, 0, limit - 1) a[i] = a[i] * b[i];
    FFT(limit, a, -1);

    For(i, 0, n + m) printf("%d ", int(a[i].x / limit + 0.5)); 

    return 0;
}

我们发现递归的效率不是很高 考虑迭代实现
盗用网上的一张图

发现最后每一位上数的下标为原下标的二进制位反转
对于$0, 1, 2, 3, 4, 5, 6, 7$
翻转后就是$0, 4, 2, 6, 1, 5, 3, 7$
把这几个数都写成二进制形式规律就很清楚了

那么迭代实现就直接加上

void FFT_Init(int limit, int *a) {
    int k = log2(limit);
    for(int j = 0; j < limit; ++ j) {
        int t = 0;
        for(int i = 0; i < k; ++ i)
            if((1 << i) & j)
                t |= 1 << (k - j - 1);
        if(t < j) swap(a[t], a[j]);
    }
}

我们来考虑一种更为优秀的基于递推的$O(n)$预处理

k = log2(limit);
for(int i = 0; i < limit; ++ i)
    rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (k - 1));

如果这个数是偶数，把它最后一位去掉后翻转，就相当于后面多了一个$0$所以删去这个$0$即可
如果是奇数那第一位是$1$，所以补上就好了

上$FFT$板子

#include<bits/stdc++.h>
#define For(i, a, b) for(register int i = a; i <= b; ++ i)

using namespace std;

const int maxn = 4e6 + 10;
const double Pi = acos(-1);

struct Complex {
    double x, y;
    Complex(double xx = 0, double yy = 0) {
        x = xx, y = yy;
    }
    Complex operator + (const Complex &T) {
        return Complex(x + T.x, y + T.y);
    }
    Complex operator - (const Complex &T) {
        return Complex(x - T.x, y - T.y);
    }
    Complex operator * (const Complex &T) {
        return Complex(x * T.x - y * T.y, x * T.y + y * T.x);
    }
}a[maxn], b[maxn];

void FFT(int n, Complex *a, int fh) {
    int k = log2(n);
    For(i, 0, n - 1) {
        int t = 0;
        For(j, 0, k - 1)
            if((1 << j) & i)
                t |= (1 << (k - j - 1));
        if(i < t) swap(a[i], a[t]);
    }
    for(int limit = 2; limit <= n; limit <<= 1) {
        double xita = 2.0 * Pi / limit;
        Complex Wn1 = Complex(cos(xita), sin(xita) * fh), W = Complex(1, 0);
        for(int j = 0; j < n; j += limit, W = Complex(1, 0))    
            for(int i = j; i < j + (limit >> 1); ++ i, W = W * Wn1) {
                Complex a1 = a[i], a2 = W * a[i + (limit >> 1)];
                a[i] = a1 + a2, a[i + (limit >> 1)] = a1 - a2;
            }
    }
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("3803.in", "r", stdin);
    freopen("3803.out", "w", stdout);
#endif
    int n, m, limit = 1, x;
    scanf("%d%d", &n, &m);
    while(limit <= (n + m))
        limit <<= 1;
    For(i, 0, n) 
        scanf("%d", &x), a[i].x = x;
    For(i, 0, m) 
        scanf("%d", &x), b[i].x = x;
    FFT(limit, a, 1), FFT(limit, b, 1);
    For(i, 0, limit - 1)
        a[i] = a[i] * b[i];
    FFT(limit, a, -1);
    For(i, 0, n + m)
        printf("%d ", int(a[i].x / limit + 0.5));
    return 0;
}

从$FFT$到$NTT$

　NTT 快速数论变换

听到这个名字就知道$NTT$这个是在数论中使用的
$FFT$已经很优秀了，那我们为什么还需要这样的东西代替它呢
因为$FFT$在复数运算的时候会有精度误差，不能保证得到答案的精确性
这个时候我们就想能不能找到一个东西代替复数的运算
这个时候$NTT$就横空出世了，首先$NTT$针对的是模意义下的运算

那么我们考虑怎么找到一个跟单位复根有着相同性质的东西
那就是原根 (本文讨论的模数为$p = qn + 1(n = 2^k)$)

原根的定义
设 $m$ 是正整数， $a$ 是整数，若 $a$ 模 $m$ 的阶等于 $\phi(m)$ ，则称 $a$ 为模 $m$ 的一个原根

阶的定义
若$a,p$互素，且$p>1$，对于$a^n≡1(mod\ p)$最小的n，我们称之为$a$模$p$的阶,记做$δ_p(a)$，例如$δ_7(2) = 3$，因为$2^3≡1(mod\ 7)$

这样子我们就可以发现原根这个东西是不唯一的

并且原根有一个很重要的定理我不会证
若$p$为素数，假设一个数$g$是$p$的原根，那么$g^i\ mod\ p(1<g<p,0<i<p)$的结果两两不同,这样子就满足了可以代不同的$x$得到多项式的点值了
上面这个东西我在网上找了很多资料也没有证明 可能可以问一下数竞的选手应该会

1.令$\omega_n = g^q$ 则$\omega_n^0, \omega_n^1, ..., \omega_n^{n - 1}$互不相同

2.$\omega_n^k = \omega_{2n}^{2k}$ 代入第一个式子即可

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

$\omega_n^{n} = g^{qn} = g^{p - 1}≡1(mod\ p)$(费马小定理可知)

$(\omega_n^\frac{n}{2})^2≡1(mod\ p)$

又因为$\omega_n^\frac{n}{2}≠1(mod\ p)$

所以$\omega_n^\frac{n}{2}≡-1(mod\ p)$

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

这样子原根就有了单位复根的所有性质了，那么就可以直接用原根代替单位复根，模意义下运算代替复数运算，就得到了可爱的$NTT$了
还有要注意的几个地方 一般$NTT$模的数是$998244353$和$1004535809$
这两个数都有一个共同的原根为$3$，$OI$中用$3$就好了,如果你还想模其他的数那么就多用几个模数做$NTT$，用中国剩余定理合并就好啦
每次分治时用的$\omega_n = g^{\frac{mod - 1}{n}}$
我的理解是为了保证$\omega_n^n = 1$因为推$\omega_n^{k + \frac{n}{2}} = -\omega_n^k$有用到过
当然逆运算的时候不要忘记把$g$换成$inv_g$就好啦

Codes

#include<bits/stdc++.h>

using namespace std;

const int N = 4e6 + 10;
const int mod = 998244353;
int a[N], b[N];

int qpow(int a, int x) {
    int ret = 1;
    for(; x; a = 1ll * a * a % mod, x >>= 1)
        if(x & 1)
            ret = 1ll * ret * a % mod;
    return ret;
}

void NTT_Init(int limit, int *a, int k) {
    for(int i = 0; i < limit; ++ i) {
        int t = 0;
        for(int j = 0; j < k; ++ j)
            if((1 << j) & i)
                t |= 1 << (k - j - 1);
        if(t < i) swap(a[t], a[i]);
    }
}

void NTT(int n, int *a, int fh) {
    int k = log2(n), inv = qpow(n, mod - 2);
    NTT_Init(n, a, k);
    for(int limit = 2; limit <= n; limit <<= 1) {

        int Wn1 = qpow(3, (mod - 1) / (limit)), w = 1;
        if(fh == -1) Wn1 = qpow(Wn1, mod - 2);

        for(int j = 0; j < n; j += limit, w = 1)
            for(int i = j; i < j + (limit >> 1); ++ i, w = 1ll * w * Wn1 % mod) {
                int a1 = a[i], a2 = 1ll * a[i + (limit >> 1)] * w % mod;
                a[i] = (a1 + a2) % mod, a[i + (limit >> 1)] = (a1 - a2 + mod) % mod;        
            }
    }
    if(fh == -1)
        for(int i = 0; i < n; ++ i)
            a[i] = 1ll * a[i] * inv % mod;
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("34.in", "r", stdin);
    freopen("34.out", "w", stdout);
#endif
    int n, m, limit = 1;
    scanf("%d%d", &n, &m);
    while(limit <= (n + m))
        limit <<= 1;
    for(int i = 0; i <= n; ++ i)
        scanf("%d", &a[i]);
    for(int i = 0; i <= m; ++ i)
        scanf("%d", &b[i]);
    NTT(limit, a, 1), NTT(limit, b, 1);
    for(int i = 0; i < limit; ++ i)
        a[i] = 1ll * a[i] * b[i] % mod;
    NTT(limit, a, -1);
    for(int i = 0; i <= n + m; ++ i)
        printf("%d ", a[i]);
    return 0;
}

终于写完了$ovo$ 这种东西对我这种蒟蒻选手来说真的理解起来不容易，感谢$menci$写的博客

引用来源
menci’s 从傅里叶变换到数论变换
menci’s FFT 学习笔记