ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

拉格朗日反演学习及其应用

2022-06-24 15:35:06  阅读:256  来源: 互联网

标签:拉格朗 frac int poly 学习 反演 1ll limit mod


拉格朗日反演

  • 多项式复合:\(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);
}

标签:拉格朗,frac,int,poly,学习,反演,1ll,limit,mod
来源: https://www.cnblogs.com/Arashimu0x7f/p/16409071.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有