[BZOJ3684] 大朋友和多叉树 [生成函数+拉格朗日反演]

传送门

如果两个多项式 F(x),G(x) 满足常数项均为 0，一次项均不为 0，并且 G(F(x)) = x，那么称 F(x)与 G(x) 互为复合逆

有结论 : 

 

 

---

关于本题，答案的生成函数为

  

另集合 S 的生成函数为 H(x)，那么有 

令 G(x) = x - H(x) , 那么 G(F(x)) = x

然后就可以多项式求逆+快速幂

#include<bits/stdc++.h>
#define N 400050
#define C 20
using namespace std;
typedef long long ll;
const int Mod = 950009857, G = 7, inv2 = (Mod+1) / 2;
ll add(ll a, ll b){ return (a + b) % Mod;}
ll mul(ll a, ll b){ return a * b % Mod;}
ll read(){
	ll x = 0; char ch = 0;
	while(!isdigit(ch)) ch = getchar();
	while(isdigit(ch)) x = add(mul(x, 10), ch-'0'), ch = getchar();
	return x;
}
ll power(ll a, ll b){ ll ans = 1;
	for(;b;b>>=1){ if(b&1) ans = mul(ans,a); a = mul(a, a);}
	return ans;
}
#define poly vector<ll> 
#define C 21
poly w[C + 1];
int n, m; ll inv[N];
int up, bit, rev[N];
void Init(int len){ up = 1, bit = 0;
	while(up < len) up <<= 1, bit++;
	for(int i = 0; i < up; i++) rev[i] = (rev[i>>1]>>1) | ((i&1) << (bit-1));
}
void prework(){
	for(int i = 1; i <= C; i++) w[i].resize(1 << (i-1));
	ll wn = power(G, (Mod - 1) / (1<<C));
	w[C][0] = 1;
	for(int i = 1; i < (1<<(C-1)); i++) w[C][i] = mul(w[C][i-1], wn);
	for(int i = C-1; i; i--)
		for(int j = 0; j < (1<<(i-1)); j++)
			w[i][j] = w[i+1][j << 1];
}
void NTT(poly &a, int flag){
	for(int i = 0; i < up; i++) if(i < rev[i]) swap(a[i], a[rev[i]]);
	for(int i = 1, l = 1; i < up; i <<= 1, l++)
		for(int j = 0; j < up; j += (i<<1))
			for(int k = 0; k < i; k++){ 
				int x = a[k + j], y = mul(w[l][k], a[k + j + i]);
				a[k + j] = add(x, y); a[k + j + i] = add(x, Mod - y);
			} 
	if(flag == -1){
		reverse(a.begin() + 1, a.begin() + up);
		for(int i = 0; i < up; i++) a[i] = mul(a[i], inv[up]);
	}
}
poly operator * (poly a, poly b){
	int len = a.size() + b.size() - 1;
	Init(len); a.resize(up); b.resize(up);
	NTT(a, 1); NTT(b, 1); for(int i = 0; i < up; i++) a[i] = mul(a[i], b[i]);
	NTT(a, -1); a.resize(len); return a;
}
poly deriv(poly a){
	for(int i = 0; i < a.size()-1; i++) a[i] = mul(a[i+1], i+1);
	a.pop_back(); return a;
}
poly integ(poly a){
	a.push_back(0);
	for(int i = a.size() - 1; i >= 1; i--) a[i] = mul(a[i-1], inv[i]);
	a[0] = 0; return a;
}
poly Inv(poly a, int len){
	poly b(1, power(a[0], Mod-2)), c;
	for(int lim = 4; lim < (len << 2); lim <<= 1){
		c = a; c.resize(lim >> 1);
		Init(lim); b.resize(up); c.resize(up); NTT(b, 1); NTT(c, 1);
		for(int i = 0; i < up; i++) b[i] = mul(b[i], add(2, Mod - mul(b[i], c[i])));
		NTT(b, -1); b.resize(lim >> 1);
	} b.resize(len); return b;
}
poly Ln(poly a, int len){ a = integ(deriv(a) * Inv(a, len)); a.resize(len); return a;} 
poly Exp(poly a, int len){
	poly b(1, 1), c; a.resize(len << 1);
	for(int lim = 2; lim < (len << 1); lim <<= 1){
		c = Ln(b, lim);
		for(int i = 0; i < lim; i++) c[i] = add(Mod - c[i], a[i]);
		c[0] = add(c[0], 1); b = b * c;
		b.resize(lim);
	} return b;
}
poly Ksm(poly a, int k, int len){
	poly b(len, 0); a = Ln(a, len); 
	for(int i = 0; i < len; i++) b[i] = mul(a[i], k);
	return Exp(b, len);
}
int main(){
	scanf("%d%d", &n, &m); poly a(n+1, 0);
	prework(); inv[0] = inv[1] = 1; 
	for(int i = 2; i <= N-50; i++) inv[i] = mul(Mod - Mod/i, inv[Mod%i]);
	a[0] = 1;
	for(int i = 0, x; i < n; i++) scanf("%d", &x), a[x-1] = Mod-1;
	poly b = Ksm(Inv(a, n), n, n);
	cout << mul(power(n, Mod - 2), b[n-1]);
	return 0; 
}