ZJT's Blog

Keep your determination.

扩展拉格朗日反演

2018-3-4 9:432018-4-6 16:8
拉格朗日反演FFT生成函数

突然发现拉格朗日反演其实还有扩展的形式。。。

定义$F^{-1}$为$F$的复合逆:

$$ F^{-1}(F(x))=x $$

$$[x^n]G(F^{-1}(x))=\frac1n[x^{-1}] ( \frac{dG(x)}{dx}\frac1{F^n(x)})$$

(这个好像叫Lagrange–Bürmann formula

注意一般$G$的表达式都同时含有$x$和$F^{-1}(x)$,我们把$F(x)$代入$x$后就能得到$G(x)$的表达式,进而求出$\frac {dG(x)}{dx}$。

感觉还是蛮有用的(大雾)

查看详细内容

[BZOJ 3684] 大朋友与多叉树

2017-10-29 22:232017-10-30 7:16
生成函数FFT拉格朗日反演

传送门

一直觉得这题是某个不可做神题,结果一看发现是个求复合逆模板题。。。大概是因为拉格朗日反演本身没那么普及?

考虑关于根节点权值生成函数$S$,满足$S=\sum_k S^{d_k}+x$,直接把右边的前一坨移到左边,就变成了一个复合逆的形式,直接拉格朗日反演就做完了。

写起来也没想象中的那么难?开始写到AC花了不到1h,过了样例就1A了。大概就是考会不会拉格朗日反演吧。。。

#include <cstdio>
#include <cstring>
#include <algorithm>
#define MAXN 400010
#define MAXW 262144
#define LL long long
using namespace std;

const LL P=950009857,G=7;

int n,m,sizew;
LL w[MAXW],inv_v[MAXN];

LL getPow(LL x,LL y){
    LL res=1;
    while(y){
        if(y&1) res=res*x%P;
        x=x*x%P;
        y>>=1;
    }
    return res;
}

void FFT(LL *a,int len,int flag){
    static int rev[MAXN],revlen;
    if(revlen!=len){
        revlen=len;
        for(int i=1;i<len;i++) rev[i]=rev[i>>1]>>1|((i&1)?(len>>1):0);
    }
    for(int i=0;i<len;i++)
        if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(int l=2;l<=len;l<<=1){
        int l2=l>>1;
        for(int i=0;i<len;i+=l){
            for(int j=0;j<l2;j++){
                LL t1=a[i+j],t2=a[i+j+l2]*w[MAXW/l*j];
                a[i+j]=(t1+t2)%P;
                a[i+j+l2]=(t1-t2)%P;
            }
        }
    }
    if(flag==-1){
        LL invn=getPow(len,P-2);
        for(int i=0;i<len;i++) a[i]=a[i]*invn%P;
        for(int i=1;i<len;i++) if(i<len-i) swap(a[i],a[len-i]);
    }
}

void getMul(LL *b,LL *a,int len){
    static LL t1[MAXN],t2[MAXN];
    for(int i=0;i<len;i++){
        t1[i]=b[i]; t2[i]=a[i];
        t1[i+len]=t2[i+len]=0;
    }
    FFT(t1,len<<1,1); FFT(t2,len<<1,1);
    for(int i=0;i<len<<1;i++) t1[i]=t1[i]*t2[i]%P;
    FFT(t1,len<<1,-1);
    for(int i=0;i<len;i++) b[i]=t1[i];
}

void getInv(LL *b,LL *a,int len){
    if(len==1){ b[0]=getPow(a[0],P-2); return; }
    static LL t1[MAXN],t2[MAXN];
    getInv(t1,a,len>>1);
    for(int i=len>>1;i<len<<1;i++) t1[i]=0;
    for(int i=0;i<len;i++) t2[i]=a[i],t2[i+len]=0; 
    FFT(t1,len<<1,1); FFT(t2,len<<1,1);
    for(int i=0;i<len<<1;i++) t1[i]=(2*t1[i]-t2[i]*t1[i]%P*t1[i])%P;
    FFT(t1,len<<1,-1);
    for(int i=0;i<len;i++) b[i]=t1[i];
}

void getLn(LL *b,LL *a,int len){
    static LL t1[MAXN],t2[MAXN],t3[MAXN];
    getInv(t1,a,len);
    for(int i=1;i<len;i++) t2[i-1]=a[i]*i%P;
    for(int i=len-1;i<len<<1;i++) t2[i]=0;
    getMul(t2,t1,len);
    b[0]=0;
    for(int i=1;i<len;i++) b[i]=t2[i-1]*inv_v[i]%P;
}

void getExp(LL *b,LL *a,int len){
    if(len==1){ b[0]=1; return; }
    static LL t1[MAXN],t2[MAXN],t3[MAXN];
    getExp(t1,a,len>>1);
    for(int i=0;i<len>>1;i++) t3[i]=t1[i],t3[i+(len>>1)]=0;
    getLn(t2,t3,len);
    for(int i=0;i<len>>1;i++) t2[i]=t2[i+(len>>1)]-a[i+(len>>1)],t2[i+(len>>1)]=0;
    getMul(t2,t3,len);
    for(int i=0;i<len>>1;i++) b[i]=t3[i],b[i+(len>>1)]=-t2[i];
}

void getPow(LL *b,LL *a,int k,int len){
    static LL t1[MAXN];
    getLn(t1,a,len);
    for(int i=0;i<len;i++) t1[i]=t1[i]*k%P;
    getExp(b,t1,len);
}

void init(){
    for(sizew=1;sizew<=n;sizew<<=1);
    w[0]=1;
    w[1]=getPow(G,(P-1)/MAXW);
    for(int i=2;i<MAXW;i++) w[i]=w[i-1]*w[1]%P;
    inv_v[1]=1;
    for(int i=2;i<MAXN;i++) inv_v[i]=-(P/i)*inv_v[P%i]%P;
}

LL a[MAXN];

void test(){
    int len=32;
    static LL t1[MAXN],t2[MAXN];
    a[0]=1;
    for(int i=1;i<len;i++) a[i]=rand()%P;
    getLn(t1,a,len);
    getExp(t2,t1,len);
    for(int i=0;i<len;i++) t2[i]=(t2[i]+P)%P;
    return;
}

void gao(){
    a[0]=1;
    while(m--){
        int x;
        scanf("%d",&x);
        a[x-1]--;
    }
    static LL t1[MAXN],t2[MAXN];
    getPow(t1,a,n,sizew);
    getInv(t2,t1,sizew);
    LL ans=t2[n-1]*inv_v[n]%P;
    printf("%lld\n",(ans+P)%P);
}

int main(){
#ifdef DEBUG
    freopen("in","r",stdin);
#endif
    scanf("%d%d",&n,&m);
    init();
    gao();
    return 0;
}

查看详细内容

拉格朗日反演

2017-10-29 21:52018-5-7 20:17
生成函数FFT拉格朗日反演

很久以前(大概去年联赛前后?)看到hgr他们在做拉格朗日反演的题,然后就去学了一个鏼爷的论文,假装自己会了。结果前几天突然想起来有这么一个东西,发现一点都不记得了。。。然后就重新学一遍咯。

一些更基础的东西

一上来就发现式子里有一个取$x^{-1}$项的系数,然后就一脸懵逼了。。。我好好的多项式,怎么就有$x^{-1}$的项了呢。。。?后来才发现这个东西根本不是在形式幂级数环下面进行的,是在形式幂级数的分式环下进行的。。。然后抽代学的一坨屎的我自然也不知道分式环是个啥,然后就只能又去学了一些更加基础的东西。

正常来说,我们遇到的幂级数都是形如$a_0+a_1x+a_2x^2+...$这样的东西的,我们把所有这些幂级数(构成一个环)记为$F[[x]]$($F$是某个域,OI里面一般要么是复数域,要么是模质数$P$下的域),称为$F$下的形式幂级数环

然而这样的东西显然不能构成一个域,比如我$1/x$就不能表示成另一个形式幂级数。这时候,我们就需要一个域,使得这个域包含整个形式幂级数环,并且正常的多项式操作依然能用,最后这个域还得尽可能小。这个东西在抽代里面就叫做分式环,一个环$R$的分式环,是由所有$a/b(a,b\in R)$定义的(这个除法只是一个形式,事实上$b$在$R$中不一定真的可逆,比如多项式$x$在$F[[x]]$中就不可逆)。虽然它叫分式环,但它实际上是个域。我们把$F[[x]]$的分式环,记为$F((x))$。这样的话,只要在$F((x))$,我们就不用担心做不了多项式除法啦(除非除数就是零多项式)!

还有一个东西,就是任何一个$F((x))$里的$A(x)/B(x)$,都能表示成$...+a_{-2}x^{-2}+a_{-1}x^{-1}+a_0+a_1x+a_2x^2+...$的形式,也就是把原来的幂级数往$x$的负数次幂那边扩了一下。这是因为,只要$B(x)\ne0$,我们就可以把$B(x)$成$x^dB'(x)$,这里$B'(x)$是某个常数项非$0$的多项式(不是$B$的导数!)因为$B'(x)$的常数项非$0$,所以$1/B'(x)$就能表示成$F[[x]]$中的一个幂级数(在$x=0$处泰勒展开)。所以我们只要把$A(x)/B'(x)$算出来(这是一个正常的幂级数),然后把整个除以$x^d$(也就是所有项的次数减掉$d$),就得到了上面的那个形式。这也解释了为什么拉格朗日反演定理里面会出现取$-1$次项系数的操作,因为整个运算根本不是在$F[[x]]$下,而是在$F((x))$下进行的,所以会存在幂级数的负数次幂。

拉格朗日反演定理

拉格朗日反演是用来求幂级数的复合逆的。所谓复合逆,就是我们有两个常数项为零且一次项非零的幂级数$f(x),g(x)$,这两个多项式满足:$$g(f(x))=x$$

我们知道$g(x)$的具体系数,现在要求$f(x)$的系数(其实知道$f$求$g$也可以,后面就知道怎么弄了)

拉格朗日反演可以用来解决这样的问题,但是它有个比较局限的地方,就是它只能用来求$f(x)$的某项系数。扯了半天,拉格朗日定理具体就是长这样:$$[x^n]f(x)=\frac1n[x^{-1}]\frac1{g(x)^n}$$

其中$[x^n]f(x)$表示取$f(x)$的$n$次项系数,同理$[x^{-1}]$就是取$-1$次项的系数了(注意只有在$F((x))$下才能这么做哦)。为了更方便计算,我们可以把烦人的分式环扔掉,直接用形式幂级数的方式来表示这个定理:$$[x^n]f(x)=\frac1n[x^{n-1}]\frac1{g'(x)^n}$$

这里$g'(x)$跟上面一样,不是$g(x)$的导数,指的是$g(x)$去掉常数项后剩下的常数非零的幂级数,即$g'(x)=g(x)/x$。如果你搞明白了前面关于$F[[x]]$和$F((x))$那一坨东西,那么这条式子跟拉格朗日反演定理一般形式的等价性就很显然了。如果你没搞明白,也没关系,因为正常使用的时候只会用到这条式子。

到这里为止,所有的操作都变成了形式幂级数下的操作,我们再也不用管烦人的$F((x))$了。可以看到,求$f(x)$的一项,只需要用到模$x^{n}$意义下的多项式求幂和求逆,这两样都可以在$O(n\log n)$里完成,总复杂度就是$O(n\log n)$了。

在证明定理之前还有一些东西

还记得之前提到不仅能知道$g$求$f$,还能知道$f$求$g$吗?

我们有一个定理,就是所有不含常数项的形式幂级数,在复合运算下构成一个群(至于为什么看下面),单位元是$e(x)=x$。由于群里面的元素和自己逆元的乘法是满足交换律的,所以$f(g(x))=x$也就等价于$g(f(x))=x$了。也就是说,我们知道$f$求$g$,和知道$g$求$f$其实是一样的!只要把$g$和$f$反过来,就变得一样了。我们后面对定理的证明也要用到这个转化。

然而至于为啥构成群,我也不会证。。。本来想上wiki查一个,结果碰巧梯子崩了,国内又搜不到相关的东西。所以。。。有没有dalao教教我这个傻逼啊。。。?

5.7更新:这个只有在常数项为零且一次项非零的情况下才成立。。之前写的东西可能都是假的== 左右逆元的存在大概是可以递推构造出来。。?

一些微小的证明

根据上面的转化,我们的$g(f(x))=x$就变成了$f(g(x))=x$。这样的话就好做一点了。

我们记$f(x)=a_0+a_1x+a_2x^2+...$,我们最终的目标就是求某个$a_n$。这样,$f(g(x))=x$就可以写成$\sum a_kg^k(x)=x$了。

我们考虑把等式两边求导,就变成了:$$\sum_{k\geq1}ka_kg^{k-1}(x)g'(x)=1$$

注意这里的$g'(x)$表示的就是$g(x)$的导数了。我们把等式两边都除以$g^n(x)$,再取$-1$次项的系数,得到:$$[x^{-1}]\sum_{k\geq1}ka_kg^{k-n-1}(x)g'(x)=[x^{-1}]\frac1{g^n(x)}$$

考虑左边的每一项,如果$k\ne n$,则$g^{k-n-1}(x)g'(x)$能表示成$\frac1{k-n}(g^{k-n})'(x)$。由于任何一个$F((x))$中的函数求导之后$-1$次项都为0,所以左边就只剩下$k=n$的那一项了:$$[x^{-1}]na_n\frac{g'(x)}{g(x)}=[x^{-1}]\frac1{g^n(x)}$$

然后我们把$\frac{g'(x)}{g(x)}$展开(具体怎么展我就懒得写了qwq),发现$-1$次项系数一定是$1$!我们把这个代回去,再变一下形,就变成了:$$a_n=\frac1n[x^{-1}]\frac1{g^n(x)}$$

这正是我们要证的拉格朗日反演定理。

关于实现?

还是放在例题里面吧。。毕竟要用到求逆+ln+exp,还是比较长的。。把这篇东西弄得太长就不好了。。。

想看的话可以看下一篇博客

最后推荐一个比较靠谱的讲拉格朗日反演的博客(我就是看的这个)。

查看详细内容