[学习笔记] 多项式全家桶系列
题外话：
本来这段时间都在洛谷写博客的...
但是因为洛谷猝不及防的改版了导致洛谷博客的搜索功能一度变的很垃圾..
所以又回来了...
等什么时候变好了再回去吧...
哦对了我去洛谷写博客的唯一理由是自己高兴。嗯，真的。
多项式全家桶系列
本来觉得这些没啥用，但是想学新东西了...那就随便学学叭
预处理单位根
upd：学了一下预处理单位根（实际上就是又学了一遍单位根到底是什么）
贴上\(\text{NTT}\)预处理单位根的部分代码：
首先是预处理部分：
w[0][0]=w[1][0]=1;
w[0][1]=ksm(3,(mod-1)/maxn);
w[1][1]=ksm((mod+1)/3,(mod-1)/maxn);
for(int i=2;i<maxn;i++)
w[0][i]=1ll*w[0][i-1]*w[0][1]%mod,w[1][i]=1ll*w[1][i-1]*w[1][1]%mod;
这个很显然
然后就是\(\text{NTT}\)里面了，直接把整个函数弄上来吧
void ntt(int *f,int g){
for(int i=1;i<lim;i++) if(i<rev[i]) swap(f[i],f[rev[i]]);
for(int mid=1;mid<lim;mid<<=1){
for(int R=mid<<1,j=0;j<lim;j+=R){
for(int k=0;k<mid;k++){
int x=f[j+k],y=1ll*w[g][maxn/R*k]*f[j+k+mid]%mod;
f[j+k]=x+y>=mod?x+y-mod:x+y,f[j+k+mid]=x-y<0?mod+x-y:x-y;
}
}
} if(g)
for(int in=ksm(lim),i=0;i<lim;i++) f[i]=1ll*f[i]*in%mod;
}
一定一定要记得 \(lim\) 先除再乘！！！！
这个是\(\text{FFT}\)的：
for(int i=0;i<2;i++)
for(int j=0;j<maxn;j++)
w[i][j]=cp(cos(2*pi/maxn*j),(i?1:-1)*sin(2*pi/maxn*j));
void fft(cp *f,int opt){
for(int i=1;i<lim;i++) if(i<rev[i]) swap(f[i],f[rev[i]]);
for(int mid=1;mid<lim;mid<<=1){
for(int R=mid<<1,j=0;j<lim;j+=R){
for(int k=0;k<mid;k++){
cp x=f[j+k],y=f[j+k+mid]*w[opt][maxn/R*k];
f[j+k]=x+y;f[j+k+mid]=x-y;
}
}
}
}
多项式求逆
Description
给定多项式 \(F(x)\)，请求出一个多项式 \(G(x)\)，满足 \(F(x)G(x)\equiv 1\; \operatorname{mod} x^n\)。对 \(998244353\) 取模。
Sol
推导过程挺简单也挺显然的...
结论就是 $$B(x)\equiv 2B'(x)-A(x)B'^2(x); \operatorname{mod} x^n$$
其中 $$A(x)B'(x)\equiv 1; \operatorname{mod} x^{\frac n2}$$
一个多项式有没有逆取决于常数项有没有逆
Code
void solveinv(int *a,int *b,int len){
if(len==1) return b[0]=ksm(a[0]),void();
solveinv(a,b,len>>1); lim=len+len;
for(int i=1;i<lim;i++) rev[i]=(rev[i>>1]>>1)|(i&1?lim>>1:0);
for(int i=0;i<len;i++) tmpa[i]=a[i];
for(int i=len;i<lim;i++) tmpa[i]=0;
ntt(tmpa,3),ntt(b,3);
for(int i=0;i<lim;i++) b[i]=1ll*b[i]*(2ll-1ll*tmpa[i]*b[i]%mod+mod)%mod;
ntt(b,g); for(int i=len;i<lim;i++) b[i]=0;
}
多项式除法
Description
给出 \(n\) 次多项式 \(F(x)\) 和 \(m\) 次多项式 \(G(x)\)，请求出多项式 \(Q(x),R(x)\)，满足:
- \(deg(Q)=n-m,deg(R)<m\)
- \(F(x)=G(x)Q(x)+R(x)\)
对 \(998244353\) 取模。
Sol
首先，对于任意多项式
显然有
也就是说通过这个可以把多项式的系数翻转。
对 \(F(x)=G(x)Q(x)+R(x)\) 进行这个操作：
注意到 \(x^{n-m+1}R^{rev}(x)\) 的最低项次数至少为 \(n-m+1\)，所以我们对式子取模：
因为 \(deg(Q)=n-m\)，所以这样是不会出问题的。
所以就可以算出 \(Q^{rev}(x)=F^{rev}(x)(G^{rev}(x))^{-1}\)
多项式求逆即可。
整个做法的关键就是消去余数。
Code
void solve(){
n=getint(),m=getint();
for(int i=0;i<=n;i++) f[i]=getint();
for(int i=0;i<=m;i++) g[i]=getint();
std::reverse(f,f+n+1),std::reverse(g,g+m+1); //翻转系数
for(int i=0;i<=n-m;i++) ing[i]=g[i];
//因为对 x^{n-m+1} 取模，所以后边的系数都应该为0
getni(ing,q,n-m+1); //求出g的逆q
for(int i=0;i<=n-m;i++) c[i]=f[i]; //后边的系数都应该为0
mul(q,c,n-m+1);
for(int i=0;i<lim;i++) c[i]=0;
std::reverse(q,q+n-m+1);
for(int i=0;i<=n-m;i++) printf("%d ",q[i]);puts("");
for(int i=n-m+1;i<lim;i++) q[i]=0;
std::reverse(f,f+n+1),std::reverse(g,g+m+1);
mul(q,g,n);
for(int i=0;i<m;i++) printf("%d ",(f[i]-q[i]+mod)%mod);
}
多项式开根
Description
给定多项式 \(A(x)\)，若存在一个多项式 \(B(x)\),满足：
- \(deg(B)\leq deg(A)\)
- \(B^2(x)\equiv A(x)\;\operatorname{mod} x^n\)
则称 \(B(x)\) 为 \(A(x)\) 在模 \(x^n\) 意义下的开根。
Sol
懒得写了...
感谢AlphaINF
补充一下，上边那个式子还可以再化简一步变成：
Code
void solvesqr(int len,int *a,int *b){
if(len==1) return b[0]=1,void();
solvesqr(len>>1,a,b);
solveinv(len,b,tmpb);
get(len);
for(int i=0;i<len;i++) tmpa[i]=a[i];
ntt(tmpb,1),ntt(tmpa,1);
for(int i=0;i<lim;i++) tmpa[i]=1ll*tmpa[i]*tmpb[i]%mod;
ntt(tmpa,-1);
for(int i=0,inv2=mod+1>>1;i<lim;i++) b[i]=1ll*(tmpa[i]+b[i])%mod*inv2%mod;
for(int i=len;i<lim;i++) b[i]=0;
for(int i=0;i<lim;i++) tmpa[i]=tmpb[i]=0;
}
三模数NTT
Description
求 \(F(x)G(x)\)，对 \(p\) 取模。\(p\) 不保证有原根。
Sol
拆成对三个模数 \(469762049,998244353,1004535809\) 分别做 \(NTT\)，它们的原根都是 \(3\)。
先合并前两个式子，得到:
设:
求出 \(x\) 在 \(\operatorname{mod}m_3\) 意义下的值：
在 \(\operatorname{mod} m_3\) 意义下， \(ym_3\) 被消掉了。所以有：
算出右边的值为 \(q\),令 \(x=km_3+q\)，代入\(ans=xM+C\):
因为 \(ans\in[0,m_1m_2m_3)\)，所以 \(k=0\)！ 所以 \(ans=qM+C\)！
Code
主要是合并部分
int get(int x,int y,int z,int mod){
ll M=p[1]*p[0];
ll tmp1=mul(1ll*x*p[1]%M,inv(p[1]%p[0],p[0]),M); //快速幂之前先取模
ll tmp2=mul(1ll*y*p[0]%M,inv(p[0]%p[1],p[1]),M);
ll tmp=(tmp1+tmp2)%M;
ll xx=1ll*(z-tmp%p[2]+p[2])%p[2]*inv(M%p[2],p[2])%p[2];
return ((xx%mod)*(M%mod)%mod+tmp%mod)%mod; //注意这里先对mod取模再乘起来
}
多项式exp，ln等学完高数马上填坑
拆系数FFT
Step1 基本介绍
设 \(F_P\) 表示多项式 \(P(x)\) DFT之后的结果。
对两个长度为 \(L\) 的多项式 \(A(x),B(x)\) 进行DFT，可以利用
对 \(P(x)\) 进行DFT，得到 \(F_p\)。
那么 \(F_q[k]=!(F_p[2L-k])\) ，（!表示取共轭）（证明见论文）
则
这就是两两合并计算DFT的方法，2次DFT优化成了1次。
主要部分的代码：
for(int i=0;i<lim;i++){
int to=i?lim-i:0;// emm这里并不知道为什么要这样写
// 嗯现在知道了 这个本质上就是 (lim-i)%lim <--- 这个是对的
b[i]=(a[i]*a[i]-!a[to]*!a[to])*cp(0,-0.25);
// 像下面这样写也可以
// tmpx[i]=(a[i]+!a[to])*cp(0.5,0);
// tmpy[i]=(a[i]-!a[to])*cp(0,-0.5);
// b[i]=tmpx[i]*tmpy[i];
}
这样搞一下我的FFT跑的就比NTT还快了！
Step2 拆系数FFT
设 \(M_0=\sqrt M\)，把多项式 \(A(x)\) 和 \(B(x)\) 拆成 \(k\times M_0+b\) 的形式，其中 \(k,b\) 都小于 \(M_0\)
所以两个多项式相乘就变成了 \(K_a(x)K_b(x)M_0^2+(K_a(x)B_b(x)+K_b(x)B_a(x))M_0+B_a(x)B_b(x)\)
然后可以把DFT从四次优化到两次，IDFT从三次优化到两次
这样只需要做四次DFT就行了！
说一下IDFT的两两合并：
假如说现在知道了多项式 \(A(x),B(x),C(x),D(x)\) 的DFT，要计算的是 \(A(x)B(x)\) 和 \(C(x)D(x)\)。令 \(DFT(A)\) 表示多项式 \(A(x)\) 的点值表达。
那么设：
我们需要对 \(dfta(x)\) 和 \(dftb(x)\) 进行IDFT。注意到这里的IDFT一定是实数，然后经过一番复杂的推导，再令：
那么 \(IDFT(p)\) 的实部除以 \(lim\) 就是 \(A(x)B(x)\)，虚部除以 \(lim\) 就是 \(C(x)D(x)\)。
代码：
void mul(int *a,int *b,int *ans,int n){
static cp p[N],q[N],c[N],d[N];
for(int i=0;i<lim;i++)
p[i]=cp(a[i]>>15,a[i]&32767),q[i]=cp(b[i]>>15,b[i]&32767);
fft(p,1),fft(q,1);
for(int i=0;i<lim;i++){
static cp aaa=cp(0.5,0),bbb=cp(0,-0.5),o=cp(0,1);
static cp ka,kb,ba,bb;
int j=(lim-i)%lim;
ka=(p[i]+!p[j])*aaa,kb=(p[i]-!p[j])*bbb;
ba=(q[i]+!q[j])*aaa;bb=(q[i]-!q[j])*bbb;
c[i]=ka*ba+kb*ba*o; d[i]=ka*bb+kb*bb*o;
} fft(c,-1),fft(d,-1);
for(int i=0;i<=n;i++){
ll A=((ll)(c[i].x+0.5))%mod,B=((ll)(c[i].y+0.5))%mod,C=((ll)(d[i].x+0.5))%mod,D=((ll)(d[i].y+0.5))%mod;
ans[i]=((A<<30)%mod+(B+C<<15)%mod+D)%mod;
}
}
多项式对数函数
Description
给定 \(n-1\) 次多项式 \(A(x)\)，求一个 \(\bmod x^n\) 下的多项式 \(B(x)\)，满足 \(B(x)\equiv \ln A(x)\;\bmod x^n\)。对 \(998244353\) 取模。
Sol
求ln还是挺傻逼的吧...
根据微积分基本知识，复合函数 \(y=f(g(x))\) 的导数和函数 \(y=f(u),u=g(x)\) 的导数间的关系为:
然后这题就是求 \(y=\ln(A(x))\)
根据上边的定义转变成 \(y'=\frac{A'(x)}{A(x)}\)
求导求逆再积分就行了
然后发现自己求逆的一点小bug就是设置长度设置为不小于 \(n\) 的 \(2\) 的整次幂就行了，不用大于 \(2n\)。
void ds(int *a,int *b,int n){
for(int i=1;i<n;i++) b[i-1]=1ll*a[i]*i%mod;b[n-1]=0;
}
void jf(int *a,int *b,int n){
for(int i=1;i<n;i++) b[i]=1ll*a[i-1]*inv(i)%mod;b[0]=0;
}
void solveln(int *a,int *b,int len){
lim=1;while(lim<len) lim<<=1;
solveinv(a,fa,lim);
lim=1;while(lim<len+len) lim<<=1;
for(int i=1;i<lim;i++) rev[i]=(rev[i>>1]>>1)|(i&1?lim>>1:0);
ds(a,da,len);
ntt(fa,3),ntt(da,3);
for(int i=0;i<lim;i++) fa[i]=1ll*fa[i]*da[i]%mod;
ntt(fa,g); jf(fa,b,len);
for(int i=0;i<lim;i++) fa[i]=da[i]=0;
}
多项式指数函数 (多项式exp)
Description
给定 \(n-1\) 次多项式 \(A(x)\)，求一个 \(\bmod x^n\) 下的多项式 \(B(x)\)，满足 \(B(x)\equiv e^{A(x)}\)。对 \(998244353\) 取模。
Sol
首先要知道牛顿迭代。 在针对一个函数的时候，牛顿迭代可以求函数零点，而针对一个多项式的时候求的就是多项式的零点。即对于一个多项式 \(G(x)\)，求满足条件的 \(G(F(x))\equiv 0\;(\bmod x^n)\) 的多项式 \(F(x)\)。
然后直接放公式就是 \(F(x)\equiv F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}\;(\bmod x^n)\)
其中 \(F_0(x)\) 满足 \(G(F_0(x))\equiv 0\;(\bmod x^{\frac n2})\)
回到多项式exp，那我们现在就是要计算 \(F(x)=e^{A(x)}\)
变形一下得
我们设 \(G(F(x))=\ln F(x)-A(x)\)，那么就是要求这个多项式的零点。如果我们把 \(F(x)\) 看做自变量，那 \(A(x)\) 就是常数，对 \(G(F(x))\) 进行求导，则 \(G'(F(x))=\frac1{F(x)}\)。
代入牛顿迭代的公式得：
Code
void solveexp(int *a,int *b,int len){
if(len==1) return b[0]=1,void();
solveexp(a,b,len>>1); solveln(b,tmpb,len); lim=len+len;
for(int i=1;i<lim;i++) rev[i]=(rev[i>>1]>>1)|(i&1?lim>>1:0);
tmpb[0]=(a[0]+1-tmpb[0]+mod)%mod;
for(int i=1;i<len;i++) tmpb[i]=(a[i]-tmpb[i]+mod)%mod;
ntt(b,3),ntt(tmpb,3);
for(int i=0;i<lim;i++) b[i]=1ll*b[i]*tmpb[i]%mod;
ntt(b,g);
for(int i=len;i<lim;i++) b[i]=tmpb[i]=0;
for(int i=0;i<len;i++) tmpb[i]=0;
}
