【瞎口胡】多项式操作
前置
快速傅里叶变换 FFT
多项式的基石操作。
快速沃尔什变换 FWT
位运算卷积。鸽了。
快速数论变换 NTT
把 FFT 搬到了模意义下，终于可以做计数问题啦。
多项式牛顿迭代
简单粗暴的推导方式。
基本操作
封装
为了学习多项式的时候更加顺手，封装板子是很有必要的，而且也方便贺。试想 CF 最后十分钟的时候切了个多项式题，然后发现把求逆和 exp 的板子分别复制过来合不到一起，那感觉心态直接爆了。
这个封装的原版是 zhouyuyang 神仙的，可以去大佬的博客找到。
我们考虑把多项式封装到 namespace 里面，这样一定程度上避免了和程序变量重名，如果没有重名的风险也可以直接 using namespace。
namespace PolyNomial{
# define ull unsigned long long
# define ll long long
# define mo MOD
# define poly std::vector <int>
const int FFTN=1<<18,MOD=998244353; // 做 NTT 的最大长度，可以调整成更大的或者更小的，要保证是 2 的幂次
int w[FFTN+10],W[FFTN+10],rev[FFTN+10]; // 小 w 和大 W 存单位根，rev 存蝴蝶变换的位置
inline void qadd(int &x,int v){ // 经典加法
x+=v,((x>=mo)&&(x-=mo));
return;
}
inline int qpow(int d,int p){ // 经典快速幂
int ans=1;
for(;p;p>>=1,d=1ll*d*d%mo) (p&1)&&(ans=1ll*ans*d%mo);
return ans;
}
inline void initW(void){ // 大 W 存的是 FFTN 次单位根
W[0]=1,W[1]=qpow(3,(mo-1)/FFTN);
for(int i=2;i<=FFTN;++i) W[i]=1ll*W[i-1]*W[1]%mo;
return;
}
inline int initRev(int len){ // 蝴蝶变换,len 是多项式**次数**（项数+1）
int n=1;
while(n<=len) n<<=1;
for(int i=0;i<n;++i) rev[i]=(rev[i>>1]>>1)|((i&1)?(n>>1):0);
return n;
}
DFT(NTT) 和 IDTT(INTT)
然后就是涉及到多项式的部分。我们希望多项式操作后的结果是一个 vector，这样操作起来方便一点，不用一会 new 个数组出来一会又 delete 回去，看着太蠢了。
但是 NTT 的时候如果直接传 vector 进去会有长度不够之类的问题，于是我们用声明好的数组 A,B 来做 DFT 和 IDFT，这样就没啥越界的问题了。
接下来考虑优化 NTT。我们注意到，取模巨大多慢，于是我们把点值数组 va[] 变成 unsigned long long 来减少取模次数，每 \(16\) 层取一次模。考虑到内层乘法只有点值乘单位根，后者是不超过 \(998244352\) 的，并且每一层下来 va[] 最多加上 \(998244352\)。考虑 \(998244352\times 998244352 \times 16\)，取 \(\log_2\) 约为 \(63.789635526\)，\(64\) 位的 ull 刚好能存下。
int A[FFTN+10],B[FFTN+10],C[FFTN+10]; // 最多的时候需要三个临时数组，一般只用 A,B
ull va[FFTN+10]; // NTT 的临时数组
inline void NTT(int *F,int len){
for(int i=0;i<len;++i) va[rev[i]]=F[i];
for(int l=1;l<len;l<<=1){ // merge 2 * l -> 1 * 2l
int siz=FFTN/(l<<1);
for(int i=0,j=0;i<l;++i,j+=siz) w[i]=W[j]; // omega_{n}^{i} = W_j
for(int s=0;s<len;s+=(l<<1))
for(int j=0;j<l;++j){
int v=va[s+j+l]*w[j]%mo;
va[s+j+l]=va[s+j]+mo-v,va[s+j]+=v;
}
if(l==(1<<15)) for(int i=0;i<len;++i) va[i]%=mo;
}
for(int i=0;i<len;++i) F[i]=va[i]%mo;
return;
}
inline void INTT(int *F,int len){
for(int i=0;i<len;++i) va[rev[i]]=F[i];
for(int l=1;l<len;l<<=1){ // 2 * l -> 1 * 2l
int siz=FFTN/(l<<1);
for(int i=0,j=FFTN;i<l;++i,j-=siz) w[i]=W[j]; // 要反着取单位根了
for(int s=0;s<len;s+=(l<<1))
for(int j=0;j<l;++j){
int v=va[s+j+l]*w[j]%mo;
va[s+j+l]=va[s+j]+mo-v,va[s+j]+=v;
}
if(l==(1<<15)) for(int i=0;i<len;++i) va[i]%=mo;
}
int invn=qpow(len,MOD-2);
for(int i=0;i<len;++i) F[i]=va[i]*invn%mo;
return;
}
调整多项式
多项式在操作之后次数可能会谜之减少，我们不希望多项式最高次项是 \(0\) 对吧？于是需要实现调整多项式的功能。
inline poly Sweep(poly a){ // 返回多项式
for(;a.size()>1&&!a.back();) a.pop_back();
return a;
}
inline void Sweeps(poly &a){ // 直接在传进来的参数上调整
for(;a.size()>1&&!a.back();) a.pop_back();
return;
}
// ps: 我们不希望空多项式存在（两个相同的多项式相减按照定义也会剩个 0），于是零次项的 0 我们不管
多项式运算
多项式加减法
计算两个多项式相加和相减的结果。
inline poly Plus(const poly &a,const poly &b){ // 常引用，省去了 copy 一份 a 和 b 的时间
// 理论上不写 const 也行，感觉写了清楚一点，一眼就能看出“这只是个参数”
int dea=a.size()-1,deb=b.size()-1;
assert(dea>=0&&deb>=0); // 说不定从哪里传了个空多项式进来，要警惕 不过不写这个也没事
poly ans(std::max(dea,deb)+1);
for(int i=0;i<=dea;++i) ans[i]=a[i];
for(int i=0;i<=deb;++i) qadd(ans[i],b[i]);
Sweeps(ans);
return ans;
}
inline poly Minus(const poly &a,const poly &b){
int dea=a.size()-1,deb=b.size()-1;
assert(dea>=0&&deb>=0);
poly ans(std::max(dea,deb)+1);
for(int i=0;i<=dea;++i) ans[i]=a[i];
for(int i=0;i<=deb;++i) qadd(ans[i],mo-b[i]);
Sweeps(ans);
return ans;
}
多项式乘法
计算两个多项式的积。
分别 NTT，然后乘起来 INTT 回去即可。
inline poly Mul(const poly &a,int k){ // 数乘
int dea=a.size();
poly ans(dea+1);
for(int i=0;i<=dea;++i) ans[i]=1ll*a[i]*k%mo;
return ans;
}
inline poly Mul(const poly &a,const poly &b){
int dea=a.size()-1,deb=b.size()-1;
assert(dea>=0&&deb>=0);
poly ans(dea+deb+1);
int blen=initRev(dea+deb); // 确定 ntt 的长度
for(int i=0;i<blen;++i) A[i]=((i<=dea)?a[i]:0),B[i]=((i<=deb)?b[i]:0); // 记得清零
NTT(A,blen),NTT(B,blen);
for(int i=0;i<blen;++i) A[i]=1ll*A[i]*B[i]%mo;
INTT(A,blen);
for(int i=0;i<=dea+deb;++i) ans[i]=A[i];
return ans;
}
多项式乘法逆
给定多项式 \(F(x)\)，求 \(F(x)\) 在模 \(x^n\) 意义下的乘法逆 \(G(x)\)，即：\(F(x)G(x) \equiv 1 \pmod{x^n}\)。
牛顿迭代得
递归，维护 \(G_0\) 和 \(F\) 做乘法即可。
inline void Getinv(int *F,int *G,int len){ // given F, calc F^{-1} = G
// G = inv(F) mod xn ; G_0 = inv(F) modxn/2 ; G = G_0(2-G_0*F)
if(len==1){ // 边界，此时
G[0]=qpow(F[0],MOD-2);
return;
}
int mlen=(len+1)>>1;
Getinv(F,G,mlen);
int blen=initRev(2*len);
for(int i=0;i<blen;++i) A[i]=((i<len)?F[i]:0),B[i]=((i<mlen)?G[i]:0);
NTT(A,blen); // A 对应 F, B 对应 G0
NTT(B,blen);
for(int i=0;i<blen;++i) A[i]=1ll*B[i]*(2+mo-1ll*B[i]*A[i]%mo)%mo;
INTT(A,blen);
for(int i=0;i<len;++i) G[i]=A[i];
return;
}
poly Getinv(const poly &a,int len){ // 外层调用只用这个函数
int *f=new int[len],*g=new int[len]; // 需要 2*len 的额外空间
for(int i=0;i<len;++i) f[i]=((i<(int)a.size())?a[i]:0),g[i]=0;
Getinv(f,g,len);
poly ans(len);
for(int i=0;i<len;++i) ans[i]=g[i];
delete[] f;
delete[] g;
return ans;
}
多项式除法
给定一个 \(n\) 次多项式 \(F(x)\) 和一个 \(m\) 次多项式 \(G(x)\) ，请求出多项式 \(Q(x)\)，\(R(x)\)，满足以下条件：
- \(Q(x)\) 次数为 \(n-m\)，\(R(x)\) 次数小于 \(m\)
- \(F(x) = Q(x)G(x) + R(x)\)
需要一点魔法。首先，\(n\) 次多项式 \(F(x)\) 把系数翻转形成的多项式是 \(x^nF(\dfrac{1}{x})\)，记作 \(F'(x)\)。可以这样想，原来次数越高的项，要乘的 \(\dfrac{1}{x}\) 就越多，在结果中的次数就越低。
实际上次数不足 \(n\) 次的多项式可以在高次项补上 \(0\)，然后这样翻转，效果也是一样的。
为了方便起见，我们下文中直接认为 \(R(x)\) 次数是 \(m-1\)。
运用魔法
于是先对 \(G'(x)\) 求逆，然后多项式乘法即可算出 \(Q'(x)\)，从而算出 \(Q(x)\)。
考虑 \(R(x)\)，它等于 \(F(x) - Q(x)R(x)\)，算出 \(Q(x)\) 后容易求得。
inline poly Rev(poly a){
std::reverse(a.begin(),a.end());
return a;
}
inline void Revs(poly &a){ // 注意这里只对多项式进行了翻转，没有去掉高位的 0，因为这个不影响操作
std::reverse(a.begin(),a.end());
return;
}
inline poly Divide(poly a,poly b){ // 计算 A / B
int dea=a.size()-1,deb=b.size()-1;
if(dea<deb) return poly(1); // 返回只有一个 0 的多项式
Revs(a),Revs(b),a.resize(dea-deb+1);
b=Mul(Getinv(b,dea-deb+1),a);
b.resize(dea-deb+1),Revs(b); // 次数是 dea-deb，于是 vector 里有 dea-deb+1 项
return b;
}
inline poly Modulo(poly a,poly b){ // 计算 A mod B
poly res=Minus(a,Mul(Divide(a,b),b));
res.resize(b.size()-1);
return res;
}
多项式开根
给定 \(F(x)\)，求 \(G(x)\) 使得 \(G^2(x) \equiv F(x) \pmod{x^n}\)。
牛顿迭代得
每层求逆即可。如果想常数小一点可以同时维护 \(G_0\) 和 \(G_0\) 的逆。另外，如果 \([x^0]F(x) \neq 1\)，要加上求二次剩余的板子。
inline void Getsqrt(int *F,int *G,int len){
if(len==1) return G[0]=1,void(); // 如果 F[0]!=1, 要用上二次剩余
int mlen=(len+1)>>1;
Getsqrt(F,G,mlen);
int *invc=new int[len];
for(int i=0;i<len;++i) invc[i]=0;
Getinv(G,invc,len); // 存 G0 的逆
int blen=initRev(2*len); // 要在模 x^{2n} 意义下做，因为 F(x) 和 1/G0(x) 都是 n-1 次多项式，相乘是 2n-2 次
for(int i=0;i<blen;++i) A[i]=((i<len)?F[i]:0),B[i]=((i<len)?invc[i]:0);
NTT(A,blen),NTT(B,blen);
for(int i=0;i<blen;++i) A[i]=1ll*A[i]*B[i]%mo;
INTT(A,blen); // 计算 F(x)/G0(x)
for(int i=0;i<len;++i) G[i]=1ll*(G[i]+A[i])*(mo+1)/2ll%mo;
delete[] invc;
return;
}
inline poly Getsqrt(const poly &a,int len){
int *f=new int[len],*g=new int[len];
for(int i=0;i<len;++i) f[i]=((i<(int)a.size())?a[i]:0),g[i]=0;
Getsqrt(f,g,len);
poly ans(len);
for(int i=0;i<len;++i) ans[i]=g[i];
delete[] f;
delete[] g;
return ans;
}
多项式对数函数（多项式 ln）
给定 \(F(x)\)，求 \(G(x)\) 使得 \(G(x) \equiv \ln F(x) \pmod{x^n}\)。
我们现在还不知道 \(\ln F(x)\) 到底是什么。所幸我们知道 \((\ln x)'=\dfrac{1}{x}\)，根据复合函数的求导法则，\((\ln F(x))' = \dfrac{1}{F(x)}\times F'(x) = \dfrac{F'(x)}{F(x)}\)。然后对它积分就可以算出 \(\ln F(x)+c\) 了，其中 \(c\) 是一个常数，因为求导再积分会丢失常数信息。
那 \(G(x)\) 的常数项是多少呢？我们知道 \(G(x)\) 是个多项式，带入 \(x=0\)，得 \(G(0)=\ln F(0) = \ln ([x^0]F(x))\)。
参考证明，容易发现 \(F(x)\) 常数项不为 \(1\) 时 \(\ln F(0)\) 是超越数，这也很符合直觉。模意义下超越数没有意义，于是我们有以下结论：当且仅当 \(F(x)\) 常数项等于 \(1\)，\(F\) 在模 \(x^n\) 意义下有对数多项式。
如果 \([x^0]F(x)=1\)，我们发现积分常数等于 \(0\)，这和求导再积分回去的常数项相同，于是我们不用刻意计算积分常数。
inline poly Int(const poly &a){ // 积分：integral(x^n) = 1/(n+1) x^{n+1}
int dea=a.size()-1;
poly ans(dea+2);
ans[0]=0;
for(int i=0;i<=dea;++i) ans[i+1]=1ll*qpow(i+1,mo-2)*a[i]%mo;
return ans;
}
inline poly Der(const poly &a){ // 求导: (x^n)' = n * x^{n-1}
int dea=a.size()-1;
poly ans(dea);
if(!dea) ans.push_back(0);
for(int i=1;i<=dea;++i) ans[i-1]=1ll*i*a[i]%mo;
return ans;
}
inline poly Getln(poly a,int len){ // 求 ln
assert(a[0]==1);
if((int)a.size()>len) a.resize(len);
poly res=Int(Mul(Der(a),Getinv(a,len)));
res.resize(len);
return res;
}
inline void Getln(int *F,int *G,int len){ // 求 F 的 ln, 存到 G 中
// 待会算 exp 的时候有用
poly a(len);
for(int i=0;i<len;++i) a[i]=F[i];
a=Getln(a,len);
for(int i=0;i<len;++i) G[i]=a[i];
return;
}
多项式指数函数（多项式 exp）
给定 \(F(x)\)，求 \(G(x)\) 使得 \(G(x) \equiv \exp F(x) \pmod{x^n}\)。
牛顿迭代得
和上面相似，求 exp 的时候要保证 \([x^0]F(x) = 0\)。
inline void Getexp(int *F,int *G,int len){
if(len==1) return G[0]=1,void();
int mlen=(len+1)>>1;
Getexp(F,G,mlen);
int *lnc=new int[len];
for(int i=0;i<len;++i) lnc[i]=0;
Getln(G,lnc,len);
int blen=initRev(2*len); // 同理，两个模 x^n 意义下的多项式相乘可能不止 n 次
for(int i=0;i<blen;++i) A[i]=((i<len)?F[i]:0),B[i]=((i<mlen)?G[i]:0),C[i]=((i<len)?lnc[i]:0); // A: F(x) B: G0(x) C: ln(G0(x))
NTT(A,blen),NTT(B,blen),NTT(C,blen);
for(int i=0;i<blen;++i) A[i]=1ll*B[i]*(1ll+mo-C[i]+A[i])%mo;
INTT(A,blen);
for(int i=0;i<len;++i) G[i]=A[i];
delete[] lnc;
return;
}
inline poly Getexp(const poly &a,int len){
int *f=new int[len],*g=new int[len];
for(int i=0;i<len;++i) f[i]=((i<(int)a.size())?a[i]:0),g[i]=0;
Getexp(f,g,len);
poly ans(len);
for(int i=0;i<len;++i) ans[i]=g[i];
delete[] f;
delete[] g;
return ans;
}
多项式快速幂
给定多项式 \(F(x)\) 和 \(k\)，求 \(G(x) \equiv F^k(x) \pmod{x^n}\)。保证 \([x^0]F(x) = 1\)。
可以类比矩阵快速幂的做法，但复杂度是 \(O(n \log^2 n)\) 的。
但这里保证 \([x^0]F(x)=1\)，而且是模 \(x^n\) 意义下，于是可以对 \(F\) 直接求 \(\ln\)，然后多项式乘上 \(k\)，再 \(\exp\) 回去。
inline poly __Getpower(poly a,int len,int k){
if((int)a.size()>len) a.resize(len);
a=Getln(a,len),a=Mul(a,k),a=Getexp(a,len);
return a;
}
设 \(F(x) = \sum \limits_{i=0}^{n-1} a_ix^i\)。对于 \(a_0\) 不为 \(1\) 的情况怎么办？如果 \(a_0\) 不是 \(0\)，那直接给每项除个 \(a_0\)，求完快速幂再给每项乘上 \(a_0^k\)。
如果 \(a_0\) 是 \(0\) 那咋办？我们可以提公因式，例如，\(F(x)=2x^3+6x^4+4x^5\)，那么就可以变成 \(F(x)=2x^3(1+3x+2x^2)\)。对后面括号里的部分快速幂，完了再乘上 \((2x^3)^k\)。当然不可能傻乎乎地去把结果和 \(2^kx^{3k}\) 做多项式乘法，直接把结果多项式往右移 \(3k\) 位，再给每个系数乘上 \(2^k\) 即可。
对于一般的情况类比即可。
inline poly Getpower(poly a,int len,int k,int pk){ // pk: k 对 phi(p) 取模的结果，方便算 a_0^k ps: 模板题的 k 大得离谱
if((int)a.size()>len) a.resize(len);
int t=0,b,dea=a.size()-1;
while(!a[t]&&t<=dea) ++t;
assert(t<=dea),b=a[t]; // 如果多项式全 0 那混进来了怪东西
poly res(len);
if(1ll*t*k>len) return res; // t*k >= len 那模 x^len 就只剩下 0 了
for(int i=0;i<=dea;++i) a[i]=((i+t<=dea)?(1ll*a[i+t]*qpow(b,mo-2)%MOD):0); // 提公因式
a=__Getpower(a,len,k); // 快速幂
int dk=qpow(b,pk),rlen=t*k; // dk: 提出来的系数
for(int i=0;i<=dea;++i) a[i]=1ll*a[i]*dk%mo,((i+rlen<=dea)&&(res[i+rlen]=a[i])); // 往右移多项式并乘上系数
return res;
}
另外，做模板题的时候有一个小细节：如果 \(k\) 非常大并且能提出来至少一个 \(x\)，那这个快速幂的结果一定是全 \(0\) 了。但 \(k\) 取完模可能变得很小，于是在读入的时候就要判断，如果在任意时刻 \(k \geq 10^7\) 且 \(a_0>0\)，那么答案就一定是全 \(0\)。
