多项式算法——快速数论变换NTT
在 FFT 中,我们使用复数计算,未免会出现精度损失,基于原根的快速数论变换 NTT 解决了这个问题。
模意义下的 n n n次单位根
在 FFT 中,我们计算出多项式在 ω n = 1 \omega^n = 1 ωn=1 的 n n n 次复数域下的 n n n 个单位根下的点值来完成 FFT ,这是因为 n n n次单位根具有下面的三个性质:
-
周期性: n n n 个单位根互不相同,且幂次具有周期性。
-
消去引理:
ω d n d k = e 2 d k π d n i = e 2 k π n i = ω n k \omega^{dk}_{dn} = e^{\frac{2 dk \pi}{dn}i} = e^{\frac{2 k \pi}{n}i} = \omega^{k}_{n} ωdndk=edn2dkπi=en2kπi=ωnk
-
折半引理:
如果 n > 0 n \gt 0 n>0是偶数,那么 n n n个 n n n次单位复数根的平方的集合,等同于 n / 2 n/2 n/2个 n / 2 n/2 n/2次单位复数根的的集合。
根据消去引理 ( ω n k ) 2 = ω n 2 k (\omega_n^k)^2=\omega^k_{\frac{n}{2}} (ωnk)2=ω2nk,如果对每个根都进行平方,那么每个不同的数正好出现两次,因为:
( ω n k + n / 2 ) 2 = ω n 2 k + n = ω n 2 k ω n n = ω n 2 k = ( ω n k ) 2 = ω n / 2 k (\omega_n^{k + n/2})^2 = \omega_n^{2k+n} = \omega_n^{2k}\omega_n^n = \omega_n^{2k} = (\omega_n^k)^2 = \omega_{n/2}^k (ωnk+n/2)2=ωn2k+n=ωn2kωnn=ωn2k=(ωnk)2=ωn/2k
我们根据这几条性质,定义模 p p p (我们只讨论素数模下的NTT)意义下的 n n n 次主单位根。
设 g g g 为 p p p 的一个原根,那么 g 0 , g 1 , … , g p − 1 g^0,g^1,\ldots,g^{p-1} g0,g1,…,gp−1 和 0 , 2 , … , p − 1 0,2,\ldots,p-1 0,2,…,p−1 构成双射关系,也就是 p p p 的一个简化剩余系,并且周期出现,满足第一条性质。
令 g n = g p − 1 n g_n = g^{\frac{p-1}{n}} gn=gnp−1 ,称 g n g_n gn 为 模 p p p 意义下的 n n n 次主单位根。
快速数论变换NTT
有了 g n g_n gn 为 模 p p p 意义下的 n n n 次主单位根,则对于消去引理:
g d n d k = g ( p − 1 ) k d d n = g n k g^{dk}_{dn} = g^{\frac{(p-1)kd}{dn}} = g^{k}_{n} gdndk=gdn(p−1)kd=gnk
对于折半引理:
( g n k + n / 2 ) 2 = g n 2 k + n = g n 2 k g n n = g n 2 k = ( g n k ) 2 = g n / 2 k (g_n^{k + n/2})^2 = g_n^{2k+n} = g_n^{2k}g_n^n = g_n^{2k} = (g_n^k)^2 = g_{n/2}^k (gnk+n/2)2=gn2k+n=gn2kgnn=gn2k=(gnk)2=gn/2k
在这里 g n n ≡ 1 m o d p g_n^n \equiv 1 \mod p gnn≡1modp 可由费马小定理得到。
对比负数域上的单位圆, 2 π n \frac{2\pi}{n} n2π 为将整个圆等分,而 p − 1 n \frac{p-1}{n} np−1 将剩余系等分。
故我们只需要替换 w n w_n wn 为 g n g_n gn 就是快速数论变换NTT。
对于逆快速数论变换NTT,也类似逆FTT,在逆FTT中,乘以单位根的倒数,并乘以 n n n 的倒数。
a j = 1 n ∑ k = 0 n − 1 y k ω n − k j a_j = \frac{1}{n}\sum_{k=0}^{n-1}y_k \omega_n^{-kj} aj=n1k=0∑n−1ykωn−kj
则对于 NTT ,我们乘以 g n j g_n^j gnj 的逆元和 n n n 的逆元即可。
要注意的几点为:
- NTT 仍存在计算损失,只不过不是浮点误差,当系数超过 p p p 的时候,那么得到的结果是模完 p p p 之后的结果。
- 对于 p p p 和 g g g 的选择,我们令 p p p 尽量是 2 2 2 的幂次的倍数加一,可以选择 p = 998244353 , g = 3 p = 998244353,g=3 p=998244353,g=3。
代码
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
#ifdef LLT_DBG
#define FR freopen("in.txt", "r", stdin)
#else
#define FR
#endif
template <ll P>
ll fpow(ll a, ll b)
{
ll res = 1;
for (; b; b >>= 1, a = (a * a) % P)
if (b & 1)
res = (res * a) % P;
return res;
}
template <ll G, ll P>
struct NTT
{
int _n;
int E;
vector<int> rev;
/**
* @brief 构建一个 NTT 计算器
*
* @param n 多项式最高项数
*/
NTT(int n)
{
_n = 1;
E = 0;
while (_n < n)
{
_n <<= 1;
E++;
}
rev.resize(_n);
// 逆位置对换
for (int i = 1; i < _n; i++)
{
rev[i] = (rev[i >> 1] >> 1) + ((i & 1) << (E - 1));
}
}
void _rNTT(ll A[], ll k)
{
for (int i = 0; i < _n; i++)
if (i < rev[i])
swap(A[i], A[rev[i]]);
for (int e = 1; e <= E; e++)
{
int m = 1 << e;
for (int i = 0; i < _n; i += m)
{
int hf = m / 2;
ll g = 1;
ll gn = fpow<P>(fpow<P>(G, (P - 1) / m), k);
for (int j = 0; j < hf; j++)
{
ll x = A[i + j];
ll y = (A[i + j + hf] * g) % P;
A[i + j] = (x + y) % P;
A[i + j + hf] = (x - y) % P;
g = (g * gn) % P;
}
}
}
}
/**
* @brief NTT 过程
*
* @param A 系数数组
*/
void doNTT(ll A[])
{
_rNTT(A, 1);
}
/**
* @brief NTT 逆过程
*
* @param A 点值数组
*/
void doINTT(ll A[])
{
ll ni = fpow<P>(_n, P - 2);
_rNTT(A, P - 2);
for (int i = 0; i < _n; i++)
{
A[i] = (A[i] * ni) % P;
A[i] = (A[i] + P) % P;
}
}
};
ll A[5000005];
ll B[5000005];
const int mod = 998244353;
void solve()
{
int n, m;
cin >> n >> m;
for (int i = 0; i <= n; i++)
{
cin >> A[i];
}
for (int i = 0; i <= m; i++)
{
cin >> B[i];
}
NTT<3, mod> ntt(n + m + 1);
ntt.doNTT(A);
ntt.doNTT(B);
for (int i = 0; i < ntt._n; i++)
{
A[i] = (A[i] * B[i]) % mod;
}
ntt.doINTT(A);
for (int i = 0; i < n + m + 1; i++)
{
cout << A[i] << " ";
}
}
int main()
{
FR;
solve();
return 0;
}