拉格朗日反演学习及其应用
拉格朗日反演
- 多项式复合：\(F(G(x))=x\)，则称\(F(x)\)与\(G(x)\)互为复合逆
- 存在条件：\([x^0]F(x)=0\)，\([x^1]F(x)\ne 0\)
- 拉格朗日反演：
- \([x^n]G(x)=\frac{1}{n}[x^{-1}](\frac{1}{F(x)})^n\)
- 但由于\([x^0]F(x)=0\)无法求逆，所以更通用的是：\([x^n]G(x)=\frac{1}{n}[x^{n-1}](\frac{x}{F(x)})^n\)
应用
BZOJ 684. 大朋友和多叉树
Problem
求有\(s\)个叶子结点，且每个结点的儿子个数属于集合\(D\)的树的个数。答案对\(950009857\)取模
Solve
设树的生成函数为\(F(x)\)，\([x^n]F(x)\)表示叶子节点个数为\(n\)的个数
则\(F(x)=x+\sum_{i\in D}F^i\)，前面的\(x\)的系数为\(1\)表示只有它自己即叶子节点，后面是它还有儿子的形成的树的个数。
形式化地，我们移动项得，\(F(x)-\sum_{i\in D}F^i=x\)
令\(G(x)=x-\sum_{i\in D}x^i\)，则\(G(F(x))=x\)
利用拉格朗日反演得到，\([x^n]F(x)=\frac{1}{n}[x^{n-1}](\frac{x}{G(x)})^n\)
注意到\(\frac{x}{G(x)}=\frac{1}{\frac{G(x)}{x}}\)，令\(G(x)^{'}=\frac{G(x)}{x}=1-\sum_{i\in D}x^{i-1}=1-H(x)\)
则\([x^n]F(x)=\frac{1}{n}[x^{n-1}](\frac{x}{G(x)})^n=\frac{1}{n}[x^{n-1}](\frac{1}{1-H(x)})^n\)
多项式求逆和快速幂即可。处理时可以令\(H(x)^{'}=1-H(x)=1+\sum_{i\in D}(mod-1)x^i\)，把减法变成加法
Code
//#pragma GCC optimize(2)
#include <bits/stdc++.h>
#define pb push_back
using namespace std;
typedef long long ll;
const int N = 5e6+10;
const int mod = 950009857;//g=7
template <typename T>void read(T &x)
{
x = 0;
register int f = 1;
register char ch = getchar();
while(ch < '0' || ch > '9') {if(ch == '-')f = -1;ch = getchar();}
while(ch >= '0' && ch <= '9') {x = x * 10 + ch - '0';ch = getchar();}
x *= f;
}
ll qpow(ll a, ll b)
{
ll res = 1;
while(b) {
if(b & 1) res = 1ll * res * a % mod;
a = 1ll * a * a % mod;
b >>= 1;
}
return res;
}
namespace Poly
{
#define mul(x, y) (1ll * x * y >= mod ? 1ll * x * y % mod : 1ll * x * y)
#define minus(x, y) (1ll * x - y < 0 ? 1ll * x - y + mod : 1ll * x - y)
#define plus(x, y) (1ll * x + y >= mod ? 1ll * x + y - mod : 1ll * x + y)
#define ck(x) (x >= mod ? x - mod : x)//取模运算太慢了
typedef vector<ll> poly;
const int G = 7;//根据具体的模数而定，原根可不一定不一样！！！
//一般模数的原根为 2 3 5 7 10 6
const int inv_G = qpow(G, mod - 2);
int RR[N], deer[2][21][N], inv[N];
void init(const int t) {//预处理出来NTT里需要的w和wn，砍掉了一个log的时间
for(int p = 1; p <= t; ++ p) {
int buf1 = qpow(G, (mod - 1) / (1 << p));
int buf0 = qpow(inv_G, (mod - 1) / (1 << p));
deer[0][p][0] = deer[1][p][0] = 1;
for(int i = 1; i < (1 << p); ++ i) {
deer[0][p][i] = 1ll * deer[0][p][i - 1] * buf0 % mod;//逆
deer[1][p][i] = 1ll * deer[1][p][i - 1] * buf1 % mod;
}
}
inv[1] = 1;
for(int i = 2; i <= (1 << t); ++ i)
inv[i] = 1ll * inv[mod % i] * (mod - mod / i) % mod;
}
inline int NTT_init(int n) {//快速数论变换预处理
int limit = 1, L = 0;
while(limit <= n) limit <<= 1, L ++ ;
for(int i = 0; i < limit; ++ i)
RR[i] = (RR[i >> 1] >> 1) | ((i & 1) << (L - 1));
return limit;
}
inline void NTT(poly &A, int type, int limit) {//快速数论变换
A.resize(limit);
for(int i = 0; i < limit; ++ i)
if(i < RR[i])
swap(A[i], A[RR[i]]);
for(int mid = 2, j = 1; mid <= limit; mid <<= 1, ++ j) {
int len = mid >> 1;
for(int pos = 0; pos < limit; pos += mid) {
int *wn = deer[type][j];
for(int i = pos; i < pos + len; ++ i, ++ wn) {
int tmp = 1ll * (*wn) * A[i + len] % mod;
A[i + len] = ck(A[i] - tmp + mod);
A[i] = ck(A[i] + tmp);
}
}
}
if(type == 0) {
for(int i = 0; i < limit; ++ i)
A[i] = 1ll * A[i] * inv[limit] % mod;
}
}
inline poly poly_mul(poly A, poly B) {//多项式乘法
int deg = A.size() + B.size() - 1;
int limit = NTT_init(deg);
poly C(limit);
NTT(A, 1, limit);
NTT(B, 1, limit);
for(int i = 0; i < limit; ++ i)
C[i] = 1ll * A[i] * B[i] % mod;
NTT(C, 0, limit);
C.resize(deg);
return C;
}
inline poly poly_inv(poly &f, int deg) {//多项式求逆
if(deg == 1)
return poly(1, qpow(f[0], mod - 2));
poly A(f.begin(), f.begin() + deg);
poly B = poly_inv(f, (deg + 1) >> 1);
int limit = NTT_init(deg << 1);
NTT(A, 1, limit), NTT(B, 1, limit);
for(int i = 0; i < limit; ++ i)
A[i] = B[i] * (2 - 1ll * A[i] * B[i] % mod + mod) % mod;
NTT(A, 0, limit);
A.resize(deg);
return A;
}
inline poly poly_exp(poly &f, int deg) {//多项式求指数
if(deg == 1)
return poly(1, 1);
poly B = poly_exp(f, (deg + 1) >> 1);
B.resize(deg);
poly lnB = poly_ln(B, deg);
for(int i = 0; i < deg; ++ i)
lnB[i] = ck(f[i] - lnB[i] + mod);
int limit = NTT_init(deg << 1);//n -> n^2
NTT(B, 1, limit), NTT(lnB, 1, limit);
for(int i = 0; i < limit; ++ i)
B[i] = 1ll * B[i] * (1 + lnB[i]) % mod;
NTT(B, 0, limit);
B.resize(deg);
return B;
}
poly poly_pow(poly f, int k) {//多项式快速幂
f = poly_ln(f, f.size());
for(auto &x : f) x = 1ll * x * k % mod;
return poly_exp(f, f.size());
}
//多项式牛顿迭代:求g(f(x))=0mod(x^n)中的f(x)
}
int main()
{
ios::sync_with_stdio(false);
cin.tie(nullptr);
Poly::init(20);//2^21 = 2,097,152,根据题目数据多项式项数的大小自由调整，注意大小需要跟deer数组同步(21+1=22)
int s,m;
read(s);
read(m);
Poly::poly H(s+2);
for(int i=1;i<=m;i++){
int x;
read(x);
H[x-1]=mod-1;
}
H[0]=1;
Poly::poly G=Poly::poly_inv(H,s+1);
G=Poly::poly_pow(G,s);
ll res=qpow(s,mod-2)*G[s-1]%mod;
printf("%lld\n",res);
return 0;
}
P3978 [TJOI2015]概率论
Problem
为了提高智商，ZJY 开始学习概率论。有一天，她想到了这样一个问题：对于一棵随机生成的 \(n\) 个结点的有根二叉树（所有互相不同构的形态等概率出现），它的叶子节点数的期望是多少呢？
Solve
这里的做法是对EI大佬做法自己的一点理解和解释
首先定义\(F(x,t)=\sum a_{n,m}x^nt^m\)，其中\(a_{n,m}\)表示一个具有\(n\)个结点和\(m\)个叶子的树的个数
那么，对于一个二叉树，一个结点可以自己成为一个叶子结点，即\(xt\)，也可以只有左子树\(F\)，也可以只有右子树\(T\)，两个都有的方案数\(F*F\)，那么就有\(F=x(t+2F+F*F)\)，把\(F\)看做是\(z\)的多项式，那么\(x=\frac{F}{t+2F+F*F}\)，令\(G(x)=\frac{x}{t+2x+x^2}\)，带入得\(G(F(x))=x\)，利用拉格朗日反演得到\([x^n]F(x)=\frac{1}{n}[x^{n-1}](\frac{x}{G(x)})^n=\frac{1}{n}[x^{n-1}](t+2x+x^2)^n\)。
由于求期望，那么我们就是求\(\frac{所有树的叶子节点之和}{所有树的个数}\)，所有树的个数显然可以用卡特兰数来计算，即\(H_n=\frac{\binom{2n}{n}}{n+1}\)，现在考虑怎么求所有树的叶子节点之和，对于一个具有\(n\)个结点和\(m\)个叶子的树，总的叶子数是\(a_{n,m}m\)，怎样才可以得到这个系数，方法是求导，即\(\frac{d}{dt}[x^n]F(x)=\frac{1}{n}[x^{n-1}]n(t+2x+x^2)^{n-1}\)，然后令\(t=1\)，得到\(\frac{d}{dt}[x^n]F(x)=[x^{n-1}](x+1)^{2(n-1)}=\binom{2n-2}{n-1}\)
所以答案就是\(\frac{\binom{2n-2}{n-1}}{\frac{\binom{2n}{n}}{n+1}}=\frac{n(n+1)}{2(n-1)}\)
Code
#include <bits/stdc++.h>
#define ll long long;
using namespace std;
int main()
{
ios::sync_with_stdio(false);
cin.tie(nullptr);
double n;
cin>>n;
double res=n*(n+1)/(2*(2*n-1));
printf("%.9lf\n",res);
}
P2767 树的数量
Problem
求\(n\)个点的\(m\)叉树的数量
Solve
令\(F(x)\)为\(m\)叉树的生成函数，那么\([x^n]F(x)\)为\(n\)个点的\(m\)叉树的数量
则\(F=x+x\sum_{i=1}^m\binom{m}{i}F^i=x(1+F)^m\)，意思就是一个结点是叶子节点用\(x\)来表示，也可以有分支，可以有\(1,\cdots,m\)个分支，用组合数算一下即可，但注意由于下面计算的都是当前结点的子树的个数，由于\(F\)是从\(x^0\)开始的，所以要乘上一个\(x\)代表目前已经选了一个父亲节点了。
那么\(x=\frac{F}{(1+F)^m}\)，令\(G(x)=\frac{x}{(1+x)^m}\)，有\(G(F(x))=x\)
由拉格朗日反演得到：\([x^n]F(x)=\frac{1}{n}[x^{n-1}](1+x)^{mn}=\frac{1}{n}\binom{nm}{n-1}\)
Code
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int mod=10007;
ll power(ll x,int y){
ll res=1;
while(y){
if(y&1) res=res*x%mod;
x=x*x%mod;
y>>=1;
}
return res;
}
int main()
{
ios::sync_with_stdio(false);
cin.tie(nullptr);
int n,m;
cin>>n>>m;
vector<ll>fac(m*n+1),infac(m*n+1);
fac[0]=1;
for(int i=1;i<=m*n;i++) fac[i]=fac[i-1]*i%mod;
infac[n*m]=power(fac[n*m],mod-2);
for(int i=n*m;i>=1;i--) infac[i-1]=infac[i]*i%mod;
ll res=fac[n*m]*infac[n-1]%mod*infac[n*m-n+1]%mod*power(n,mod-2)%mod;
cout<<res<<'\n';
}
浙公网安备 33010602011771号