码迷,mamicode.com
首页 > 其他好文 > 详细

@总结 - 11@ 拉格朗日反演与复合逆

时间:2019-09-26 23:05:53      阅读:181      评论:0      收藏:0      [点我收藏+]

标签:algorithm   math   mp3   +=   deb   algo   list   printf   转换   


@0 - 参考资料@

参考博客-1
zjt 的博客地址挂掉了,所以这里暂时没有。

@1 - 问题描述@

我们知道,有时候可以使用牛顿迭代法求 G(F(x)) = 0 的某个根 F0(x)。

但是有时候函数 G 是一个非常复杂的多项式。
这个时候使用牛顿迭代法,每一层的复杂度可能要做 n 次乘法与加法(n 为当前多项式长度)。
于是时间复杂度将会急剧退化,甚至不如暴力。

但假如 G(F(x)) = x,其中 G, F 都为多项式函数,在已知 F 的情况下,我们可以使用拉格朗日反演求出 G 的某一项。

同时,假如 G(F(x)) = H(x),我们也可以使用扩展拉格朗日反演求出在已知 F 的情况下求 G 的某一项。

@2 - 拉格朗日反演@

@2.1 - 解法及证明@

形式化地给出问题:

若已知多项式函数 F(x),且复合函数 G(F(x)) 满足 G(F(x)) = x,求 G(x) 的第 n 项。

此时可以看作是 x 在复合函数 G·F 的作用下保持不变,所以我们称 G 为 F 的复合逆(复合运算下的逆元)。
不过对于 F,G 是有要求的:要求两个多项式常数项为 0 且一次项系数不为 0。不然复合逆是不存在的。

有一个结论:若 F 为 G 的复合逆,则 G 为 F 的复合逆。即 F(G(x)) = x 意味着 G(F(x)) = x。
所以:我们称 F 与 G 互为复合逆。
这个技巧可以在已知 F 的情况下,将 F(G(x)) = x 转为 G(F(x)) = x,有助于我们解决一些题目。

开始描述解法:
\(G(x) = \sum_{i=1}a_i*x^i\),那么就有:
\[\sum_{i=1}a_i*F^i(x) = x\]

