多项式简陋入门
多项式全家桶
然而并没有多点求值,快速插值,转下降/上升幂,复合,复合逆
疯狂多项式,v我50
namespace efX_poly
{
const int maxlen=(1<<23)+1,maxSqrt=1e5+1;
inline int add(int x,int y,int m) { x+=y; return x>=m?x-m:x; }
inline int sub(int x,int y,int m) { x-=y; return x<0?x+m:x; }
inline int mul(int x,int y,int m) { ll z=1ll*x*y; return z>=m?z%m:z; }
inline void Add(int &x,int y,int m) { x=add(x,y,m); }
inline void Sub(int &x,int y,int m) { x=sub(x,y,m); }
inline void Mul(int &x,int y,int m) { x=mul(x,y,m); }
inline int qpow(int a,int x,int m)
{
int r=1;
while(x)
{
if(x&1) Mul(r,a,m);
Mul(a,a,m),x>>=1;
} return r;
} inline int _inv(int a,int m) { return qpow(a,m-2,m); }
namespace set_mod
{
static vector<int> d;
inline void init(int p)
{
int tmp=p-1;
if(!d.empty()) d.clear();
for(int i=2;1ll*i*i<=tmp;++i)
if(tmp%i==0)
{
d.emplace_back(i);
while(tmp%i==0) tmp/=i;
}
if(tmp>1) d.emplace_back(tmp);
}
inline bool check(int x,int p)
{
for(int v:d)
if(qpow(x,(p-1)/v,p)==1)
return false;
return true;
}
inline int getg(int p)
{
init(p);
irep(i,2,p)
if(check(i,p))
return i;
return -1;
}
}
int r[maxlen];
inline void make_rev(int lim,int l)
{
static int his=-1;
if(his==lim) return;
for(int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
his=lim;
}
inline pii getLimit(int n)
{
int lim=1,l=0;
while(lim<n) lim<<=1,++l;
return make_pair(lim,l);
}
int inv[maxlen];
int tmp[maxlen],tmpb[maxlen];
inline void view_arr(int *a,int lim) { cerr<<"## "; for(int i=0;i<lim;++i) {cerr<<a[i]<<' ';} cerr_out(); }
class poly
{
public:
const int p,g,gi;
private:
int G[2][maxSqrt+1],GI[2][maxSqrt+1];
inline void init()
{
memset(inv,0,sizeof inv);//can't use it any more
G[0][0]=GI[0][0]=1;
for(int i=1;i<=maxSqrt;++i)
{
G[0][i]=mul(G[0][i-1],g,p);
GI[0][i]=mul(GI[0][i-1],gi,p);
}
G[1][0]=GI[1][0]=1;
for(int i=1;i<=maxSqrt;++i)
{
G[1][i]=mul(G[1][i-1],G[0][maxSqrt],p);
GI[1][i]=mul(GI[1][i-1],GI[0][maxSqrt],p);
}
}
inline int very_quick_pow(bool o,int x)
{
if(o) return mul(G[1][x/maxSqrt],G[0][x%maxSqrt],p);
return mul(GI[1][x/maxSqrt],GI[0][x%maxSqrt],p);
}
public:
inline void NTT(int *a,int lim,int op)
{
for(int i=0;i<lim;++i) if(i<r[i]) swap(a[i],a[r[i]]);
for(int mid=1,wn;mid<lim;mid<<=1)
{
wn=very_quick_pow(op==1,(p-1)/(mid<<1));
for(int j=0,w;j<lim;j+=(mid<<1))
{
w=1;
for(int k=0,x,y;k<mid;++k,w=mul(w,wn,p))
x=a[j+k],
y=mul(w,a[j+k+mid],p),
a[j+k]=add(x,y,p),
a[j+k+mid]=sub(x,y,p);
}
}
if(op==-1)
{
int iv=_inv(lim,p);
for(int i=0;i<lim;++i) Mul(a[i],iv,p);
}
}
inline void clear(int *a,int st,int ed) { for(int i=st;i<ed;++i) a[i]=0; }
inline void arrSub(int *a,int *b,int lim) { for(int i=0;i<lim;++i) Sub(a[i],b[i],p); }
inline void arrSub(int *a,int c,int lim) { for(int i=0;i<lim;++i) Sub(a[i],c,p); }
inline void arrAdd(int *a,int *b,int lim) { for(int i=0;i<lim;++i) Add(a[i],b[i],p); }
inline void arrAdd(int *a,int c,int lim) { for(int i=0;i<lim;++i) Add(a[i],c,p); }
inline void Reverse(int *a,int *b,int n) { for(int i=0;i<n;++i) b[i]=a[n-i-1]; }
inline void arrMul(int *a,int *b,int lim) { for(int i=0;i<lim;++i) Mul(a[i],b[i],p); }
inline void arrMul(int *a,int c,int lim) { for(int i=0;i<lim;++i) Mul(a[i],c,p); }
inline void Inv(int *a,int *b,int n)
{
static int tmp[maxlen];//开成全局就gg了
if(n==1) { b[0]=_inv(a[0],p); return; }
Inv(a,b,(n+1)>>1);
int l=0,lim=1;
while(lim<(n<<1)) lim<<=1,++l;
make_rev(lim,l);
memcpy(tmp,a,sizeof(int)*n);
for(int i=n;i<lim;++i) tmp[i]=0;
NTT(tmp,lim,1),NTT(b,lim,1);
for(int i=0;i<lim;++i) b[i]=mul(sub(2,mul(b[i],tmp[i],p),p),b[i],p);
NTT(b,lim,-1);//tmp就不用管了,以后用的时候会刷新的
for(int i=n;i<lim;++i) b[i]=0;
}
inline void Sqrt(int *a,int *b,int n)
{
if(n==1) { b[0]=1; return; }
if(!inv[2]) inv[2]=_inv(2,p);
Sqrt(a,b,(n+1)>>1);
memset(tmp,0,sizeof(int)*(n<<1));
Inv(b,tmp,n);
int l=0,lim=1;
while(lim<(n<<1)) lim<<=1,++l;
make_rev(lim,l);
memcpy(tmpb,a,sizeof(int)*n);
clear(tmpb,n,lim);
NTT(tmpb,lim,1);
NTT(tmp,lim,1);
NTT(b,lim,1);
for(int i=0;i<lim;++i)
Add(b[i],mul(tmpb[i],tmp[i],p),p),
Mul(b[i],inv[2],p);
NTT(b,lim,-1);
clear(b,n,lim);
}//be very very careful!
inline void Int(int *a,int *b,int n)
{
for(int i=n-1;i;--i)
{
if(!inv[i]) inv[i]=_inv(i,p);
b[i]=mul(a[i-1],inv[i],p);
} b[0]=0;
}
inline void Der(int *a,int *b,int n)
{
for(int i=1;i<n;++i) b[i-1]=mul(a[i],i,p);
b[n-1]=0;//attention to the bound
}
inline void Ln(int *a,int *b,int n)
{
int l=0,lim=1;
while(lim<(n<<1)) lim<<=1,++l;
make_rev(lim,l);
memset(tmp,0,sizeof(int)*lim);
memset(tmpb,0,sizeof(int)*lim);
memcpy(tmp,a,sizeof(int)*n);
Der(a,tmp,n);
Inv(a,tmpb,n);
NTT(tmp,lim,1),NTT(tmpb,lim,1);
for(int i=0;i<lim;++i) tmp[i]=mul(tmp[i],tmpb[i],p);
NTT(tmp,lim,-1);
Int(tmp,tmp,n);
memcpy(b,tmp,sizeof(int)*n);
}
inline void Exp(int *a,int *b,int n)
{
if(n==1) { b[0]=1; return; }
Exp(a,b,(n+1)>>1);
Ln(b,tmp,n);
int l=0,lim=1;
while(lim<(n<<1)) lim<<=1,++l;
make_rev(lim,l);
for(int i=0;i<n;++i) tmp[i]=sub(a[i],tmp[i],p);
clear(tmp,n,lim);
clear(b,n,lim);
tmp[0]++;
NTT(tmp,lim,1);
NTT(b,lim,1);
arrMul(b,tmp,lim);
NTT(b,lim,-1);
clear(b,n,lim);
}
inline void Pow(int *a,int k,int *b,int n)
{
static int tmp[maxlen];
memset(tmp,0,sizeof(int)*(n<<1));
Ln(a,tmp,n);
arrMul(tmp,k,n);
Exp(tmp,b,n);
}//it's not totally safe
inline void Mod(int *a,int *b,int *P,int *Q,int n,int m)
{
--n,--m;//强行钦定使用左闭右开区间
static int ar[maxlen],br[maxlen];
memset(ar,0,sizeof(int)*(n));
memset(br,0,sizeof(int)*(n));
Reverse(a,ar,n+1);
Reverse(b,br,m+1);
irep(i,n-m+1,m) br[i]=0;
Inv(br,P,n-m+1);
pii lm=getLimit(2*n-m+1);//ar*p=br^-1
make_rev(lm.first,lm.second);
NTT(ar,lm.first,1),
NTT(P,lm.first,1);
arrMul(ar,P,lm.first);
NTT(ar,lm.first,-1);
Reverse(ar,P,n-m+1);//所以p不需要ntt回来
lm=getLimit(n+1);
make_rev(lm.first,lm.second);
clear(P,n-m+1,lm.first);
memcpy(ar,b,sizeof(int)*lm.first);
NTT(P,lm.first,1),
NTT(ar,lm.first,1);
arrMul(ar,P,lm.first);
NTT(P,lm.first,-1);//p is needed
NTT(ar,lm.first,-1);
for(int i=0;i<lm.first;++i) Q[i]=sub(a[i],ar[i],p);
clear(Q,m,lm.first);
}//[a]=n,[b]=m,a%b=p...q
poly():p(998244353),g(3),gi(332748118) { init(); }
poly(int mod):p(mod),g(set_mod::getg(mod)),gi(_inv(g,mod)) { init(); }
};
}
\(\text{poly}\) 提供了一个变幻模数的接口,不过使得常数变大了一倍,如果没有必要写任意模数的话应该拆开
函数的规范是 传入参数,传出参数,数组长度
比如求 \(g(x)=f^-1(x)\pmod {x^n}\) 就是 Inv(f,g,n)
,不难注意到这样写我们的数组是左闭右开的,事实上,我们要求所有函数这样(所以取模里面硬减了一下来)
更详细的约定是
-
NTT(a,b,c)
,a
为需要变幻的多项式,b
表示次数,c
为 \(-1\) 时进行逆变换,其他时候进行正变幻 -
Inv(a,b,c)
,a
为需要求逆的多项式,b
表示逆元(在 \(\mathbb{Z}_p\) 的意义下),c
为次数,之后默认最后一个参数表示数列长度(也意味着是在 \(\bmod {x^n}\) 的意义下) -
clear(a,b,c)
,清空a
在 \([l,r)\) 的值 -
arrAdd(a,b,c)
,如果b
为数列,则将对应位相加,如果b
为常数,则a
的每一位加上b
-
arrSub(a,b,c)
,同上,加法变为减法 -
arrMul(a,b,c)
,同上,加法变为乘法 -
Sqrt(a,b,c)
,求解 \(b\equiv\sqrt a\pmod{x^n}\) -
Der(a,b,c)
,求解 \(b=a^\prime\) -
Int(a,b,c)
,求解 \(b=\int a\),常数项置 \(0\) -
Ln(a,b,c)
,求解 \(b=\ln a\),应当有 \(a_0=1\) -
Exp(a,b,c)
,求解 \(b=e^x\),应当有 \(a_0=0\) -
Pow(a,k,b,c)
,求解 \(b=a^k\),应当有 \(a_0=1\) -
Mod(a,b,p,q,n,m)
,a
的长度为 \(n\),b
的长度为 \(m\),求解 \(c=\lfloor\dfrac{a}{b}\rfloor,d=a\bmod b\),就是小学学的大除法
原理都不难,主要就是多项式牛顿迭代
多项式牛迭是为了解决下面这个问题
我们期望用 \(G(x)\) 表示或者求出 \(F(x)\),当然,它们都是多项式
考虑一个倍增的做法
我们假设已经求出了 \(F_0(x)\) 满足
显然有
把 \(G(F(x))\) 看做关于 \(F(x)\) 的函数,在 \(F_0(x)\) 处展开
可得
我们只关心 \(F(x)\) 在 \(\bmod x^n\) 的意义下的值,而显然有 \(F(x)-F_0(x)\equiv0\pmod{x^\frac{n}{2}}\)
也就是说,\((F(x)-F_0(x))^2\equiv0\pmod {x^n}\)
于是只需要保留泰勒展开的前两项,即有
最后 \(n=1\) 的边界一般都可以转化为数论问题,然后就可以倍增了
复杂度一般都是 \(O(n\log n)\),用主定理看一眼,都是 \(T(n)=O(n\log n)+T(\dfrac{n}{2})\),当然,还有一些这里没有的操作,有些就不是这个复杂度,如果你写半在线卷积或者分治 \(\text{NTT}\) 的话也不是
符号化方法
实例
或许先看一些例子对于理解理论是有好处的
我们期望用生成函数解决组合问题,于是需要关联生成函数和组合意义,假定需要求解的组合对象构成的集合为 \(\mathcal{A}\),我们定义一个函数把组合对象映射到非负整数,这样就完成了连接生成函数和组合概念的第一步
我们规定 \(\alpha\in\mathcal{A}\),\(|\alpha|\) 为其大小函数,表示从组合对象到自然数的映射
比如我们考虑求解 0/1 背包问题拿总重量为 \(V_0\) 的方案数。假设对于单个物品,它的重量为 \(V\),那么我们可以直接把重量当做这个物品的大小函数,如果没有选择,那么大小自然是 \(0\)。于是我们可以写出单个物品的生成函数 \(1+x^V\),考虑每个组合对象(这里就是每个物品),它们的组合类应该是 \(\mathcal{A}_i=1+x^{V_i}\),最后对于整个组合对象,也就是我们要求解的问题,应当有
其中,\(\times\) 表示笛卡尔积,可以感性的理解为把两个组合对象直接拼接起来,而对应的代数含义就是生成函数的卷积
引入中性对象 \(\xi=1\) 和元对象 \(\circ=x\),可以类比自然数中的 \(0\) 和 \(1\)
形式化地,我们用 \(\mathcal{A}\cong \mathcal{B}\) 表示 \(\mathcal{A}\) 和 \(\mathcal{B}\) 是同构的组合对象,显然我们有
SEQ
终于开始有用的东西了...
我们考虑 \(\operatorname{SEQ}(S)\) 生成了包括 \(\varnothing\) 在内的所有用 \(S\) 中的元素构成的序列
比如 \(\text{SEQ}(\{a,b\})=\{\varnothing,\{a\},\{b\},\{a,a\},\{a,b\},\{b,a\},\{b,b\}\cdots\}\)
很显然,我们可以发现 \(\text{SEQ}\) 只是把若干个集合中选出任意元素有序地拼接在一起,那么考虑把这个自然语言表述的过程符号化
回想一下上面的过程,我们推导出 \(\text{SEQ}(\mathcal{A})\) 的封闭形式是 \(\dfrac{1}{1-\mathcal{A}(x)}\) 的过程只依赖于代数技巧,这启发我们这种函数可以应用于更多的场景,而不是局限于某种特定的组合概念(这也是符号化方法的优点,它简洁地指明了不同组合意义的共性)
求有序有根树的生成函数
树一定有根,儿子之间的顺序有意义
之后的处理方式有很多,个人比较推荐像链接里那样用拉格朗日反演
MSET
比起 \(\text{SEQ}\),无序的情形似乎更加常见,于是就有了 \(\text{MSET}(S)\),其包含所有可能的无序的组合
比如 \(\text{MSET}(\{a,b\})=\{\varnothing,\{a\},\{b\},\{a,a\},\{a,b\},\{b,b\}\cdots\}\)
注意到 \(\text{SEQ}\) 中出现了 \(\{a,b\},\{b,a\}\),而在 \(\text{MSET}\) 中,它们被认为是一样的
为了求解 \(\text{MSET}\),我们考虑递推,递推就是用我们已经知道的去求解不知道的,我们注意到 \(\text{MSET}(\{a\})=\text{SEQ}(\{a\})=\{\varnothing,\{a\},\{a.a\},\cdots\}\)
于是可以有
累乘的形式当然不能算作封闭形式,考虑高中数学常见的套路先取 \(\ln\) 再取 \(\exp\),设 \(c_i\) 表示 \(|a_k|=i\) 的 \(a\) 的数量,那么
利用麦克劳林公式 \(\ln(1+z)=\sum\limits_{i\ge 1}\dfrac{(-1)^{i-1}z^i}{i}\)
此时,我们发现 \(S(x)=-\sum\limits_{i\ge 1}c_i\ln(1-x^i)\)
整理得
较为简单的例子
\(\text{MSET}\) 的一个应用实例是求解数的拆分,如果考虑用 DP 解决,比如 \(f_{i,j}\) 表示用到了 \(i\),和为 \(j\) 的方案数
那么不难列出 \(f_{i,j}=\sum f_{i-1,j-ki}\)
发现这和 \(\text{MSET}\) 的递推式没有区别,于是要求的就是 \(\text{MSET}(\mathbb{N})\)
更为直接的方式是直接写出目标的生成函数 \(A(x)=\prod\limits_{i\ge 1}\dfrac{1}{1-x^i}\),不过这样不够机械,有种组合意义的感觉
更为困难的例子
求大小为 \(n\) 的无编号无根树的数量 \(a_n\)
先求无编号有根树 \(\mathcal{T}\)
无编号有根树的子树也是无编号有根树,因为没有编号,所以儿子之间的顺序不重要,这完全等价于儿子间存在严格的偏序关系
于是
由论文中可知,无编号无根树的生成函数是 \(T(z)-\frac 12T(z)^2+\frac12T(z^2)\)
PSET
定义 \(\text{PSET}(\mathcal{S})\) 为 \(\mathcal{S}\) 的所有子集构成的集合
比如 \(\text{PSET}(\{a,b\})=\{\varnothing,\{a\},\{b\},\{a,b\}\}\)
考虑 \(\text{PSET}(\{a\})=1+x^{|a|}\)
由定义知
CYC
我们定义 \(CYC(\mathcal{C})\) 表示用 \(\mathcal{C}\) 中元素可以表示的不同的环排列的集合,其中 \(\mathcal{C}\) 一般是可重集合
有
或者
证明有点难,我看懂了再补,咕咕咕咕咕咕