FFT & NTT

导言

FFT/NTT在acm竞赛中常被用于处理多项式乘法优化问题，并常常与生成函数等数学内容一同考察，笔者在此整理了一些自己对FFT的理解，希望对读者有所帮助

多项式

表达

对于一个在实数域上阶多项式，要唯一确定这样一个多项式，我们可以使用两种方法，一是如表达式一样使用序列记录的系数来唯一确定这一矩阵，记为系数表示法，二可利用插值多项式唯一性的性质通过构造个不同点来唯一确定这一多项式，记为点值表示法。

乘法运算

在这里我们重点关注n阶多项式的乘法，考虑两个阶多项式和 ,求乘积。直接采用系数表示法暴力相乘可得记，有整体时间复杂度

如果使用点值表示法，则可以在已经知道个和对应的情况下，只需要就可求出目标多项式的相关点值，从而唯一确定,此时在处值为

复数的引入

点值的选取

根据上文，点值表示法的优越在于知道点值的情况下可大幅提高运算效率，可如果暴力求点值，每个点需要次操作，进而使整体复杂度退化为，故使用点值表示法加速多项式乘法的关键在于选取一组特殊点，使得我们能够快速求取多项式在这组点上的对应点值

欧拉公式和复数的单位根

如图在复平面上，以原点为圆心，半径为构造单位圆，用条过原点直线讲圆等分并交出个点。由欧拉公式可知,对有

单位根的特殊性质

引入复平面的原因就在于单位根的一些特殊性质可提高计算效率，这里直接给出，读者往下便可体会到这些性质的作用

性质一：

性质二:

性质三：

FFT求点值

对阶多项式，按下标奇偶分为和，显然有，将、分别代入。前者得由性质一可化简得记为。后者代入得由性质二化简得记为。

由和可知求解的问题最终可以层层递归至求故可以用的复杂度求出和，最终可以用的复杂度求出所有

蝴蝶变换

考虑到每次拆分函数都会使用大量空间和时间，因此考虑优化此过程对于系数序列 $，不足则高位补$ ，手动操作FFT拆分函数的过程不难发现，对于系数，记其经过全部拆分后的位置为始终有，其中为的二进制表示字串。因此我们可直接算出最终系数使时间复杂度真正成为对于二进制长度为的数,尝试求其二进制反转后的数，显然0反转后仍为0，1反转后为,对于,有，即可递推求出全部

傅里叶逆变换（IDFT）

现在我们得到了的点值表示，需要将点值表示转化为系数表示。

使用矩阵表示FFT的过程，得到

不妨记为，要得到,只需要在等式左右两行左乘

在实际进行矩阵运算之前，不妨回到的一些性质，由性质三可知，考察公式因式分解得当时，有，即更进一步，有由此构造矩阵可得即

故对傅里叶变换后的多项式在做一次选取根为的傅里叶变换即可获得原多项式

FFT模板

NTT

由于FFT中存在复数和大量的三角函数运算且不支持取模，在整数域内取点的NTT应运而生，FFT的加速运用了复数单位根的性质，同样,NTT选取的根也应当具有类似性质。

原根

前置知识：欧拉函数,欧拉定理对一对互质正整数，且设满足的最小n为a模p的阶，记为。对所有，所有结果不同，证明:不妨设存在使则有，与阶定义矛盾

对正整数，若模的阶且，则称为的一个原根，注意为欧拉函数对正整数,整数为其原根的充要条件为且为的一组剩余化简系，由阶的性质和原根定义可直接证明

对质数，存在原根设有，

对存在如下性质：

性质一：性质二：性质三：其中性质一可直接由原根性质得出，下面给出性质二和三的证明性质二：对，由定义可得又同时，根据原根定义有故得性质三：

将原根和单位根对比，不难发现原根具有与单位根相似性质，因此原本通过代入的FFT就可以通过换成代入变成NTT，最后逆变换的过程则可以视为代入在下的逆元

关于原根的求取，一般采用直接枚举。

模板

#include <bits/stdc++.h>
typedef long long ll;
const int N = 4000000 + 100, g = 3, invg = 332748118, mod = 998244353;
//998244353的一个原根为3且998244353-1=2^23*119，3在模998244353意义下的逆元为332748118
using namespace std;
int n, m, r[N]; //rev[i]为i的二进制翻转
long long a[N], b[N];
ll qpow(ll a, ll k) //快速幂，a为底数，k为指数
{
    ll res = 1;
    while (k)
    {
        if (k & 1)
            res = res * a % mod;
        a = a * a % mod;
        k >>= 1;
    }
    return res;
}
void NTT(ll *a,int type,int len)
{
    for(int i = 0;i < len;++i)
        if(i < r[i])
            swap(a[i],a[r[i]]);
    for(int mid = 1;mid < len;mid <<= 1)
    {
        ll wn = 1;
        if(type == 1)
            wn = qpow(g,(mod - 1) / (mid << 1));
        else
            wn = qpow(invg,(mod - 1) / (mid << 1));
        for(int i = 0;i < len;i += (mid << 1))
        {
            for(int j = 0,w0 = 1;j < mid;++j,w0 = 1ll * w0 * wn % mod)
            {
                ll x = a[i + j], y = w0 * a[i + j + mid] % mod;
                a[i + j] = (x + y) % mod;
                a[i + j + mid] = (x - y + mod) % mod; 
            }
        }
    }
}
int main()
{
    string s1,s2;
    cin >> s1 >> s2;
    // cout << s1 << endl << s2 << endl;
    n = s1.length() - 1;
    m = s2.length() - 1;
    for(int i = 0;i <= n;++i)
        a[i] = s1[n - i] - '0';
    for(int i = 0;i <= m;++i)
        b[i] = s2[m - i] - '0';
    int len = 1 << max((int)ceil(log2(n + m)), 1);
    for(int i = 0;i < len;++i)
        r[i] = (r[i >> 1] >> 1) | ((i & 1) << (max((int)ceil(log2(n + m)), 1) - 1));
    NTT(a,1,len),NTT(b,1,len);
    for(int i = 0;i < len;++i)
        a[i] = a[i] * b[i] % mod;
    NTT(a,-1,len);
    ll inv = qpow(len,mod-2);
    for (int i = 0; i <= len; i++)
        a[i] = a[i] * inv % mod;
    for (int i=0; i < len; i++)
    {
        if (a[i]>=10) 
        {
            a[i+1] += a[i] / 10; 
            a[i] %= 10;
        }
    }
    while(!a[len] && len) len--;
    while(a[len] >= 10)
    {
        a[len + 1] += a[len] / 10;
        a[len++] %= 10;
    }
    for(int i = len;i >= 0;--i)
        cout << a[i];
    return 0;	
}

练习

P1919 【模板】A*B Problem 升级版（FFT 快速傅里叶变换）

P3803 【模板】多项式乘法（FFT）

2020中国大学生程序设计竞赛（CCPC） - 网络选拔赛 M - Residual Polynomial

CF1824B2