先粘个 \(\rm NTT\) 和 \(\rm FFT\) 的 板子。
inline void times(LL *f,LL *g,int n,int lim){
int kn=initr(n); NTT(f,kn,1); NTT(g,kn,1);
for(int i=0;i<kn;i++) f[i]=f[i]*g[i]%Mod;
NTT(f,kn,-1); NTT(g,kn,-1);
clr(f,lim,kn);
}
分治 \(NTT/FFT\)
大概就是维护下面这种较简单半在线卷积。
给定序列 \(g_{1\dots n - 1}\)，求序列 \(f_{0\dots n - 1}\)。
其中 \(f_i=\sum_{j=1}^if_{i-j}g_j\)，边界为 \(f_0=1\)。
运用 \(\rm CDQ\) 分治的思想，将区间 \([ L,R )\) 的求解 转化为 \([L,mid)\) 和 \([mid,R)\) 的求解。
先求解 \([L,mid)\) 相当于求出了左半部分的 \(f\) ,然后做一次 \(NTT\) 将左边的贡献加到右边，不断分治求解。
模板代码：
void DivNTT(int L,int R){
if(L==R-1) return ;
int mid=(L+R)/2,p=R-L; DivNTT(L,mid);
cpy(f2,f+L,p/2); cpy(g2,g,p); clr(f2,p/2,2*p); clr(g2,p,2*p);
times(f2,g2,2*p);
for(int i=mid;i<R;i++) (f[i]+=f2[i-L])%=Mod;
DivNTT(mid,R);
}
具体分治过程应视题目而定。
更加真实的半在线卷积：
已知 \(F[0]=G[0]=1,F[1]=1\)。
令 \(G[n]\) 表示 \(F_{1..n}\) 的前缀和。
\(F_i=\sum_{j=1}^iF_{i-j}G_j\)
这样按照我们原先的分治方法是不行的。
正确策略：
-
\(L=0\) 时
\(f[0,mid] * G[0,mid] \to f[mid,r)\)
-
\(L\not=0\) 时
\(\left.\begin{matrix} G[0,mid]*f[0,r-L)\\ F[L,mid]*G[0,r-L) \end{matrix}\right\}\to f[mid,R)\)
多项式求逆
给定多项式 \(F(x)\) ，求多项式 \(G(x)\) 使得 \(F(x) * G(x)\equiv 1 (\bmod \space x^n)\)。
求解：
采用倍增的过程。
我们已经求出了 \(F(x)* P(x) \equiv 1(\bmod \space x^{\frac{n}{2}})\)
因为 \(F(x) * G(x)\equiv 1 (\bmod \space x^n)\)。
所以 \(F(x) * G(x)\equiv 1 (\bmod \space x^{\frac{n}{2}})\)。
\(F(x) * (G(x)-P(x))\equiv 0 (\bmod \space x^{\frac{n}{2}})\)。
\(G(x)-P(x) \equiv 0 (\bmod \space x^{\frac{n}{2}})\)
平方得 \(G(x)^2+P(x)^2-2G(x)P(x)\equiv 0(\bmod \space x^n)\)
移项，同乘 \(F(x)\) ,\(G(x)\equiv 2P(x)-F(x)P(x)^2(\bmod \space x^n)\)
不断倍增求出 \(G\) 。
void invp(LL *f,int n){
static LL g[M+5],h[M+5],ff[M+5];int p=2;
g[0]=Pow(f[0],Mod-2);
for(;p<(n<<1);p<<=1){
int pn=p<<1; cpy(ff,f,p); cpy(h,g,p);
initr(pn); NTT(ff,pn,1); NTT(h,pn,1);
for(int i=0;i<pn;i++) g[i]=((2ll-ff[i]*h[i]%Mod)*h[i]%Mod+Mod)%Mod;
NTT(g,pn,-1); clr(g,p,pn);
}
cpy(f,g,n); clr(g,0,p); clr(h,0,p); clr(ff,0,p);
}
多项式开根
同样地倍增。
给定多项式 \(F(x)\) ，求多项式 \(G(x)\) 使得 $ G(x)^2\equiv F(x) (\bmod \space x^n)$。
多项式 \(F\) 的只有常数项时，用二次剩余解出 \(G\) 的常数项。蒟蒻不会。
在模板题中，常数项为 \(1\)，直接令 \(G\) 的常数项为 \(1\).
已知 \(P(x)^2\equiv A(x) (\bmod \space x^{\frac{n}{2}})\)
\(G(x)\equiv P(x) (\bmod \space x^{\frac{n}{2}})\)
\(G(x)-P(x)\equiv 0 (\bmod \space x^{\frac{n}{2}})\)
平方，\(G(x)^2-2G(x)P(x)+P(x)^2\equiv 0 (\bmod x^{n})\)
\(G(x)\equiv \frac{G(x)^2+P(x)^2}{2P(x)}\equiv \frac{F(x)+P(x)^2}{2P(x)}\equiv \frac{P(x)+\frac{F(x)}{P(x)}}{2}(\bmod \space x^n)\)
void Sqrtp(LL *f,int n){
static LL ff[M+5],gg[M+5],g[M+5]; memset(g,0,sizeof(g));
g[0]=1; int p=2;
for(;p<(n<<1);p<<=1) {
int pn=p<<1; cpy(ff,f,p); cpy(gg,g,p);
invp(gg,p); initr(pn);
NTT(ff,pn,1); NTT(gg,pn,1);
for(int i=0;i<pn;i++) ff[i]=ff[i]*gg[i]%Mod;
NTT(ff,pn,-1);
for(int i=0;i<p;i++) g[i]=(g[i]+ff[i])%Mod*inv2%Mod;
clr(g,p,pn);
}
cpy(f,g,n); clr(g,0,p); clr(ff,0,p); clr(gg,0,p);
}
多项式除法
给定 \(n\) 次多项式 \(F(x)\) 和 \(m\) 次多项式 \(G(x)\) 。求多项式 \(Q(x)\) 和 \(R(x)\)。
满足 \(F(x)=G(x)* Q(x)+R(x)\) 。
其中 \(Q(x)\) 的次数为 \(n-m\)。 \(R(x)\) 的次数为 \(m-1\)
定义 \(n\) 次多项式反转操作 \(F^T(x)=x^nF(x^{-1})\)。
例子：如 \(F(x)=3x^2+2x+1,F^T(x)=x^3F(x^{-1})=x^3+2x^2+3\)。
因为 \(F(x)=G(x)* Q(x)+R(x)\) ，
所以 \(F(x^{-1})=G(x^{-1})* Q(x^{-1})+R(x^{-1})\)
同时乘 \(x^n\)，\(x^nF(x^{-1})=x^nQ(x^{-1})* G(x^{-1})+x^nR(x^{-1})\)
得 \(F^T(x)=Q^T(x)* G^T(x)+x^{n-m+1} R^T(x)\)。
得到 \(Q^T(x)=F^T(x)* G^T(x)^{-1} (\bmod \space x^{n-m+1})\)
求出 \(Q^T(x)\) 并反转得到 \(Q(x)\)。
\(R(x)=F(x)-Q(x)* G(x)\)。
可能是对的代码？？？
inline void rev(LL *f,int n){
for(int i=0,j=n-1;i<j;i++,j--) swap(f[i],f[j]);
}
inline void mof(LL *f,LL *g,int n,int m){
static LL q[M+5],r[M+5];
int L=n-m+1;
rev(g,m); cpy(q,g,L); rev(g,m); invp(q,L);
rev(f,n); cpy(r,f,L); rev(f,n); times(q,r,L+L-1,L);
rev(q,L); prta(g,m+L-1);
times(g,q,m+L-1,m-1);
for(int i=0;i<m-1;i++) g[i]=(f[i]-g[i]+Mod)%Mod;
cpy(f,q,L); clr(f,L,n); clr(q,0,m); clr(r,0,L);
}
一些看不懂的东西：
简单微积分：
定义多项式 \(F(x)=\sum_{i=0}^na_ix^i\)。
-
定义复合多项式 \(F(G(x))=\sum\limits_{i=0}^na_iG(x)^i\)。
-
定义 \(F(x)\) 的求导: $F'(x)=\sum\limits_{i=0}^{n-1}(i+1)a_{i+1} x^i $。
-
定义 \(F(x)\) 的积分: \(\int(F(x))=\sum\limits_{i=1}^n\frac{a_{i-1}x^i}{i}\)。
-
复合多项式的求导：根据链式法则得，若\(A(x)=F(G(x))\)，同时对两边得主元 \(x\) 求导，则有 \(A'(x)=F'(G(x))G'(x)\)
inline void Getdao(LL *f,int n){
for(int i=1;i<n;i++) f[i-1]=f[i]*i%Mod;
f[n-1]=0;
}
inline void Getjf(LL *f,int n){
for(int i=n-1;i>=0;i--) f[i]=f[i-1]*inv[i]%Mod;
f[0]=0;
}
看得出来，积分和求导其实是互逆的过程。
多项式牛顿迭代 。
已知函数 \(G\) 且 \(G(F(x))=0\)，求多项式 \(F(\bmod \space x^n)\)。
定义 $F_* $ 表示 \(F(mod \space x^\frac{n}{2})\)。
结论：\(F(x)=F_* (x)-\frac{G(F_* (x))}{G'(F_* (x))}\)。
如此倍增的复杂度为 \(O(n\log n+O(n/2))=O(n\log n)\)。
多项式对数函数和指数函数（多项式 ln）
若 \(e^x=a\)，则 \(\ln(a)=x\)。
若 \(e^x=a\), 则 \(\exp(x)=a\)
由麦克劳林级数定义：
\[\ln(x)=-\sum_{i=1}\frac{(1-x)^i}{i} \]\[\exp(x)=\sum_{i=0}\frac{x^i}{i!} \]一些性质：\(\exp(ln(x))=x,\exp(ln(F(x))+\ln(G(x)))=F(x)* G(x)\)。
给定多项式 \(F(x)\) 求 \(\ln(F(x))\) 和 \(\exp(F(x))\)。
ln
给定多项式 \(F\) 求多项式 \(G(x)\equiv \ln (F(x))(\bmod \space x^n)\)。保证 \(a_0=1\)
\(\ln'(x)=\frac{1}{x}\)
对 \(G(x)\equiv \ln(F(x))\) 两边同时求导。
\(G'(x)\equiv \ln'(F(x))F'(x)\equiv \frac{F'(x)}{F(x)}(\bmod \space x^n)\) ，再用积分求出 \(G=\int(G'(x))\)。
inline void lnp(LL *f,int n){
static LL g[M+5];
cpy(g,f,n); invp(g,n);
Getdao(f,n); times(f,g,n+n-1,n);
Getjf(f,n); clr(g,0,n);
}
exp
给定多项式 \(F\) 求多项式 \(G(x)\equiv \exp (F(x))(\bmod \space x^n)\)。保证 \(a_0=0\)
利用牛顿迭代求 \(G\)，这里用 \(B\) 表示。 回忆式子 \(B(x)=B_* (x)-\frac{G(B_* (x))}{G'(B_* (x))}\)。
根据 \(e^{F(x)}= B(x)\) 得 \(\ln(B(x))-F(x) = 0\)。构造 $G(B_* (x))=\ln(B_* (x))-F(x) $ ，所以 \(G'(B_* (x))=\ln'(B_* (x))=B_* (x)^{-1}\)。
\(B(x)=B_* (x)-\frac{\ln(B_* (x))-F(x)}{B_x ^{-1}}=(1-\ln(B_* (x))+F(x)) * B_* (x)\)。
inline void exp(LL *f,int n){
static LL g[M+5],gg[M+5]; int p=2;
g[0]=1;
for(;p<(n<<1);p<<=1){
cpy(gg,g,p); lnp(gg,p);
for(int i=0;i<p;i++) gg[i]=(f[i]-gg[i]+Mod)%Mod;
gg[0]=(gg[0]+1)%Mod; times(g,gg,p*2,p);
}
cpy(f,g,p); clr(g,0,p); clr(gg,0,p);
}
多项式快速幂
给定 \(n-1\) 次多项式 \(F(x)\) 求\(F(x)^k (\bmod \space x^n)\)。
暴力 \(\rm NTT\)
对于常见的题目中，一般是加速 \(dp\) 的求解。我们可以采用 \(O(n\log n\log k)\) 。即写成一般快速幂的形式,用 \(times\) 函数来让多项式相乘。
这种做法还有以下优点：
见例题： 有 \(n\) 个数 \(a_{1...n}\)，构造一个序列 \(b_{1..T} (b_i\in a)\) ，求 \(S_b\equiv p\bmod m\) 。我们有一个矩阵加速做法，但是为复杂度 \(p^3\log T\)。
但是我们思考以下，这个矩阵相乘的形式为： \(F'[k]=\sum\limits_{i+j\equiv k}F[i]\times G[j]\)。
这是一个加法卷积的形式，我们可以做一遍加法卷积之后再让 \(F'[k]+=F'[k+m]\)。
然后用暴力的倍增求快速幂即可。但是 \(O(n\log n)\) 的做法只能维护普通的 \(F(x)^k (\bmod \space x^n)\) 的形式。
利用 \(\ln\) 和 \(\exp\)。
模板题弱化版：保证常数项为 \(1\)。
\(F^k(x)\equiv e^{\ln (F(x))\times k} (\bmod \space x^n)\)
见 \(Ei\) 大神的博客可以知道 \(k\) 可以直接对模数取模。
inline void polypow(LL *f,int n,int k){
lnp(f,n);
for(int i=0;i<n;i++) f[i]=f[i]*k%Mod;
exp(f,n);
}
强化版 ：
不保证常数项为 \(1\)。
我们可以取一个常数 \(c\) 。对用 \(c\) 除多项式 \(F\)，使得常数项为 \(1\)。
那么 \(F^k(x)=F_2^k(x)\times c^k\)。
但是这里，\(f_0=0\) 会出问题。
所以我们找到次数最小的非零项 设为 \(c\times x^{px}\)
左边这个多项式的快速幂可以用上面的代码求出。
在这个结果前面接 \(px\times k\) 个 \(0\) 即可。
注意，这里的 \(c^k\) 根据费马小定理，\(k\) 对模数减一取模。
inline void polypow(LL *f,int n,int k1,int k2){
int px=0; while(f[px]==0 && px<n) px++;
if(px==n) return ;
if(1ll*k2*px>=n) {clr(f,0,n); return ;}
LL iv1=Pow(f[px],Mod-2),x=Pow(f[px],k2);
for(int i=px;i<n;i++) f[i]=f[i]*iv1%Mod;
lnp(f+px,n-px);
for(int i=px;i<n;i++) f[i]=f[i]*k1%Mod;
exp(f+px,n-px);
for(int i=px;i<n;i++) f[i]=f[i]*x%Mod;
static LL g[M+5];
for(int i=k2*px;i<n;i++) g[i]=f[px+i-k2*px];
cpy(f,g,n); clr(g,0,n);
}