隐藏
「bzoj3684」大朋友和多叉树 - 生成函数+多项式求导+拉格朗日反演+多项式k次幂(多项式lnp+多项式exp) | Bill Yang's Blog

路终会有尽头,但视野总能看到更远的地方。

0%

「bzoj3684」大朋友和多叉树 - 生成函数+多项式求导+拉格朗日反演+多项式k次幂(多项式lnp+多项式exp)

题目大意

    我们的大朋友很喜欢计算机科学,而且尤其喜欢多叉树。对于一棵带有正整数点权的有根多叉树,如果它满足这样的性质,我们的大朋友就会将其称作神犇的:点权为$1$的结点是叶子结点;对于任一点权大于$1$的结点$u$,$u$的孩子数目$deg[u]$属于集合$D$,且$u$的点权等于这些孩子结点的点权之和。
    给出一个整数$s$,你能求出根节点权值为$s$的神犇多叉树的个数吗?请参照样例以更好的理解什么样的两棵多叉树会被视为不同的。
    我们只需要知道答案关于$950009857$($453*2^{21}+1$,一个质数)取模后的值。


题目分析

设$f_i$表示根结点权值为$i$的神犇多叉树个数,$F(x)$为$f_i$的生成函数,$C(x)$
为$S$的生成函数,则有:

令$G(F(x))=x$,因此$F(x)$是$G(x)$的复合逆。
根据拉格朗日反演,可知:

其中$[x^n]F(x)$表示$F(x)$的第$n$次项系数。

然后就可以使用多项式求逆与多项式$k$次幂计算了。

多项式$k$次幂使用$lnp+exp$转化后比快速幂快。


代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
#include<algorithm>
#include<iostream>
#include<iomanip>
#include<cstring>
#include<cstdlib>
#include<climits>
#include<vector>
#include<cstdio>
#include<cmath>
#include<queue>
using namespace std;

inline const int Get_Int() {
int num=0,bj=1;
char x=getchar();
while(x<'0'||x>'9') {
if(x=='-')bj=-1;
x=getchar();
}
while(x>='0'&&x<='9') {
num=num*10+x-'0';
x=getchar();
}
return num*bj;
}

const int maxn=262144+5;
const int mod=950009857,g=7;

void check(int &x) {
if(x>=mod)x-=mod;
if(x<0)x+=mod;
}

void add(int &x,int v) {
x+=v;
check(x);
}

int Quick_Pow(int a,int b) {
int sum=1;
for(; b; b>>=1,a=1ll*a*a%mod)if(b&1)sum=1ll*sum*a%mod;
return sum;
}

int get_inv(int x) {
return Quick_Pow(x,mod-2);
}

int inv[maxn];

struct NumberTheoreticTransform {
int n,rev[maxn];
int omega[maxn],iomega[maxn];

void init(int n) {
this->n=n;
int x=Quick_Pow(g,(mod-1)/n);
omega[0]=iomega[0]=1;
for(int i=1; i<n; i++) {
omega[i]=1ll*omega[i-1]*x%mod;
iomega[i]=get_inv(omega[i]);
}
int k=log2(n);
for(int i=0; i<n; i++) {
int t=0;
for(int j=0; j<k; j++)if(i&(1<<j))t|=(1<<(k-j-1));
rev[i]=t;
}
}

void transform(int* a,int* omega) {
for(int i=0; i<n; i++)if(i<rev[i])swap(a[i],a[rev[i]]);
for(int len=2; len<=n; len*=2) {
int mid=len>>1;
for(int* p=a; p!=a+n; p+=len)
for(int i=0; i<mid; i++) {
int t=1ll*omega[n/len*i]*p[mid+i]%mod;
p[mid+i]=p[i]-t,check(p[mid+i]);
add(p[i],t);
}
}
}

void dft(int* a) {
transform(a,omega);
}

void idft(int* a) {
transform(a,iomega);
int x=inv[n];
for(int i=0; i<n; i++)a[i]=1ll*a[i]*x%mod;
}
} ntt;

void polynomial_inverse(const int* a,const int n,int* b) {
if(n==1) {
b[0]=get_inv(a[0]);
return;
}
polynomial_inverse(a,n>>1,b);
int p=n<<1;
static int x[maxn];
copy(a,a+n,x),fill(x+n,x+p,0);
ntt.init(p),ntt.dft(x),ntt.dft(b);
for(int i=0; i<p; i++)b[i]=1ll*b[i]*((2-1ll*x[i]*b[i]%mod+mod)%mod)%mod;
ntt.idft(b),fill(b+n,b+p,0);
}

void Multiply(const int* a1,const int n1,const int* a2,const int n2,int* ans) {
int n=1;
while(n<n1+n2)n<<=1;
static int c1[maxn],c2[maxn];
copy(a1,a1+n1,c1),fill(c1+n1,c1+n,0);
copy(a2,a2+n2,c2),fill(c2+n2,c2+n,0);
ntt.dft(c1);
ntt.dft(c2);
for(int i=0; i<n; i++)c1[i]=1ll*c1[i]*c2[i]%mod;
ntt.idft(c1);
for(int i=0; i<n1+n2-1; i++)ans[i]=c1[i];
}

void polynomial_lnp(const int* a,int n,int* b) {
static int tmp[maxn];
polynomial_inverse(a,n,tmp);
for(int i=0; i<n-1; i++)b[i]=1ll*a[i+1]*(i+1)%mod;
b[n-1]=0;
Multiply(b,n,tmp,n,b);
fill(b+n,b+(n<<1),0),fill(tmp,tmp+n,0);
for(int i=n-1; i>=1; i--)b[i]=1ll*b[i-1]*inv[i]%mod;
b[0]=0;
}

void polynomial_exp(const int *a,int n,int* b) {
if(n==1) {
b[0]=1;
return;
}
polynomial_exp(a,n>>1,b);
int p=n<<1;
static int x[maxn];
polynomial_lnp(b,n,x);
for(int i=0; i<n; i++)x[i]=-x[i]+a[i],check(x[i]);
add(x[0],1);
Multiply(b,n,x,n,b),fill(b+n,b+p,0);
}

int n,m;
int G[maxn],inv_G[maxn],ln_inv_G[maxn],inv_G_n[maxn];

int main() {
n=Get_Int();
m=Get_Int();
G[0]=1;
for(int i=1; i<=m; i++)G[Get_Int()-1]=mod-1;
int p=1;
while(p<n)p<<=1;
inv[1]=1;
for(int i=2; i<=(p<<1); i++)inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
ntt.init(p<<1);
polynomial_inverse(G,p,inv_G);
fill(inv_G+n,inv_G+p,0);
polynomial_lnp(inv_G,p,ln_inv_G);
fill(ln_inv_G+n,ln_inv_G+p,0);
for(int i=0; i<p; i++)ln_inv_G[i]=1ll*ln_inv_G[i]*n%mod;
polynomial_exp(ln_inv_G,p,inv_G_n);
printf("%d\n",1ll*inv_G_n[n-1]*inv[n]%mod);
return 0;
}
姥爷们赏瓶冰阔落吧~