多项式学习笔记(二) NTT
前置芝士:
原根:满足 \(g^{\varphi(p)-1} \equiv 1 \pmod p\) , 且 \(g^0,g^1,g^2...g^{\varphi(p)-1}\) 互不相同,就称 \(g\) 为 \(p\) 的一个原根
NTT
在FFT的时候,我们一个比较大的问题就是精度处理。如果说题目要求系数对大质数取模的话,那么FFT可能在计算中溢出,这时候我们就需要另外一种算法 \(NTT\) (其实是 FFT 的变形)。
在FFT 中,我们单位根满足以下性质:
- \(w_{n}^{0} = w_{n}^{n} = 1\)
- \(w_{n}^{0},w_{n}^{1},w_{n}^{2}....w_{n}^{n}\) 互不相同。
- \(w_{n}^{i} = w_{2*n}^{2*i}\)
- \(w_{n}^{k+{n\over 2}} = -w_{n}^{k}\)
- \(\displaystyle\sum_{i=0}^{n-1} (w_{n}^{k})^{i} = 0\) (\(k\nmid n\))
然后我们找到了单位根一个很好的替代品原根,我们设 \(w_{n}^{i} = (g^{p-1\over n})^i\) 。带入之后可以发现它是满足以上性质的。
性质1: \(w_{n}^{0} = w_{n}^{n} = 1\)
很简单, \(w_{n}^{0} = (g^{p-1\over n}) ^0 = 1\) , \(w_{n}^{n} = (g^{p-1\over n})^n = g^{p-1} = 1\) 。
性质2: \(w_{n}^{0},w_{n}^{1},w_{n}^{2}....w_{n}^{n}\) 互不相同。
根据原根的性质可得: \(g^{0},g^1....g^{p-1}\) 都互不相同,所以 性质2得证。
性质3: \(w_{2*n}^{2*i} = w_{n}^{i}\)
\(w_{2*n}^{2*i} = (g^{p-1\over 2n}) ^{2i} = g^{{p-1\over 2n} \times 2i} = (g^{p-1\over n})^i\)
\(w_{n}^{i} = (g^{p-1\over n})^i\) .
所以 \(w_{2*n}^{2*i} = w_{n}^{i}\) ,性质三得证。
性质4: \(w_{n}^{k+{n\over 2}} = -w_{n}^{k}\)
\(w_{n}^{k+{n\over 2}} = (g^{p-1\over n})^{k+{n\over 2}} = (g^{p-1\over n})^k (g^{p-1\over n})^{n\over 2} = (g^{p-1\over n})^{k} g^{p-1\over 2}\)
有一个有意思的性质就是 \(g^{p-1\over 2} = -1\)
\(\because (g^{p-1\over 2})^2 = 1\)
\(\therefore g^{p-1\over 2} = 1 或 g^{p-1\over 2} = -1\)
\(\because g^{p-1\over 2} \neq g^{p-1}\)
\(\therefore g^{p-1\over 2} = -1\)
所以 \(w_{n}^{k+{n\over 2}} = -(g^{p-1\over n})^k = -w_{n}^{k}\)
性质5: \(\displaystyle\sum_{i=0}^{n-1} (w_{n}^{k})^i = 0\) \(k\nmid n\)
根据等比数列的求和公式原式等价于 \(1-(w_{n}^{k})^{n}\over 1-w_{n}^{k}\)
分母 \(1-(w_{n}^{k})^n = 1-(g^{p-1\over n})^{kn} = 1-(g^{p-1})^k = 0\)
性质五,得证。
NTT 和FFT 的代码其实是差不多的。
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
#define int long long
const int N = 5e6+10;
const int p = 998244353;
int n,m,len,tim,a[N],b[N],Rev[N];
inline int read()
{
int s = 0,w = 1; char ch = getchar();
while(ch < '0' || ch > '9'){if(ch == '-')w = -1; ch = getchar();}
while(ch >= '0' && ch <= '9'){s = s * 10 + ch - '0'; ch = getchar();}
return s * w;
}
int ksm(int a,int b)
{
int res = 1;
for(; b; b >>= 1)
{
if(b & 1) res = res * a % p;
a = a * a % p;
}
return res;
}
int inv(int n){return ksm(n,p-2);}
void NTT(int *a,int flag)
{
for(int i = 0; i < len; i++)
{
if(i < Rev[i]) swap(a[i],a[Rev[i]]);
}
for(int h = 1; h < len; h <<= 1)
{
int wn = ksm(3,(p-1)/(h<<1));//原根
if(flag == -1) wn = inv(wn);
for(int j = 0; j < len; j += (h<<1))
{
int w = 1;
for(int k = 0; k < h; k++)
{
int u = a[j+k];
int v = w * a[j+h+k] % p;
a[j+k] = (u + v) % p;
a[j+h+k] = (u-v+p) % p;
w = w * wn % p; //wni
}
}
}
}
signed main()
{
n = read(); m = read();
for(int i = 0; i <= n; i++) a[i] = read();
for(int i = 0; i <= m; i++) b[i] = read();
len = 1, tim = 0;
while(len <= n+m) len <<= 1, tim++;
for(int i = 0; i < len; i++) Rev[i] = (Rev[i>>1]>>1) | ((i&1)<<(tim-1));
NTT(a,1); NTT(b,1);
for(int i = 0; i < len; i++) a[i] = a[i] * b[i] % p;
NTT(a,-1);
for(int i = 0; i <= n+m; i++) printf("%lld ",a[i]*inv(len) % p);
printf("\n");
return 0;
}