多项式
_ANIG_
·
·
个人记录
多项式
写多项式题像嗑药,出多项式题像贩毒。
写多项式总结叫...?
FFT
如何计算多项式的卷积?
对于一般的系数表示法,相乘是困难的。但是对于点值,相乘的简单的。
并且 n+1 个点可以确定唯一的 n 次多项式。
所以有一个思路:
-
转点值表示法。
-
点值相乘。
-
转系数表示法。
第二个操作很简单,第一,三个操作较为麻烦。
有这么一个神奇的函数:
struct fs{
double w1,w2;
friend fs operator+(fs a,fs b){
return(fs){a.w1+b.w1,a.w2+b.w2};
}
friend fs operator-(fs a,fs b){
return(fs){a.w1-b.w1,a.w2-b.w2};
}
friend fs operator*(fs a,fs b){
return(fs){a.w1*b.w1-a.w2*b.w2,a.w1*b.w2+a.w2*b.w1};
}
};
int r[10000005];
const double Pi=3.1415926535;
void Fast_Fast_TLE(fs*p,int n,int op){
int bt=1;
while((1<<bt)<n)bt++;
for(int i=0;i<n;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<bt-1);
if(i<r[i])swap(p[i],p[r[i]]);
}
for(int mid=1;mid<n;mid<<=1){
fs tmp=(fs){cos(Pi/mid),op*sin(Pi/mid)};
for(int i=0;i<n;i+=mid<<1){
fs om=(fs){1,0};
for(int j=0;j<mid;j++,om=om*tmp){
fs x=p[i+j],y=om*p[i+j+mid];
p[i+j]=x+y;
p[i+j+mid]=x-y;
}
}
}
}
函数的第三项,传进去 1,就是系数转点值,传进去 -1,就是点值转系数。
最上面是一个维护虚数的结构体,支持加减乘操作。
降低复杂度的原因是虚数平面内单位圆上的点的相反数可以通过单位根的若干次幂得到,从而减少运算。
NTT
FFT 无法进行取模,且有爆精度的可能性。对于一个能表示为 a×2^k+1 (k 较大)的质数,可以使用 NTT。
常见的 NTT 质数有 998244353,1004535809,469762049。
取模数的原根(以上三个质数的原根都是 3),原根有与单位复根相似的性质,所以可以代替单位复根,过程与 FFT 类似。
void ntt(int*p,int n,int op){
int bt=1;
while((1<<bt)<n)bt++;
for(int i=0;i<n;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<bt-1);
if(i<r[i])swap(p[i],p[r[i]]);
}
for(int mid=1;mid<n;mid<<=1){
int tmp=pows(op==1?3:332748118,(mods-1)/(mid<<1));
for(int i=0;i<n;i+=mid<<1){
int om=1;
for(int j=0;j<mid;j++,om=om*tmp%mods){
int x=p[i+j],y=om*p[i+j+mid]%mods;
p[i+j]=(x+y)%mods,p[i+j+mid]=(x-y+mods)%mods;
}
}
}
}
F_n=\sum\limits_{i=0}^nf_ig_{n-i}。
可以发现,把 f,g 的每一项作为多项式的系数,f,g 的乘积对应项的系数就是 F,直接 NTT 就能求 F。
实际做题中,往往把含 i 的分在一边,含 n-i 的分在一边,使式子转化为此形式。有时候也可以翻转数组来达到效果。
第二类斯特林数·行
令 f_{n,m} 表示第 n 行第 m 列的斯特林数。
通项公式:f_{n,m}=\sum\limits_{i=0}^m\frac{(-1)^{m-i}i^n}{i!(m-i)!}。
把含 i 的分一组,含 m-i 的分一组。
原式 =\sum\limits_{i=0}^m\frac{i^n}{i!}×\frac{(-1)^{m-i}}{(m-i)!}。
令 f1(i)=\frac{i^n}{i!},f2(i)=\frac{(-1)^i}{i!}。
则:f_{n,m}=(f1×f2)(m)。
NTT 卷积即可。
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int mods=167772161,n1=3,n2=55924054;
int pows(int a,int b){
if(b==0)return 1;
int res=pows(a,b>>1);
res=res*res%mods;
if(b&1)res=res*a%mods;
return res;
}
int r[1000005],jc[1000005];
void fft(int*p,int n,int op){
int bt=1;
while((1<<bt)<n)bt++;
for(int i=0;i<n;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<bt-1);
if(i<r[i])swap(p[i],p[r[i]]);
}
for(int mid=1;mid<n;mid<<=1){
int tmp=pows(op==1?n1:n2,(mods-1)/(mid<<1));
for(int i=0;i<n;i+=mid<<1){
int om=1;
for(int j=0;j<mid;j++,om=om*tmp%mods){
int x=p[i+j],y=p[i+j+mid]*om%mods;
p[i+j]=(x+y)%mods,p[i+j+mid]=(x-y)%mods;
}
}
}
}
int gets(int x){
int res=1;
while((1<<res)<=x)res++;
return 1<<res;
}
int n,f1[1000005],f2[1000005];
signed main(){
cin>>n;
jc[0]=1;
for(int i=1;i<=n;i++)jc[i]=jc[i-1]*i%mods;
for(int i=0;i<=n;i++)f1[i]=pows(-1,i)*pows(jc[i],mods-2),f2[i]=pows(i,n)*pows(jc[i],mods-2)%mods;
fft(f1,gets(2*n),1);
fft(f2,gets(2*n),1);
for(int i=0;i<=gets(2*n);i++)f1[i]=f1[i]*f2[i]%mods;
fft(f1,gets(2*n),-1);
int ny=pows(gets(2*n),mods-2);
for(int i=0;i<=n;i++){
printf("%lld ",(f1[i]*ny%mods+mods)%mods);
}
}
Tyvj1953 Normal
对于每个点分开考虑。这个点对答案的贡献是期望在多少次递归后被删除。
假设当前考虑的节点是根节点。
如果删除树上一个节点,树会分成若干个连通块。我们要考虑根节点所在的连通块。显然,这个连通块是被删除节点的父节点所处的连通块。也就是说删一个点等价于删除这个点的子树。
这样,题目就转化为:每次删除一个未被删除的节点的子树,求删除根节点的期望次数。
此时,题目就和 Game on tree 一模一样。
考虑每个点的贡献。如果这个点在操作序列中在所有祖先的前面,就对答案有贡献。概率就是祖先的个数的倒数,也就是深度的倒数。
答案就是所有节点的深度的倒数和,也就是所有点到根结点的距离的倒数和。
把每个点作为根节点考虑,把答案加起来,可以发现答案就是任意两个点的距离加 1 的倒数和。
考虑点分治。
求穿过根节点的路径长度加 1 的倒数和。
预处理子树中每个节点的深度,答案就是任意两个不在同一个子树内的点的深度和加 1 的倒数和。可以容斥。答案就是任意两个点的深度和加 1 的倒数和减去每个子树内深度和加 1 的倒数和。
问题转化为:\sum\limits_{i=1}^n\sum\limits_{j=1}^n\frac1{dep_i+dep_j+1}。
统计深度为 i 的点的数量 c_i,答案为:
\sum\limits_{i=1}^m\sum\limits_{j=1}^m\frac1{i+j+1}c_ic_j
枚举 $i+j$,答案为
$\sum\limits_{k=1}^{2m}\frac1{k+1}\sum\limits_{i=1}^kc_ic_{k-i}$。
右边的东西可以直接 FFT。
时间复杂度 $O(n\log^2n)$。
```cpp
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int mods=998244353;
int pows(int a,int b){
if(b==0)return 1;
int res=pows(a,b>>1);
res=res*res%mods;
if(b&1)res=res*a%mods;
return res;
}
const int n1=3,n2=pows(3,mods-2);
int r[1000005];
void fft(int *p,int n,int op){
int bt=1;
while((1<<bt)<n)bt++;
for(int i=0;i<n;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<bt-1);
if(i<r[i])swap(p[i],p[r[i]]);
}
for(int mid=1;mid<n;mid<<=1){
int tmp=pows(op==1?n1:n2,(mods-1)/(mid<<1));
for(int i=0;i<n;i+=mid<<1){
int om=1;
for(int j=0;j<mid;j++,om=om*tmp%mods){
int x=p[i+j],y=p[i+j+mid]*om%mods;
p[i+j]=(x+y)%mods,p[i+j+mid]=(x-y)%mods;
}
}
}
}
int gets(int x){
int res=1;
while((1<<res)<=x)res++;
return 1<<res;
}
double res;
int n,siz[100005],zd[100005],zx,mk[100005],bk[100005],dep[100005],c[100005],cr[100005],md,mdx;
vector<int>p[100005],jl;
void dfs1(int x,int y){
zd[x]=0;
mk[x]=1;siz[x]=1;
for(int i=0;i<p[x].size();i++){
int c=p[x][i];
if(mk[c]||bk[c])continue;
dfs1(c,y);
siz[x]+=siz[c];
zd[x]=max(zd[x],siz[c]);
}
zd[x]=max(zd[x],y-siz[x]);
if(zx==0||zd[zx]>zd[x])zx=x;
mk[x]=0;
}
void dfs2(int x){
mk[x]=1;siz[x]=1;
c[dep[x]]++;cr[dep[x]]++;
md=max(md,dep[x]);
mdx=max(mdx,dep[x]);
jl.push_back(x);
for(int i=0;i<p[x].size();i++){
int c=p[x][i];
if(mk[c]||bk[c])continue;
dep[c]=dep[x]+1;
dfs2(c);
siz[x]+=siz[c];
}
mk[x]=0;
}
void solve(int x){
if(siz[x]==1)return;
bk[x]=1;md=0;siz[x]=1;
for(int i=0;i<p[x].size();i++){
int c=p[x][i];
if(bk[c])continue;
dep[c]=1;mdx=0;
dfs2(c);
int m=gets(2*mdx);
fft(cr,m,1);
for(int i=0;i<=m;i++)cr[i]=cr[i]*cr[i]%mods;
fft(cr,m,-1);
int ny=pows(m,mods-2);
for(int i=0;i<=m;i++)cr[i]=cr[i]*ny%mods,cr[i]=(cr[i]+mods)%mods;
for(int i=1;i<=m;i++)res-=1.0/(i+1)*cr[i];
for(int i=0;i<=m;i++)cr[i]=0;
siz[x]+=siz[c];
}
int m=gets(2*md);
c[0]=1;
fft(c,m,1);
for(int i=0;i<=m;i++)c[i]=c[i]*c[i]%mods;
fft(c,m,-1);
int ny=pows(m,mods-2);
for(int i=0;i<=m;i++)c[i]=c[i]*ny%mods,c[i]=(c[i]+mods)%mods;
for(int i=1;i<=m;i++)res+=1.0/(i+1)*c[i];
for(int i=0;i<=m;i++)c[i]=0;
for(int i=0;i<p[x].size();i++){
int c=p[x][i];
if(bk[c])continue;
zx=0;
dfs1(c,siz[x]);
solve(zx);
}
}
signed main(){
cin>>n;
for(int i=1;i<n;i++){
int x,y;
scanf("%lld%lld",&x,&y);
x++,y++;
p[x].push_back(y);
p[y].push_back(x);
}
solve(1);
printf("%.4lf",res+n);
}
```
# 任意模数 NTT
有的毒瘤题会让你取模一个非 NTT 模数。此时需要任意模数 NTT。
一般取三个 NTT 模数 $p1,p2,p3$,分别算出来模这三个数的结果。
用 exCRT 算出模 $p1\times p2\times p3$ 的结果。
大部分情况下,$p1\times p2\times p3$ 是大于真实系数的。所以算出来的结果就是实际的系数。把这个结果模题目要求的模数即可。
具体的,设:
$x\equiv a\pmod{p1}
x\equiv b\pmod{p2}
x\equiv c\pmod{p3}
先考虑前两个方程,设 x=kp1+a,则 kp1+a\equiv b\pmod{p2}。
kp1\equiv b-a\pmod{p2}
k\equiv (b-a)p1^{-1}\pmod{p2}
快速幂求逆元即可。设 x'=kp1+a,x=Km+x'(m=p1p2)。
再考虑第三个方程:
Km+x'\equiv c\pmod{p3}
K\equiv (c-x')\times m^{-1}\pmod{p3}
再求逆元解出来 K,答案就是 Km+x'。
#include <bits/stdc++.h>
using namespace std;
#define int long long
int pows(int a,int b,int mods){
if(b<0)return pows(pows(a,mods-2,mods),-b,mods);
if(b==0)return 1;
int res=pows(a,b>>1,mods);
res=res*res%mods;
if(b&1)res=res*a%mods;
return res;
}
const int mods1=998244353,mods2=1004535809,mods3=469762049,ny1=pows(mods1,mods2-2,mods2),ny2=pows(mods1*mods2%mods3,mods3-2,mods3);
int mods;
struct node{
int a,b,c,d;
friend node operator+(node a,node b){return (node){(a.a+b.a)%mods1,(a.b+b.b)%mods2,(a.c+b.c)%mods3,(a.d+b.d)%mods};}
friend node operator-(node a,node b){return (node){(a.a-b.a)%mods1,(a.b-b.b)%mods2,(a.c-b.c)%mods3,(a.d-b.d)%mods};}
friend node operator*(node a,node b){return (node){a.a*b.a%mods1,a.b*b.b%mods2,a.c*b.c%mods3,a.d*b.d%mods};}
friend node operator*(node a,int b){return (node){a.a*b%mods1,a.b*b%mods2,a.c*b%mods3,a.d*b%mods};}
void crt(){
a=(a+mods1)%mods1,b=(b+mods2)%mods2,c=(c+mods3)%mods3;
int k1=(b-a)*ny1%mods2,x=k1*mods1+a,k2=(c-x%mods3+mods3)%mods3*ny2%mods3;
d=x%mods+k2*mods1%mods*mods2%mods;
d=(d%mods+mods)%mods;
}
void rst(){a=b=c=d;}
}f1[1000005],f2[1000005],f[1000005];
node pows(node a,int b){
return (node){pows(a.a,b,mods1),pows(a.b,b,mods2),pows(a.c,b,mods3),pows(a.d,b,mods)};
}
int r[1000005];
void fft(node*p,int n,int op){
int bt=1;
while((1<<bt)<n)bt++;
for(int i=0;i<n;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<bt-1);
if(i<r[i])swap(p[i],p[r[i]]);
}
for(int mid=1;mid<n;mid<<=1){
node tmp=(node){pows(3,op*(mods1-1)/(mid<<1),mods1),pows(3,op*(mods2-1)/(mid<<1),mods2),pows(3,op*(mods3-1)/(mid<<1),mods3)};
for(int i=0;i<n;i+=mid<<1){
node om=(node){1,1,1,1};
for(int j=0;j<mid;j++,om=om*tmp){
node x=p[i+j],y=om*p[i+j+mid];
p[i+j]=x+y,p[i+j+mid]=x-y;
}
}
}
}
int gets(int x){
int res=1;
while((1<<res)<=x)res++;
return 1<<res;
}
int n,m,p1[1000005],p2[1000005];
signed main(){
cin>>n>>m>>mods;
for(int i=0;i<=n;i++)scanf("%lld",&p1[i]),f1[i]=(node){p1[i],p1[i],p1[i],p1[i]};
for(int i=0;i<=m;i++)scanf("%lld",&p2[i]),f2[i]=(node){p2[i],p2[i],p2[i],p2[i]};
int mm=gets(n+m);
fft(f1,mm,1);fft(f2,mm,1);
for(int i=0;i<=mm;i++)f[i]=f1[i]*f2[i];
fft(f,mm,-1);
node ny=pows((node){mm,mm,mm,mm},-1);
for(int i=0;i<=mm;i++)f[i]=f[i]*ny,f[i].crt();
for(int i=0;i<=n+m;i++)printf("%lld ",(f[i].d+mods)%mods);
}
Evaluation
设答案为 F_k。
F_k=\sum\limits_{i=0}^{n-1}a_i(b\times c^{2k}+d)^i
=\sum\limits_{i=0}^{n-1}a_i\sum\limits_{j=0}^iC_j^i(b\times c^{2k})^jd^{i-j}
=\sum\limits_{i=0}^{n-1}a_i\sum\limits_{j=0}^i\frac{i!(b\times c^{2k})^jd^{i-j}}{j!(i-j)!}
令负数的阶乘的逆元为 0。
F_k=\sum\limits_{i=0}^{n-1}a_i\sum\limits_{j=0}^{n-1}\frac{i!(b\times c^{2k})^jd^{i-j}}{j!(i-j)!}
=\sum\limits_{j=0}^{n-1}\sum\limits_{i=0}^{n-1}\frac{a_ii!(b\times c^{2k})^jd^{i-j}}{j!(i-j)!}
=\sum\limits_{j=0}^{n-1}\frac{(b\times c^{2k})^j}{j!}\sum\limits_{i=0}^{n-1}\frac{a_ii!d^{i-j}}{(i-j)!}
右边的东西可以 NTT 预处理,设为 f_j。
F_k=\sum\limits_{j=0}^{n-1}\frac{(b\times c^{2k})^jf_j}{j!}
设 $\frac{b^jf_j}{j!}=g_j$。
$F_k=\sum\limits_{j=0}^{n-1}g_jc^{2jk}$。
这个貌似不好处理。有一个很妙的转化:
$F_k=\sum\limits_{j=0}^{n-1}g_jc^{j^2+k^2-(k-j)^2}=c^{k^2}\sum\limits_{j=0}^{n-1}\frac{g_jc^{j^2}}{c^{(k-j)^2}}
这就又是卷积的形式,直接 NTT。
还要任意模数,细节很多。
#include <bits/stdc++.h>
using namespace std;
#define int long long
int pows(int a,int b,int mods){
if(b<0)return pows(pows(a,mods-2,mods),-b,mods);
if(b==0)return 1;
int res=pows(a,b>>1,mods);
res=res*res%mods;
if(b&1)res=res*a%mods;
return res;
}
const int mods1=998244353,mods2=1004535809,mods3=469762049,ny1=pows(mods1,mods2-2,mods2),ny2=pows(mods1*mods2%mods3,mods3-2,mods3),mods=1e6+3;
struct node{
signed a,b,c,d;
friend node operator+(node a,node b){return (node){(a.a+b.a)%mods1,(a.b+b.b)%mods2,(a.c+b.c)%mods3,(a.d+b.d)%mods};}
friend node operator-(node a,node b){return (node){(a.a-b.a)%mods1,(a.b-b.b)%mods2,(a.c-b.c)%mods3,(a.d-b.d)%mods};}
friend node operator*(node a,node b){return (node){1ll*a.a*b.a%mods1,1ll*a.b*b.b%mods2,1ll*a.c*b.c%mods3,1ll*a.d*b.d%mods};}
friend node operator*(node a,int b){return (node){1ll*a.a*b%mods1,1ll*a.b*b%mods2,1ll*a.c*b%mods3,1ll*a.d*b%mods};}
void crt(){
a=(a+mods1)%mods1,b=(b+mods2)%mods2,c=(c+mods3)%mods3;
int k1=(b-a)*ny1%mods2,x=k1*mods1+a,k2=(c-x%mods3+mods3)%mods3*ny2%mods3;
d=x%mods+k2*mods1%mods*mods2%mods;
d=(d%mods+mods)%mods;
a=b=c=d;
}
}f1[600005],f2[600005],f[600005],jc[600005],ny[600005],g[600005];
node pows(node a,int b){
return (node){pows(a.a,b,mods1),pows(a.b,b,mods2),pows(a.c,b,mods3),pows(a.d,b,mods)};
}
int r[1000005];
void fft(node*p,int n,int op){
int bt=1;
while((1<<bt)<n)bt++;
for(int i=0;i<n;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<bt-1);
if(i<r[i])swap(p[i],p[r[i]]);
}
for(int mid=1;mid<n;mid<<=1){
node tmp=(node){pows(3,op*(mods1-1)/(mid<<1),mods1),pows(3,op*(mods2-1)/(mid<<1),mods2),pows(3,op*(mods3-1)/(mid<<1),mods3)};
for(int i=0;i<n;i+=mid<<1){
node om=(node){1,1,1,1};
for(int j=0;j<mid;j++,om=om*tmp){
node x=p[i+j],y=om*p[i+j+mid];
p[i+j]=x+y,p[i+j+mid]=x-y;
}
}
}
}
int gets(int x){
int res=1;
while((1<<res)<=x)res++;
return 1<<res;
}
int n,b,c,d,p[1000005];
signed main(){
cin>>n>>b>>c>>d;
int m=gets(3*n);
jc[0]=ny[0]=(node){1,1,1,1};
for(int i=1;i<=n+1;i++)jc[i]=jc[i-1]*i,ny[i]=pows(jc[i],-1);
for(int i=0;i<n;i++)scanf("%lld",&p[i]);
for(int i=0;i<=n;i++)f1[i].d=p[i]*jc[i].d%mods,f1[i].a=f1[i].b=f1[i].c=f1[i].d;
for(int i=0;i<=n+1;i++)f2[n-i+1].d=1ll*pows(d,i,mods)*ny[i].d%mods,f2[n-i+1].a=f2[n-i+1].b=f2[n-i+1].c=f2[n-i+1].d;
fft(f1,m,1);fft(f2,m,1);
for(int i=0;i<=m;i++)f[i]=f1[i]*f2[i];
fft(f,m,-1);
node nys=pows((node){m,m,m,m},-1);
for(int i=0;i<=m;i++)f[i]=f[i]*nys,f[i].crt();
for(int i=0;i<n;i++)g[i].d=1ll*ny[i].d*f[n+i+1].d%mods*pows(b,i,mods)%mods,g[i].a=g[i].b=g[i].c=g[i].d;
memset(f1,0,sizeof(f1));
memset(f2,0,sizeof(f2));
memset(f,0,sizeof(f));
for(int i=0;i<n;i++)f1[i].d=1ll*g[i].d*pows(c,i*i,mods)%mods,f1[i].a=f1[i].b=f1[i].c=f1[i].d;
for(int i=-n+1;i<=n;i++)f2[n+i].d=pows(c,-i*i,mods),f2[n+i].a=f2[n+i].b=f2[n+i].c=f2[n+i].d;
fft(f1,m,1);fft(f2,m,1);
for(int i=0;i<=m;i++)f[i]=f1[i]*f2[i];
fft(f,m,-1);
for(int i=0;i<=m;i++)f[i]=f[i]*nys,f[i].crt();
for(int i=0;i<n;i++){
int res=1ll*pows(c,i*i,mods)*f[n+i].d%mods;
printf("%lld\n",res);
}
}
多项式求逆
递归求解 f 的逆元。
设 hf\equiv 1\pmod{\frac n2}
要求 gf\equiv 1\pmod n。
所以 g-h\equiv 0\pmod{ \frac n2}。
(g-h)^2\equiv 0\pmod n
g^2-2gh+h^2\equiv0\pmod n
g^2f-2ghf+h^2f\equiv 0\pmod n
g-2h+h^2f\equiv 0\pmod n
于是得到了 $g$ 的递推式。
复杂度 $O(n\log n)$。
## [ [集训队作业2013] 城市规划 ](https://www.luogu.com.cn/problem/P4841)
求 $n$ 个点的无向连通图的数量。
设 $f_n$ 为 $n$ 个点无向连通图的数量,$g_n$ 为 $n$ 个点无向图的数量。
显然,$g_n=2^{\frac{n\times(n-1)}2}$。
考虑容斥,$f_n$ 就是 $n$ 个点的无向图数量减去 $n$ 个点不连通图的数量。
后面的可以枚举 $1$ 号点所在的连通块大小,剩余的图就是一个任意的无向图。
$f_n=g_n-\sum\limits_{i=1}^{n-1}C_{n-1}^{i-1}f_ig_{n-i}$。
$g_n-\sum\limits_{i=1}^{n-1}C_{n-1}^{i-1}f_ig_{n-i}-f_n=0
g_n-\sum\limits_{i=1}^{n}C_{n-1}^{i-1}f_ig_{n-i}=0
令 f_0=0。
g_n-\sum\limits_{i=0}^{n}C_{n-1}^{i-1}f_ig_{n-i}=0
g_n=\sum\limits_{i=0}^{n}C_{n-1}^{i-1}f_ig_{n-i}
把组合数展开。
g_n=\sum\limits_{i=0}^n\frac{(n-1)!f_ig_{n-i}}{(i-1)!(n-i)!}
g_n=(n-1)!\sum\limits_{i=0}^n\frac{f_i}{(i-1)!}\times\frac{g_{n-i}}{(n-i)!}
\frac{g_n}{(n-1)!}=\sum\limits_{i=0}^n\frac{f_i}{(i-1)!}\times\frac{g_{n-i}}{(n-i)!}
令 F_x=\frac{f_x}{(x-1)!},G_x=\frac{g_i}{i!},H_x=\frac{g_i}{(i-1)!}。
则 H=F\times G。
直接求逆即可。
```cpp
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int mods=1004535809,n1=3,n2=334845270;
int pows(int a,int b){
if(b==0)return 1;
int res=pows(a,b>>1);
res=res*res%mods;
if(b&1)res=res*a%mods;
return res;
}
int r[1000005],f[1000005],n,jc[1000005],ny[1000005],ps[1000005],g[1000005],p[1000005];
void fft(int*p,int n,int op){
int bt=1;
while((1<<bt)<n)bt++;
for(int i=0;i<n;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<bt-1);
if(i<r[i])swap(p[i],p[r[i]]);
}
for(int mid=1;mid<n;mid<<=1){
int tmp=pows(op==1?n1:n2,(mods-1)/(mid<<1));
for(int i=0;i<n;i+=mid<<1){
int om=1;
for(int j=0;j<mid;j++,om=om*tmp%mods){
int x=p[i+j],y=p[i+j+mid]*om%mods;
p[i+j]=(x+y)%mods,p[i+j+mid]=(x-y)%mods;
}
}
}
}
int gets(int x){
int res=1;
while((1<<res)<=x)res++;
return 1<<res;
}
int C(int n,int m){
if(n<m)return 0;
return jc[n]*ny[m]%mods*ny[n-m]%mods;
}
void solve(int n){
if(n==1){
g[0]=pows(p[0],mods-2);
g[1]=pows(p[1],mods-2);
return;
}
solve(n+1>>1);
memset(f,0,sizeof(f));
for(int i=0;i<=n;i++)f[i]=p[i];
fft(f,gets(2*n),1);
fft(g,gets(2*n),1);
for(int i=0;i<=gets(2*n);i++)g[i]=(2*g[i]-f[i]*g[i]%mods*g[i]%mods)%mods;
fft(g,gets(2*n),-1);
int ny=pows(gets(2*n),mods-2);
for(int i=0;i<=gets(2*n);i++)g[i]=g[i]*ny%mods;
for(int i=n+1;i<=gets(2*n);i++)g[i]=0;
}
signed main(){
cin>>n;
ny[0]=1;
jc[0]=1;
for(int i=1;i<=n;i++)jc[i]=jc[i-1]*i%mods,ny[i]=pows(jc[i],mods-2);
for(int i=0;i<=n;i++)p[i]=pows(2,i*(i-1)/2)*ny[i]%mods,ps[i]=pows(2,i*(i-1)/2)*ny[i-1]%mods;
solve(n);
fft(g,gets(2*n),1);
fft(ps,gets(2*n),1);
memset(f,0,sizeof(f));
for(int i=0;i<=gets(2*n);i++)f[i]=g[i]*ps[i]%mods;
fft(f,gets(2*n),-1);
cout<<(f[n]*pows(gets(2*n),mods-2)%mods*jc[n-1]%mods+mods)%mods;
}
```
# 其他操作
## 多项式 ln
设 $h=\ln f$。
则 $h'=\frac{f'}{f}$。
求导,求逆,积分即可。
## 多项式 exp
设原式为 $p(x)
设 h(f(x))=\ln f(x)-p(x)
牛顿迭代得
f(x)=f_0(x)-\frac{h(f_0(x))}{h'(f_0(x))}
f(x)=f_0(x)-\frac{\ln f_0(x)-p(x)}{\frac{1}{f_0(x)}}
f(x)=f_0(x)-f_0(x)(\ln f_0(x)-p(x))
f(x)=f_0(x)(1-\ln f_0(x)+p(x))
直接递归求解。复杂度 O(n\log n)。
多项式快速幂
设 f=p^n
\ln f=n\ln p
e^{\ln f}=e^{n\ln p}
f=e^{n\ln p}
复杂度 O(n\log n)。
多项式开根
设 g^2\equiv p\pmod \frac n2
f^2\equiv p\pmod n
f-g\equiv 0\pmod \frac n2
f^2-2fg+g^2\equiv0\pmod n
p-2fg+g^2\equiv 0\pmod n
f\equiv\frac{p+g^2}{2g}\pmod n
复杂度 O(n\log n)
The Child and Binary Tree
设 f_n 为权值为 n 的二叉树的个数,c_i 表示 i 是否在集合内。
枚举根节点的权值。
f_n=\sum\limits_{i=0}^{n}c_i\sum\limits_{j=0}^{n-i}f_jf_{n-i-j}+[n==0]
令 g=f^2。
所以 $f=c\times g+1$。
$f=c\times f^2+1
cf^2-f+1=0
f=\frac{1-\sqrt{1-4c}}{2c}=\frac{(1-\sqrt{1-4c})(1+\sqrt{1-4c})}{2c(1+\sqrt{1-4c})}=\frac{4c}{2c+2c\sqrt{1-4c}}=\frac{2}{1+1\sqrt{1-4c}}
开根解出来即可。
其他应用
套矩阵乘法
PolandBall and Many Other Balls
设 f_{n,k} 为前 n 个球分 k 组的方案数。
设 $f_n(x)$ 为第 $n$ 行代表的多项式。
$f_n=(x+1)f_{n-1}+xf_{n-2}$。
这就是矩阵乘法的形式。直接矩阵快速幂。
```cpp
#include <bits/stdc++.h>
using namespace std;
const int mods=998244353;
int r[1000005];
int pows(int a,int b){
if(b==0)return 1;
int res=pows(a,b>>1);
res=1ll*res*res%mods;
if(b&1)res=1ll*res*a%mods;
return res;
}
const int n1=3,n2=pows(3,mods-2);
void fft(int*p,int n,int op){
int bt=1;
while((1<<bt)<n)bt++;
for(int i=0;i<n;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<bt-1);
if(i<r[i])swap(p[i],p[r[i]]);
}
for(int mid=1;mid<n;mid<<=1){
int tmp=pows(op==1?n1:n2,(mods-1)/(mid<<1));
for(int i=0;i<n;i+=mid<<1){
int om=1;
for(int j=0;j<mid;j++,om=1ll*om*tmp%mods){
int x=p[i+j],y=1ll*p[i+j+mid]*om%mods;
p[i+j]=(x+y)%mods,p[i+j+mid]=(x-y)%mods;
}
}
}
}
int gets(int x){
int res=1;
while((1<<res)<=x)res++;
return 1<<res;
}
int n,k;
struct node{
int p[(1<<17)+5];
friend node operator*(node a,node b){
int m=gets(2*k);
fft(a.p,m,1);fft(b.p,m,1);
for(int i=0;i<=m;i++)a.p[i]=1ll*a.p[i]*b.p[i]%mods;
fft(a.p,m,-1);
int ny=pows(m,mods-2);
for(int i=0;i<=m;i++)a.p[i]=1ll*a.p[i]*ny%mods;
for(int i=k+1;i<=m;i++)a.p[i]=0;
return a;
}
friend node operator+(node a,node b){
for(int i=0;i<=k;i++)a.p[i]+=b.p[i],a.p[i]%=mods;
return a;
}
};
struct jz{
node w[2][2];
int n,m;
friend jz operator*(jz a,jz b){
jz c;
c.n=a.n,c.m=b.m;
for(int i=0;i<a.n;i++){
for(int j=0;j<b.m;j++){
memset(c.w[i][j].p,0,sizeof(c.w[i][j].p));
for(int k=0;k<a.m;k++){
c.w[i][j]=c.w[i][j]+a.w[i][k]*b.w[k][j];
}
}
}
return c;
}
}zy,emp,cs;
jz pows(jz a,int b){
jz res=emp,nw=a;
while(b){
if(b&1)res=res*nw;
nw=nw*nw;
b>>=1;
}
return res;
}
signed main(){
cin>>n>>k;
emp.n=emp.m=2;
zy.n=zy.m=2;
cs.n=1,cs.m=2;
emp.w[0][0].p[0]=emp.w[1][1].p[0]=1;
zy.w[0][1].p[1]=1;
zy.w[1][1].p[0]=1;
zy.w[1][1].p[1]=1;
zy.w[1][0].p[0]=1;
cs.w[0][1].p[0]=1;
cs.w[1][1].p[0]=1;
cs.w[1][1].p[1]=1;
jz res=cs*pows(zy,n);
for(int i=1;i<=k;i++)cout<<(res.w[0][1].p[i]+mods)%mods<<" ";
}
```
## 倍增 FFT
## [ Transforming Sequence ](https://www.luogu.com.cn/problem/CF623E)
转化题意,一个长度为 $k$ 的数组,$n$ 次操作,每次选择至少一个为 $0$ 的地方变成 $1$,问最后有多少种操作序列。
设 $f_{n,k}$ 为长度为 $n$ 的序列,最终有 $k$ 位为 $1$ 的方案数。
最终答案就是 $\sum\limits_{i=0}^kC_k^if_{n,i}$。
显然 $n>k$ 无解,直接输出 $0$。
考虑 $f_{a+b,k}$。
枚举这 $k$ 个 $1$,有多少个是前 $a$ 个数中的,此时后 $b$ 个数每个数多了 $i$ 种选择,所以贡献还要乘上 $2^{ib}$。
$f_{a+b,k}=\sum\limits_{i=0}^kC_k^i2^{ib}f_{a,i}f_{b,k-i}$。
这个可以直接卷起来。用类似快速幂的方法,每次从 $f_{\frac n2}$ 转移到 $f_n$,时间复杂度 $O(k\log^2 k)$。还要写任意模数。
```cpp
#include <bits/stdc++.h>
using namespace std;
#define int long long
int pows(int a,int b,int mods){
if(b<0)return pows(pows(a,mods-2,mods),-b,mods);
if(b==0)return 1;
int res=pows(a,b>>1,mods);
res=res*res%mods;
if(b&1)res=res*a%mods;
return res;
}
const int mods1=998244353,mods2=1004535809,mods3=469762049,mods=1e9+7,n1=pows(mods1,mods2-2,mods2),n2=pows(mods1*mods2%mods3,mods3-2,mods3);
struct node{
int a,b,c;
friend node operator+(node a,node b){
return (node){(a.a+b.a)%mods1,(a.b+b.b)%mods2,(a.c+b.c)%mods3};
}
friend node operator-(node a,node b){
return (node){(a.a-b.a)%mods1,(a.b-b.b)%mods2,(a.c-b.c)%mods3};
}
friend node operator*(node a,node b){
return (node){a.a*b.a%mods1,a.b*b.b%mods2,a.c*b.c%mods3};
}
void crt(){
a=(a+mods1)%mods1,b=(b+mods2)%mods2,c=(c+mods3)%mods3;
int k1=(b-a)%mods2*n1%mods2,x=k1*mods1+a,k2=(c-x%mods3+mods3)%mods3*n2%mods3;
a=b=c=((x%mods+k2*mods1%mods*mods2%mods)%mods+mods)%mods;
}
}f1[400005],f2[400005];
int r[400005];
void fft(node*p,int n,int op){
int bt=1;
while((1<<bt)<n)bt++;
for(int i=0;i<n;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<bt-1);
if(i<r[i])swap(p[i],p[r[i]]);
}
for(int mid=1;mid<n;mid<<=1){
node tmp=(node){pows(3,op*(mods1-1)/(mid<<1),mods1),pows(3,op*(mods2-1)/(mid<<1),mods2),pows(3,op*(mods3-1)/(mid<<1),mods3)};
for(int i=0;i<n;i+=mid<<1){
node om=(node){1,1,1};
for(int j=0;j<mid;j++,om=om*tmp){
node x=p[i+j],y=om*p[i+j+mid];
p[i+j]=x+y,p[i+j+mid]=x-y;
}
}
}
}
int gets(int x){
int res=1;
while((1<<res)<=x)res++;
return 1<<res;
}
node pows(node a,int b){
return (node){pows(a.a,b,mods1),pows(a.b,b,mods2),pows(a.c,b,mods3)};
}
void mul(int*a,int*b,int*g,int n){
int m=gets(2*n);
memset(f1,0,sizeof(f1));
memset(f2,0,sizeof(f2));
for(int i=0;i<=n;i++)f1[i]=(node){a[i],a[i],a[i]};
for(int i=0;i<=n;i++)f2[i]=(node){b[i],b[i],b[i]};
fft(f1,m,1);fft(f2,m,1);
for(int i=0;i<=m;i++)f1[i]=f1[i]*f2[i];
fft(f1,m,-1);
node ny=pows((node){m,m,m},-1);
for(int i=0;i<=m;i++)f1[i]=f1[i]*ny,f1[i].crt(),g[i]=f1[i].a;
}
int g1[400005],g2[400005],ny[400005],jc[400005];
void solve(int n,int k,int*p){
if(n==1){
for(int i=1;i<=k;i++)p[i]=1;
return;
}
solve(n>>1,k,p);
memset(g1,0,sizeof(g1));
memset(g2,0,sizeof(g2));
for(int i=0;i<=k;i++)g1[i]=p[i]*pows(2,i*(n>>1),mods)%mods*ny[i]%mods,g2[i]=p[i]*ny[i]%mods;
memset(p,0,sizeof(p));
mul(g1,g2,p,k);
for(int i=0;i<=k;i++)p[i]=p[i]*jc[i]%mods;
if(n&1){
memset(g1,0,sizeof(g1));
memset(g2,0,sizeof(g2));
for(int i=0;i<=k;i++)g1[i]=p[i]*pows(2,i,mods)%mods*ny[i]%mods,g2[i]=(i!=0)*ny[i]%mods;
memset(p,0,sizeof(p));
mul(g1,g2,p,k);
for(int i=0;i<=k;i++)p[i]=p[i]*jc[i]%mods;
}
}
int C(int a,int b){
return jc[a]*ny[b]%mods*ny[a-b]%mods;
}
int n,k,p[1000005],res;
signed main(){
cin>>n>>k;
if(n>k){
cout<<0;
return 0;
}
jc[0]=ny[0]=1;
for(int i=1;i<=k;i++)jc[i]=jc[i-1]*i%mods,ny[i]=pows(jc[i],mods-2,mods);
solve(n,k,p);
for(int i=1;i<=k;i++)res+=C(k,i)*p[i]%mods,res%=mods;
cout<<res;
}
```
## 分治 FFT
## [Tree Coloring ](https://www.luogu.com.cn/problem/CF1613F)
考虑容斥。
用 $f_i$ 表示选定了 $i$ 组不合法的方案数。考虑此时的方案数。
选定了 $i$ 组不合法,也就是有 $i$ 个点,满足它的权值是父节点权值减 $1$。也就是确定了父节点的权值就能确定这个点的权值。不妨把这两个点合并。这样,权值自由的点数减少到 $n-i$。考虑颜色序列,$w_i$ 表示颜色为 $i$ 的点的编号。那么就相当于确定了 $k$ 个连续的块。可以理解为把 $n-k$ 个有长度的木块排列组合,方案数 $(n-i)!$。
答案为 $\sum\limits_{i=0}(-1)^if_i(n-i)!$。
考虑求 $f_i$。每个非叶子节点,最多有一组不合法。也就是我们可以从非叶子节点内选择 $i$ 个点,再在每个节点的儿子节点内选一个点。这就类似第 $i$ 个节点有 $c_i$ 种,要选出来 $k$ 个,求方案数。($c_i$ 表示 $i$ 的子结点个数)
这样就转化成一个背包问题:每个点的体积为 $1$,种类为 $c_i$,求装满一个容量为 $k$ 的背包的方案数。
令 $dp_{n,k}$ 表示考虑编号 $1-n$ 中的点,目前选了 $k$ 个点的方案数。
$dp_{n,k}=c_ndp_{n-1,k-1}+dp_{n-1,k}$。
看作 $n$ 个多项式,则 $dp_n=(c_nx+1)dp_{n-1}$。
也就是 $dp_n=\prod\limits_{i=1}^n(c_ix+1)$。
直接乘是 $n^2$ 的。
可以分治,每层要卷的长度是这层的区间长度。总要卷的长度就是 $n\log n$。
总复杂度 $O(n\log^2n)$。
要使用动态开空间的数组,需要写 vector NTT。
```cpp
#include <bits/stdc++.h>
using namespace std;
const int mods=998244353;
int r[5000005];
int pows(int a,int b){
if(b==0)return 1;
int res=pows(a,b>>1);
res=1ll*res*res%mods;
if(b&1)res=1ll*res*a%mods;
return res;
}
const int n1=3,n2=pows(3,mods-2);
void fft(vector<int>&p,int n,int op){
int bt=1;
while((1<<bt)<n)bt++;
for(int i=0;i<n;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<bt-1);
if(i<r[i])swap(p[i],p[r[i]]);
}
for(int mid=1;mid<n;mid<<=1){
int tmp=pows(op==1?n1:n2,(mods-1)/(mid<<1));
for(int i=0;i<n;i+=mid<<1){
int om=1;
for(int j=0;j<mid;j++,om=1ll*om*tmp%mods){
int x=p[i+j],y=1ll*p[i+j+mid]*om%mods;
p[i+j]=(x+y)%mods,p[i+j+mid]=(x-y)%mods;
}
}
}
}
int gets(int x){
int res=1;
while((1<<res)<=x)res++;
return 1<<res;
}
vector<int>e[2000005],q,dp;
int mk[2000005],n,res,dxs[2000005],jc[2000005],ny[2000005],sm;
void dfs(int x){
mk[x]=1;
for(int i=0;i<e[x].size();i++){
int c=e[x][i];
if(mk[c])continue;
dfs(c);
dxs[x]++;
}
if(dxs[x])sm++,q.push_back(dxs[x]);
}
void mul(vector<int>&a,vector<int>&b,vector<int>&c,int n){
int m=gets(2*n);
for(int i=a.size();i<=m;i++)a.push_back(0);
for(int i=b.size();i<=m;i++)b.push_back(0);
c.clear();
for(int i=0;i<=m;i++)c.push_back(0);
fft(a,m,1);fft(b,m,1);
for(int i=0;i<=m;i++)c[i]=1ll*a[i]*b[i]%mods;
fft(c,m,-1);
int ny=pows(m,mods-2);
for(int i=0;i<=m;i++)c[i]=1ll*c[i]*ny%mods;
for(int i=n+1;i<=m;i++)c.pop_back();
}
void solve(int l,int r,vector<int>&g){
if(l==r){
g[0]=1;
g[1]=q[l-1];
return;
}
vector<int>g1,g2;
int mid=l+r>>1;
int m=min((r-l+1)*2,n);
for(int i=0;i<=4*m;i++)g1.push_back(0),g2.push_back(0);
solve(l,mid,g1);solve(mid+1,r,g2);
mul(g1,g2,g,m);
}
signed main(){
cin>>n;
jc[0]=ny[0]=1;
for(int i=1;i<=n;i++)jc[i]=1ll*jc[i-1]*i%mods,ny[i]=pows(jc[i],mods-2);
for(int i=1;i<n;i++){
int x,y;
scanf("%d%d",&x,&y);
e[x].push_back(y);
e[y].push_back(x);
}
dfs(1);
for(int i=0;i<=4*n;i++)dp.push_back(0);
solve(1,sm,dp);
for(int i=0;i<=sm;i++){
int ans=dp[i];
res+=1ll*pows(-1,i)*ans*jc[n-i]%mods;
res%=mods;
}
cout<<(res+mods)%mods;
}
```
## [有标号二分图计数](https://www.luogu.com.cn/problem/P7364)
很容易想到一个思路。
枚举左部点的个数,两个部分之间任意连边。
$f_n=\sum\limits_{i=0}^nC_n^i2^{i(n-i)}$。
但是会重复。
对于多个连通块,每个连通块都有两种染色方式,只应统计一次。
但是,如果只有一个连通块,就只会多统计一次,直接除以 $2$ 即可。
所以可以用类似统计连通无向图数的方法统计连通二分图数。
现在就变成了大小为 $i$ 的块有 $a_i$ 种,构成 $n$ 共有多少种方式。
枚举连通块个数。
答案为 $\sum\limits_{1\le i}\frac{f^i}{i!}=e^f$。
```cpp
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int mods=998244353,s2=116195171;
int pows(int a,int b){
if(b==0)return 1;
int res=pows(a,b>>1);
res=res*res%mods;
if(b&1)res=res*a%mods;
return res;
}
const int n1=3,n2=pows(3,mods-2),ny2=pows(s2,mods-2);
int n,r[1000005],jc[1000005],ny[1000005],g[1000005],f[1000005],gs[1000005],fs[1000005],rs[1000005];
void fft(int*p,int n,int op){
int bt=1;
while((1<<bt)<n)bt++;
for(int i=0;i<n;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<bt-1);
if(i<r[i])swap(p[i],p[r[i]]);
}
for(int mid=1;mid<n;mid<<=1){
int tmp=pows(op==1?n1:n2,(mods-1)/(mid<<1));
for(int i=0;i<n;i+=mid<<1){
int om=1;
for(int j=0;j<mid;j++,om=om*tmp%mods){
int x=p[i+j],y=om*p[i+j+mid]%mods;
p[i+j]=(x+y)%mods,p[i+j+mid]=(x-y)%mods;
}
}
}
}
int gets(int x){
int res=1;
while((1<<res)<=x)res++;
return 1<<res;
}
void mul(int*a,int*b,int*g,int n){
int m=gets(2*n);
fft(a,m,1);fft(b,m,1);
for(int i=0;i<=m;i++)g[i]=a[i]*b[i]%mods;
fft(g,m,-1);
int ny=pows(m,mods-2);
for(int i=0;i<=m;i++)g[i]=g[i]*ny%mods;
}
int ff[1000005];
void inv(int*p,int*g,int n){
if(n==1){
g[0]=pows(p[0],mods-2);
return;
}
inv(p,g,n+1>>1);
int m=gets(2*n);
for(int i=0;i<=n;i++)ff[i]=p[i];
for(int i=n;i<=m;i++)ff[i]=0;
fft(g,m,1);fft(ff,m,1);
for(int i=0;i<=m;i++)g[i]=(2*g[i]-ff[i]*g[i]%mods*g[i]%mods)%mods;
fft(g,m,-1);
int nys=pows(m,mods-2);
for(int i=0;i<=n;i++)g[i]=g[i]*nys%mods;
for(int i=n;i<=m;i++)g[i]=0;
}
int C(int a,int b){
return jc[a]*ny[b]%mods*ny[a-b]%mods;
}
signed main(){
cin>>n;
jc[0]=ny[0]=1;
for(int i=1;i<=n;i++)jc[i]=jc[i-1]*i%mods,ny[i]=pows(jc[i],mods-2);
for(int i=1;i<=n;i++)f[i]=mods-pows(-1,i+1)*ny[i]%mods*pows(ny2,i*i)%mods;
f[0]=1;
inv(f,g,n+1);
for(int i=0;i<=n;i++)g[i]=g[i]*jc[i]%mods*pows(s2,i*i)%mods;
g[0]=1;
for(int i=0;i<=n;i++)gs[i]=g[i]*ny[i]%mods;
for(int i=1;i<=n;i++)g[i]=g[i]*ny[i-1]%mods;
g[0]=0;
inv(gs,fs,n+1);
mul(fs,g,rs,n);
for(int i=1;i<=n;i++){
rs[i]=rs[i]*jc[i-1]%mods;
cout<<(rs[i]+mods)%mods<<endl;
}
}
```
## [付公主的背包](https://www.luogu.com.cn/problem/P4389)
$f_{n,k}=\sum\limits_{i=0}f_{n-1,k-iV}$。
即 $f_n=f_{n-1}(1+x^{V_n}+x^{2V_n}+...)=f_{n-1}\sum\limits_{i=0}x^{iV_n}$。
答案就是 $\prod\limits_{i=1}^n\sum\limits_{j=0}x^{jV_i}$。
令 $S=\sum\limits_{j=0}x^{jV_i}
x^{V_i}S=\sum\limits_{j=1}x^{jV_i}
x^{V_i}S-S=-1
S(x^{V_i}-1)=-1
$ans=\prod\limits_{i=1}^n\frac{1}{1-x^{V_i}}
\ln ans=\ln\prod\limits_{i=1}^n\frac{1}{1-x^{V_i}}
=\sum\limits_{i=1}^n\ln\frac{1}{1-x^{V_i}}
=\sum\limits_{i=1}^n\ln(1-x^{V_i})^{-1}
=-\sum\limits_{i=1}^n\ln(1-x^{V_i})
设 f(x)=1-x^{V_i}
\ln' f(x)=\frac{f'(x)}{f(x)}
=\frac{-V_ix^{V_i-1}}{1-x^{V_i}}
=-V_ix^{V_i-1}\sum\limits_{j=0}x^{jV_i}
=-\sum\limits_{j=0}V_ix^{jV_i+V_i-1}
由于 ax^b 积分的结果为 \frac{a}{b+1}x^{b+1}
\ln f(x)=-\sum\limits_{j=0}\frac{V_i}{jV_i+V_i}x^{jV_i+V_i}
=-\sum\limits_{j=0}\frac{1}{j+1}x^{(j+1)V_i}
最后 exp 就行了
# 恒等式
$\frac{1}{1-x^k}=1+x^k+x^{2k}+...
\frac{1-x^{k+1}}{1-x}=1+x+x^2+...+x^k
\frac{1}{(1-x)^k}=\sum\limits_{i=0}C(i+k-1,k-1)x^i