两边求导,得到:
\[\sum_{i=1}i*a_i*F^{i-1}(x)*F'(x) = 1\]

为了得到 G 的第 n 项 \(a_n\),我们将两边除以 \(F^n(x)\),得到:
\[\sum_{i=1}i*a_i*F^{i-n-1}(x)*F'(x) = \frac{1}{F^n(x)}\]

(这个地方注意一点:我们此时函数并不一定为多项式——我们将函数推广到了负整数幂)
对于幂函数 \(F^a(x)\),当 a ≠ 0 时,它的导数为 \(a*F^{a-1}(x)*F'(x)\)
反过来,当 b ≠ -1 时,\(F^{b}(x)*F'(x)\) 的积分为 \(\frac{F^{b+1}(x)}{b+1}\)

于是对上面的式子进行进一步转换:
\[n*a_n*\frac{F'(x)}{F(x)}+\sum_{i=1}^{i\not = n}\frac{i}{i-n}*a_i*(F^{i-n}(x))' = \frac{1}{F^n(x)}\]

等式左右对应项系数相等。假如我们取 \(x^{-1}\) 的系数,则等式左边只留下 \(n*a_n*\frac{F'(x)}{F(x)}\) 这一部分(另一部分因为幂函数求导不会出现 \(x^{x-1}\) 这一项)。

于是就有:
\[[x^{-1}](n*a_n*\frac{F'(x)}{F(x)}) = [x^{-1}](\frac{1}{F^n(x)})\]

其中\([x^{-1}]\) 表示取 \(x^{-1}\) 的系数。

然后我们考虑一下 \(\frac{F'(x)}{F(x)}\)\(x^{-1}\) 的系数:
\[ \begin{align}\frac{F'(x)}{F(x)}&=\frac{a_1+2a_2x+3a_3x^2+\cdots}{a_1x+a_2x^2+\cdots}\&=\frac{a_1+2a_2x+3a_3x^2+\cdots}{a_1x}\cdot \frac{1}{1+\left(\frac{a_2}{a_1}x+\frac{a_3}{a_1}x^2+\cdots\right)}\&=\left(x^{-1}+\frac{2a_2}{a_1}+\cdots\right)\left(1-x\left(\frac{a_2}{a_1}+\frac{a_3}{a_1}x+\cdots\right)\right) \end{align} \]

总之经过一系列代数变形得到 \([x^{-1}]\frac{F'(x)}{F(x)} = 1\)
于是就有:
\[n*a_n =[x^{-1}](\frac{1}{F^n(x)})\]
\[a_n = \frac{1}{n}*[x^{-1}](\frac{1}{F^n(x)}) = \frac{1}{n}*[x^{n-1}](\frac{1}{(\frac{F(x)}{x})^n})\]

跑一个多项式逆元与多项式 exp + ln 求 k 次幂即可。

@2.2 - 典例@

bzoj3684:大朋友和多叉树
题目点这里查看

根据题目给出的定义,设权值为 i 的结点为根的树方案数为 {ai},设 f(x) 为 {ai} 的生成函数。
则有:
\[f(x) = \sum_{k\in D}f^k(x) + x\]

其中最后一个 x 表示叶子结点的方案数。
稍微变形就可以得到复合逆的基本形式 \(g(f(x)) = x\),直接复合逆即可。
代码:

#include<cstdio>
#include<algorithm>
using namespace std;
const int MOD = 950009857;
const int MAXN = 400000;
const int G = 5;
int pow_mod(int b, int p) {
    int ret = 1;
    while( p ) {
        if( p & 1 ) ret = 1LL*ret*b%MOD;
        b = 1LL*b*b%MOD;
        p >>= 1;
    }
    return ret;
}
struct poly{
    int pw[20 + 5], ipw[20 + 5];
    poly() {
        for(int i=0;i<=21;i++)
            pw[i] = pow_mod(G, (MOD-1)/(1<<i)), ipw[i] = pow_mod(pw[i], MOD-2);
    }
    void debug(int *A, int n) {
        for(int i=0;i<n;i++)
            printf("%d ", A[i]);
        puts(""), puts("");
    }
    void clear(int *A, int l, int r) {
        for(int i=l;i<r;i++)
            A[i] = 0;
    }
    void copy(int *A, int *B, int n) {
        for(int i=0;i<n;i++)
            A[i] = B[i];
    }
    int length(int n) {
        int len; for(len = 1; len < n; len <<= 1);
        return len;
    }
    void ntt(int *A, int n, int type) {
        for(int i=0,j=0;i<n;i++) {
            if( i < j ) swap(A[i], A[j]);
            for(int k=(n>>1);(j^=k)<k;k>>=1);
        }
        for(int i=1;(1<<i)<=n;i++) {
            int s = (1<<i), t = (s>>1);
            int u = (type == 1) ? pw[i] : ipw[i] ;
            for(int j=0;j<n;j+=s) {
                for(int k=0,p=1;k<t;k++,p=1LL*p*u%MOD) {
                    int x = A[j+k], y = 1LL*p*A[j+k+t]%MOD;
                    A[j+k] = (x + y)%MOD, A[j+k+t] = (x + MOD - y)%MOD;
                }
            }
        }
        if( type == -1 ) {
            int iv = pow_mod(n, MOD-2);
            for(int i=0;i<n;i++)
                A[i] = 1LL*A[i]*iv%MOD;
        }
    }
    int tmp1[MAXN + 5], tmp2[MAXN + 5];
    void poly_mul(int *A, int *B, int *C, int n, int m, int k) {
        int len = length(n + m - 1);
        clear(tmp1, 0, len), clear(tmp2, 0, len);
        for(int i=0;i<n;i++) tmp1[i] = A[i];
        for(int i=0;i<m;i++) tmp2[i] = B[i];
        ntt(tmp1, len, 1), ntt(tmp2, len, 1);
        for(int i=0;i<len;i++)
            C[i] = 1LL*tmp1[i]*tmp2[i]%MOD;
        ntt(C, len, -1);
    }
    int tmp3[MAXN + 5];
    void poly_inv(int *A, int *B, int n) {
        if( n == 1 ) {
            B[0] = pow_mod(A[0], MOD-2);
            return ;
        }
        poly_inv(A, B, (n+1)>>1);
        poly_mul(A, B, tmp3, n, (n+1)>>1, n);
        for(int i=0;i<n;i++)
            tmp3[i] = (MOD - tmp3[i])%MOD;
        tmp3[0] = (tmp3[0] + 2)%MOD;
        poly_mul(tmp3, B, B, n, (n+1)>>1, n);
    }
    void poly_int(int *A, int *B, int n) {
        for(int i=n-1;i>0;i--)
            B[i] = 1LL*pow_mod(i, MOD-2)*A[i-1]%MOD;
        B[0] = 0;
    }
    void poly_der(int *A, int *B, int n) {
        for(int i=0;i<n-1;i++)
            B[i] = 1LL*(i+1)*A[i+1]%MOD;
        B[n-1] = 0;
    }
    int tmp4[MAXN + 5];
    void poly_ln(int *A, int *B, int n) {
        poly_inv(A, tmp4, n);
        poly_der(A, B, n);
        poly_mul(B, tmp4, B, n, n, n);
        poly_int(B, B, n);
    }
    int tmp5[MAXN + 5];
    void poly_exp(int *A, int *B, int n) {
        if( n == 1 ) {
            B[0] = 1;
            return ;
        }
        poly_exp(A, B, (n+1)>>1), clear(B, (n+1)>>1, n);
        poly_ln(B, tmp5, n);
        for(int i=0;i<n;i++)
            tmp5[i] = (A[i] + MOD - tmp5[i])%MOD;
        tmp5[0] = (tmp5[0] + 1)%MOD;
        poly_mul(tmp5, B, B, n, (n+1)>>1, n);
    }
    int tmp6[MAXN + 5];
    void poly_pow(int *A, int *B, int n, int k) {
        poly_ln(A, tmp6, n);
        for(int i=0;i<n;i++)
            tmp6[i] = 1LL*k*tmp6[i]%MOD;
        poly_exp(tmp6, B, n);
    }
}op;
int s, m, f[MAXN + 5], g[MAXN + 5];
int main() {
    scanf("%d%d", &s, &m);
    for(int i=1;i<=m;i++) {
        int x; scanf("%d", &x);
        f[x] = (MOD - 1);
    }
    f[1] = 1;
    for(int i=0;i<s;i++)
        f[i] = f[i+1];
    op.poly_pow(f, f, s, s);
    op.poly_inv(f, g, s);
    printf("%lld\n", 1LL*g[s-1]*pow_mod(s, MOD-2)%MOD);
}

@3 - 扩展拉格朗日反演@

@3.1 - 解法及证明@

形式化地给出问题:

若已知多项式函数 F(x),H(x),且复合函数 G(F(x)) 满足 G(F(x)) = H(x),求 G(x) 的第 n 项。

基本上和拉格朗日反演的推导过程是一样的,最后得到:
\[n*a_n =[x^{-1}](\frac{H'(x)}{F^n(x)})\]
\[a_n = \frac{1}{n}*[x^{n-1}](\frac{H'(x)}{(\frac{F(x)}{x})^n})\]

不过这里有一个常用的变换形式:
\(x' = F(x)\),则 \(x = F^{-1}(x')\)。这意味着我们可以将复合函数内部的某个函数“置换”出来。

@3.2 - 典例@

咕。

@总结 - 11@ 拉格朗日反演与复合逆

标签:algorithm   math   mp3   +=   deb   algo   list   printf   转换   

原文地址:https://www.cnblogs.com/Tiw-Air-OAO/p/11233318.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!