多项式
系数表示法与点值表示法
系数表示法：\(f(x)=a_0x^0+a_1x^1+a_2x^2+a_3x^3+a_4x^4+a_5x^5\)
点值表示法：一个 \(n\) 次多项式由 \(n+1\) 个点 \((x_i,f(x_i))\) 表示。
为何 \(n+1\) 个点可以确定一个 \(n\) 次多项式？
写成矩阵形式：
当 \(x_i\) 确定且互不相同时，中间那个矩阵是范德蒙德矩阵，矩阵可逆，所以点值表示法和系数表示法可以相互转换。
FFT
FFT 的基本思想
- 系数表示法 \(\rightarrow\) 点值表示法（DFT）
- 将两个多项式点值相乘（点积）
- 点值表示法 \(\rightarrow\) 系数表示法（IDFT）
点的选择
我们当然可以随意选择 \(n+m+1\) 个点，然后代入求值，但是那样是 \(n^2\) 的，我们需要寻找一些特殊的点来简化计算。
我们选取 \(k=2^p\) 个点，满足 \(2^p\geq n+m+1\)，点 \(x_i={w_k}^i\)，其中 \({w_k}^1\) 表示复平面上将单位圆 \(k\) 等分后逆时针数第 \(i\) 个的向量。
一些需要了解的数学知识，可以通过复数的运算或者在坐标系上旋转证明（下面公式中所有 \(i\) 均为虚数单位）
\(e^{i\pi}=-1\)（欧拉公式）
\(e^{\theta\pi}=\cos\theta+i\sin\theta\)
\({w_k}^p=\cos({p\over k}\times 2\pi)+i\sin({p\over k}\times 2\pi)\)
\({w_k}^{p+{k\over 2}}=-{w_k}^p\)
若 \(w_k^p=\cos\theta+i\sin\theta\)，那么\({1\over w_k^p}=\cos\theta-i\sin\theta\)
求出点值
先看看一个点 \(x\) 如何求值：
奇偶分离
设 \(g(x)=a_0x^0+a_2x^1\)，\(h(x)=a_1x^0+a_3x^1\)
对于多个点的情况，有：
所以我们可以通过计算 \(g(x),h(x)\) 来得到 \(f(x)\)，因此复杂度优化为 \(O(n\log n)\)。
非递归实现
逐层处理，递推出 rev 数组表示反转二进制位后的结果。
IDFT
将 \(f(x)=g(x^2)+xh(x^2)\) 改为 \(f(x)=g(x^2)+{1\over x}h(x^2)\)，最后把结果除以 \(2^p\) 即可。
NTT
原根
设原根为 \(g\)，模数为 \(p\)。
-
原根的性质
对于 \(i\in [0,p),g^i \bmod p\) 互不相同
-
原根的判定
\(\gcd(g,p)=1\)，且 \(g^i\equiv 1\) 的最小解为 \(\varphi(p)\)（欧拉定理）
-
原根存在定理
原根存在当且仅当 \(p=2,4,c^a,2c^a\)，其中 \(c\) 为奇质数，\(a\in N^{*}\)
-
判定方法
找出 \(\varphi(p)\) 的质因数 \(p_i\)，如果所有 \(p_i\) 都满足 \(g^{\varphi(p)\over p_i}\not\equiv 1\)，那么 \(g\) 是原根。
-
寻找最小原根
最小的原根在 \(O(n^{0.25})\) 级别，直接暴力即可。
单位根的改变
和 FFT 类似，NTT 也需要寻找单位根。
设质数 \(p=a\times 2^b+1\)，原根为 \(g\)。
我们发现 \(g^{(p-1)\over k}\) 具有和 \(w_k\) 相同的性质，可以拿来做单位根。
将 FFT 中所有运算在剩余系下进行即可。
常见模数
\(998244353=479\times 2^{21}+1,g=3\)（注意：该模数能够处理的最大数据范围是 \(n+m\leq 2\times 10^6\)）
\(1004535809=7\times 17\times 2^{23}+1,g=3\)
多项式求逆
采用倍增实现，递归公式为：
\(g\prime(x)\) 为新的答案数组，次数为 \(n-1\)，\(g(x)\) 为原来的答案数组。
常数
计算常数时应当把 \(n\) 视为 \(2^{\lceil \log n\rceil}\)，常数按照取模次数算。
FFT \(n=10^6,m=10^6\) 用时 \(1s\)
NTT 折算 \(4.5\) 倍常数（取 \((n+m)\log(n+m)\) 的倍数计算），\(n=2\times 10^6,m=2\times 10^6\) 用时 \(1s\)
多项式求逆 折合 \(27\) 倍常数，\(n=10^6\) 用时 \(1s\)。
NTT 和爆乘比较
\(n\times n \rightarrow 2n\)，\(n=180\) 时二者相当。
\(n\times n \rightarrow n\)，\(n=360\) 时二者相当。（使用最新的优化，可以做到 \(n=500\) 时相当）
\(n\times n \rightarrow 2n\) 和 三模 NTT，\(n=450\) 时二者相当。
\(n^{1.584}\) 乘法 和 三模 NTT，\(n=2048\) 时二者相当。
\(n^2\) 的多项式操作
https://www.cnblogs.com/tzcwk/p/dxs-sqr.html
代码
- FFT
#include<bits/stdc++.h>
using namespace std;
#define int long long
struct cp{
double x,y;
};
cp operator +(const cp&n1,const cp&n2){return {n1.x+n2.x,n1.y+n2.y};}
cp operator -(const cp&n1,const cp&n2){return {n1.x-n2.x,n1.y-n2.y};}
cp operator *(const cp&n1,const cp&n2){return {n1.x*n2.x-n1.y*n2.y,n1.y*n2.x+n1.x*n2.y};}
cp operator /(const cp&n1,const int&x){return {n1.x/x,n1.y/x};}
const int N=2.2e6+5;
const double pi=acos(-1);
int n,m;
cp a[N],b[N],c[N];
cp pw[N];
int rev[N];
void dft(cp *a,int op,int n){
for(int i=0;i<n;i++){
if(i<rev[i]) swap(a[i],a[rev[i]]);
}
for(int k=1;k*2<=n;k<<=1){
for(int i=0;i<k;i++){
pw[i]={cos(i*pi/k),sin(i*pi/k)};
if(op) pw[i].y=-pw[i].y;
}
for(int i=0;i<n;i+=k*2){
for(int j=0;j<k;j++){
cp g=a[i+j],h=pw[j]*a[i+j+k];
a[i+j]=g+h;
a[i+j+k]=g-h;
}
}
}
if(op){
for(int i=0;i<n;i++) a[i]=a[i]/n;
}
}
void fft(cp *a,cp *b,cp *c,int n,int m){
int P=1,D=0;
while(P<n+m+1) P<<=1,D++;
for(int i=0;i<P;i++){
rev[i]=(rev[i>>1]>>1)|((i&1)<<(D-1));
}
dft(a,0,P);dft(b,0,P);
for(int i=0;i<P;i++) c[i]=a[i]*b[i];
dft(c,1,P);
}
signed main(){
ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
cin>>n>>m;
for(int i=0;i<=n;i++){
int x;cin>>x;a[i]={x,0};
}
for(int i=0;i<=m;i++){
int x;cin>>x;b[i]={x,0};
}
fft(a,b,c,n,m);
for(int i=0;i<=n+m;i++) cout<<(int)(c[i].x+0.5)<<' ';
return 0;
}
- NTT
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int mod=998244353;
const int N=1.1e6+5;
int qpow(int a,int x){
int sum=1;
while(x){
if(x&1) sum=sum*a%mod;
a=a*a%mod;
x>>=1;
}return sum;
}
inline int Add(int x){return x>=mod? x-mod: x;}
inline int Dec(int x){return x<0? x+mod: x;}
inline int Ba(int x){
const static __int128 bar=(__int128(1)<<64)/mod;
return x-(x*bar>>64)*mod;
}
namespace NTT{
const int g[2]={3,qpow(3,mod-2)};
int h[N],rev[N];
int P,D;
void dft(int a[],int op,int n){
for(int i=0;i<n;i++){
if(i<rev[i]) swap(a[i],a[rev[i]]);
}
for(int k=1;k*2<=n;k<<=1){
int bs=qpow(g[op],(mod-1)/k/2);
h[0]=1;
for(int i=1;i<k;i++) h[i]=h[i-1]*bs%mod;
for(int i=0;i<n;i+=k*2){
for(int j=i;j<i+k;j++){
int A=a[j],B=Ba(h[j-i]*a[j+k]);
a[j]=Add(A+B);
a[j+k]=Dec(A-B);
}
}
}
if(op){
int iv=qpow(n,mod-2);
for(int i=0;i<n;i++) a[i]=(a[i]*iv%mod+mod)%mod;
}
}
void ntt(int a[],int b[],int c[],int n,int m){
for(int i=0;i<=n;i++) a[i]=(a[i]%mod+mod)%mod;
for(int i=0;i<=m;i++) b[i]=(b[i]%mod+mod)%mod;
P=1;D=0;
while(P<=n+m+1) P<<=1,D++;
for(int i=0;i<P;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(D-1));
dft(a,0,P);dft(b,0,P);
for(int i=0;i<P;i++) c[i]=a[i]*b[i]%mod;
dft(c,1,P);
}
}
int a[N],b[N],c[N];
int n;
signed main(){
NTT::ntt(a,b,c,n,n);
return 0;
}
- \(n^{1.584}\) 乘法
const int S=4e7+5;
int PP[S];
int *make(int n){
static int now=0;
if(n<0){now=0;return 0;}
memset(PP+now,0,sizeof(int)*n);
now+=n;
assert(now<S);
return PP+now-n;
}
void cdq(int *a,int *b,int *c,int n){
static unsigned long long cc[100];
if(n<=16){
for(int i=0;i<n*2;i++) cc[i]=0;
for(int i=0;i<n;i++){
for(int j=0;j<n;j++) cc[i+j]+=a[i]*b[j];
}
for(int i=0;i<n*2;i++) c[i]+=cc[i]%mod;
return;
}
int m=n/2;
int *f3=make(m),*f4=make(m);
int *g1=make(n),*g3=make(n);
for(int i=0;i<m;i++) f3[i]=Add(a[i]+a[i+m]),f4[i]=Add(b[i]+b[i+m]);
cdq(a,b,g1,m);cdq(f3,f4,c+m,m);cdq(a+m,b+m,g3,m);
for(int i=0;i<n;i++) c[i]+=g1[i],c[i+n]+=g3[i],c[i+m]+=-g1[i]-g3[i];
}
void cdq(){
for(int i=0;i<n*2;i++) c[i]=0;
make(-1);
cdq(a,b,c,P);
for(int i=0;i<n*2;i++) c[i]=(c[i]%mod+mod)%mod;
}
浙公网安备 33010602011771号