多项式杂谈(更新中)
前言
最近在学习多项式,稍微开个坑吧。
一些比较好的学习资料、模板:
我个人习惯写指针,在我的模板中基本没有 vector
,可能是我觉得 vector
常数大吧。
本文以存板子为主,一些边界条件可以参照洛谷模板题的输入格式。
快速傅里叶变换 FFT
题目大意:给定两个多项式 \(f(x),g(x)\),求 \(f(x) * g(x)\)。即求:\(h[i]=\sum f[j]\times g[i-j]\)。
这题显然有一个 \(O(n^2)\) 的暴力,不过我们有一个 \(O(n\log n)\) 的优(du)秀(liu)算法,常数大就是了。
一些前置知识:复数,单位根的性质。
我们知道 \(n\) 个点可以确定一个 \(n-1\) 次多项式,即用 \(n\) 个点 \((x_i,y_i)\) 来表示一个多项式。我们称之为点值表示法。而两个点值多项式乘起来只需要将对应点值相乘即可,即 \((x_i,y1_i), (x_i,y2_i)\to (x_i,y1_i\times y2_i)\),时间复杂度是 \(O(n)\) 的。
注意 \(n\) 次多项式和 \(m\) 次多项式相乘要选 \(n+m+1\) 个点值(多项式的系数是 \(n+m\) 次的)。
接下来我们的问题就是如何将系数多项式和点值多项式互相转化。我们将系数 \(\to\) 点值的过程称为 DFT
(离散傅里叶变换),将点值 \(\to\) 系数的过程称为 IDFT
(离散傅里叶逆变换)。
如果我们暴力带入点值,时间复杂度还是 \(O(n^2)\) 的,不可接受。但是我们可以通过带入一些特殊的点值来优化这个过程。下面不妨假设 \(n\) 为 \(2\) 的次幂。设一个多项式 \(A(x)\)。
\[ {\displaystyle} A(x)=\sum _{i=0} ^ {n-1} a_ix^i \]
将 \(x^i\) 按 \(i\) 的奇偶性分成两部分,并将后一部分的 \(x\) 提出来,有:
\[ {\displaystyle} A(x)=\sum _{i=0} ^{\lfloor \frac n2 \rfloor-1} a_{2i}x^{2i}+x\sum _{i=0} ^{\lfloor \frac n2 \rfloor-1} a_{2i+1}x^{2i} \]
再设两个多项式 \(A_1(x),A_2(x)\)。
\[A_1(x)=a_0x^0+a_2x^2+\cdots \\ A_2(x)=a_1x^0+a_3x^2+\cdots \\ A(x)=A_1(x^2)+xA_2(x^2) \]
设 \(k< \frac n2\),将 \(\omega_n^k\) 和 \(\omega_n^{k+n/2}\) 带入。
\[A(\omega_n^k)=A_1((\omega_n^k)^2)+\omega_n^kA_2((\omega_n^k)^2) \\ =A_1(\omega_{n/2}^k)+\omega_n^kA_2(\omega_{n/2}^k) \\ A(\omega_n^{k+n/2})=A_1((\omega_n^{k+n/2})^2)+\omega_n^{k+n/2}A_2((\omega_n^{k+n/2})^2) \\ =A_1(\omega_{n/2}^k)-\omega_n^kA_2(\omega_{n/2}^k) \]
注意到上下两个式子只有符号不同,也就是说,我们带入 \(\omega_n^k\) 求出 \(A_1(\omega_n^k)\) 和 \(A_2(\omega _n^k)\) 就能得到两个点值。
这样我们就可以递归求解了,时间复杂度 \(T(n)=2T(\frac n2)+O(n)=O(n\log n)\)。
接下来解决 IDFT
。
我们得到了点值 \(s[k]=\sum \limits _ {i=0} ^ {n-1} f[i]\times (\omega _ n ^ k) ^ i\)。要求的多项式 \(n\times p[k]=\sum \limits _ {i=0} ^{n-1} s[i]\times (\omega_n^{n-k})^i\)。证明带入即可。这样 IDFT
就和 DFT
一样了,唯一的区别就是带入了逆单位根。
这样我们就可以写出递归版的 FFT
了,但是常数极大,可能模板题都过不了。我们需要些非递归版的 FFT
。
但是在 FFT
的时候需要将位置分奇偶,这样和原来的位置就不一样了。我们可以先将原序列变换好,放在最终的位置上,再一步一步向上合并。
注意到每个数最后的位置为其二进制翻转后的得到的位置。
struct Complex {
double x,y;
inline Complex(double a=0,double b=0) {x=a,y=b;}
inline Complex operator + (const Complex a) const {return Complex(x+a.x,y+a.y);}
inline Complex operator - (const Complex a) const {return Complex(x-a.x,y-a.y);}
inline Complex operator * (const Complex a) const {return Complex(x*a.x-y*a.y,x*a.y+y*a.x);}
}b[N],a[N];
int rev[N],f[N],g[N];
inline void FFT(Complex *a,int len,int f) {
int k=0;
while ((1<<k)<len) ++k;
for (int i=0;i<len;i++) {
rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
if (i<rev[i]) swap(a[i],a[rev[i]]);
}
for (int l=2;l<=len;l<<=1) {
Complex wn(cos(2*pi*f/l),sin(2*pi*f/l));
for (int j=0;j<len;j+=l) {
Complex w(1.0,0.0);
for (int k=j;k<j+l/2;k++) {
Complex t=w*a[k+l/2];
a[k+l/2]=a[k]-t,a[k]=a[k]+t,w=w*wn;
}
}
}
}
inline void DFT(Complex *a,int len) {
FFT(a,len,1);
}
inline void IDFT(Complex *a,int len) {
FFT(a,len,-1);
for (int i=0;i<=len;i++) a[i].x/=(double)len;
}
inline void Mul(Complex *a,Complex *b,int n,int m) {
int len=1; while (len<=(n+m)) len<<=1;
for (int i=n;i<len;i++) a[i]=Complex(0,0);
for (int i=m;i<len;i++) b[i]=Complex(0,0);
DFT(a,len),DFT(b,len);
for (int i=0;i<=len;i++) a[i]=a[i]*b[i];
IDFT(a,len);
}
快速数论变换 NTT
FFT
要算一大堆三角函数常数巨大,还可能会出现精度问题导致 WA
。这个时候需要 NTT
。
题目大意:给定两个多项式 \(f(x),g(x)\),求 \(f(x) * g(x)\)。即求:\(h[i]=\sum f[j]\times g[i-j]\)。对 \(998244353\) 取模。
这道模板题中最终系数保证不会大于 \(998244353\),所以可以视为对 \(998244353\) 取模。
前置知识:原根。
FFT
之所以快是因为 \(\omega\) 有优秀的性质,实际上原根也有同样的性质。相当于将 \(\omega_n^1\) 变成 \(g^{\frac {p-1}n}\)。其它部分都差不多。
但是需要注意:设模数为 \(p=r\times 2^k+1\)。NTT
只能处理 \(n\leq 2^k\) 的数据,比如模数为 \(998244353=119\times 2^{23}+1\),此时原根为 \(3\)。但当模数比较毒瘤,比如 \(10^9+7\),此时 \(k=1\) ,NTT
不能再使用。具体可以参照这篇博客。
int rev[N],a[N],b[N],c[N];
inline void NTT(int *a,int len,int f) {
int k=0;
while ((1<<k)<len) ++k;
for (int i=0;i<len;i++) {
rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
if (i<rev[i]) swap(a[i],a[rev[i]]);
}
for (int l=2;l<=len;l<<=1) {
int m=l>>1,w=Pow(G,(mod-1)/l);
if (f) w=Pow(w,mod-2);
for (int i=0;i<len;i+=l) {
int wk=1;
for (int j=0;j<m;j++,wk=1LL*wk*w%mod) {
int p=a[i+j],q=1LL*wk*a[i+j+m]%mod;
a[i+j]=(p+q)%mod,a[i+j+m]=(p-q+mod)%mod;
}
}
}
}
inline void DFT(int *a,int len) {
NTT(a,len,false);
}
inline void IDFT(int *a,int len) {
NTT(a,len,true); int invn=Pow(len,mod-2);
for (int i=0;i<len;i++) a[i]=1LL*a[i]*invn%mod;
}
多项式求逆
题目大意,给定多项式 \(f(x)\) ,求一个多项式 \(g(x)\) ,满足 \(f(x)*g(x)\equiv 1\pmod {x^n}\)。
主要的思想是倍增。假设我们已知 \(h(x)*f(x)\equiv 1 \pmod {x^{n/2}}\)。
有 \(g(x)-h(x)\equiv 1 \pmod {x^{n/2}}\)。
两边平方,有: \(g^2(x)-2g(x)h(x)+h^2(x)\equiv 1 \pmod {x^{n}}\)。
两边同乘 \(f(x)\),有:\(g(x)-2h(x)+h^2(x)f(x)\equiv 1 \pmod {x^n}\)。
所以有:\(g(x)=2h(x)-h^2(x)f(x)+1\pmod {x^n}\)。
倍增求解 \(g(x)\) 即可。
inline void PolyInv(int n,int *a,int *b) {
if (n==1) return b[0]=Pow(a[0],mod-2),void();
PolyInv((n+1)>>1,a,b); int len=1;
while (len<(n<<1)) len<<=1;
for (int i=0;i<n;i++) c[i]=a[i];
for (int i=n;i<len;i++) c[i]=0;
DFT(c,len),DFT(b,len);
for (int i=0;i<len;i++)
b[i]=1LL*(2-1LL*c[i]*b[i]%mod+mod)%mod*b[i]%mod;
IDFT(b,len);
for (int i=n;i<len;i++) b[i]=0;
}
多项式除法&取模
多项式对数函数
多项式指数函数
多项式快速幂
多项式多点求值
多项式快速插值
一些卡常技巧和提示
多项式题一般比较难写,而且常数大。下面给出一些卡常技巧和提示。
-
在
Mul(a,b)
的时候最好将a
和b
先放到临时变量中。 -
注意变量不要重名,
static
可以很好的解决这个问题,并优化内存,但是可能会增大常数。 -
在做一些多项式操作前一定要想清楚边界范围,并记得清空数组。这点可能使用
vector
会比较方便。 -
在做一些多项式操作前一定要记得清空数组。
-
使用取模优化
-
加一些
inline
和register
。 -
static
好像会变慢? -
使用一些操作减少
FFT
次数,比如三次变两次的技巧。 -
预处理可能要用到的所有原根次幂,这样在求原根次幂的时候回访问一段连续内存加快速度。具体可以如下:
int GPow[2][20][N]; inline void InitG() { for (int p=1;p<=19;p++) { int buf1=Pow(G,(Mod-1)/(1<<p)); int buf0=Pow(invG,(Mod-1)/(1<<p)); GPow[1][p][0]=GPow[0][p][0]=1; for (int i=1;i<(1<<p);i++) GPow[1][p][i]=1LL*GPow[1][p][i-1]*buf1%Mod, GPow[0][p][i]=1LL*GPow[0][p][i-1]*buf0%Mod; } }
同时将
NTT
改成这样inline void NTT(int *a,int len,int f) { int k=0; while ((1<<k)<len) ++k; for (Re int i=0;i<len;++i) if (i<rev[i]) swap(a[i],a[rev[i]]); for (Re int l=2,cnt=1;l<=len;l<<=1,cnt++) { int m=l>>1; for (Re int i=0;i<len;i+=l) { int *buf=GPow[f^1][cnt]; for (Re int j=0;j<m;j++,buf++) { int p=a[i+j],q=1LL*(*buf)*a[i+j+m]%Mod; a[i+j]=chk(p+q),a[i+j+m]=chk(p-q+Mod); } } } }
结语
咕咕咕