zoukankan      html  css  js  c++  java
  • 斯特林数学习笔记

    斯特林数学习笔记


    好像全忘了,那就再开一个博客吧。。

    第一类斯特林数

    (egin{bmatrix} n\mend{bmatrix})表示将(n)个不同元素构成(m)个圆排列的方案数。

    递推式为(egin{bmatrix} n\mend{bmatrix}=egin{bmatrix} n-1\m-1end{bmatrix}+egin{bmatrix} n-1\mend{bmatrix}(n-1))

    (egin{bmatrix} n-1\m-1end{bmatrix})就是第(n)个元素新开一个圆排列,否则将第(n)个元素插进前面任意一个元素的后面。

    有两个式子:

    (x^{overline n}=sum_{i=0}^nx^iegin{bmatrix} n\iend{bmatrix})

    (x^{underline n}=sum_{i=0}^n(-1)^{n-i}x^iegin{bmatrix} n\iend{bmatrix})

    求一行

    https://www.luogu.org/problemnew/show/P5408

    套上面的上升幂式子,可以直接分治ntt,(O(nlog^2n))

    还有一种类似快速幂的做法:

    用类似快速幂来推这玩意儿,要支持快速做两种操作:

    (x^{overline n} ightarrow x^{overline{n+1}})
    (x^{overline n} ightarrow x^{overline{2n}})

    第一个直接暴力乘上去就好了。

    第二个列一下式子:

    (x^{overline{2n}}=x^{overline{n}}(x+n)^{overline{n}}),那么要求((x+n)^{overline{n}})

    套上面的式子:((x+n)^{overline{n}}=sum_{i=0}^n(x+n)^iegin{bmatrix} n\iend{bmatrix})

    ((x+n)^{overline{n}}=sum_{i=0}^negin{bmatrix} n\iend{bmatrix}sum_{j=0}^ix^jn^{i-j}C_{i}^{j})

    化成NTT可以接受的形式:

    (=sum_{j=0}^nfrac{x^j}{n^jj!}sum_{i=j}^negin{bmatrix} n\iend{bmatrix}n^ii! cdot frac{1}{(i-j)!})

    因为你已经求出来了(n)以内的斯特林数,所以可以求出((x+n)^{overline{n}}),然后乘起来可以得到(x^{overline{2n}})

    #include<bits/stdc++.h>
    #define il inline
    #define vd void
    #define mod 167772161
    #define poly std::vector<int>
    typedef long long ll;
    il ll gi(){
        ll x=0,f=1;
        char ch=getchar();
        while(!isdigit(ch))f^=ch=='-',ch=getchar();
        while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
        return f?x:-x;
    }
    il int pow(int x,int y){
        int ret=1;
        while(y){
            if(y&1)ret=1ll*ret*x%mod;
            x=1ll*x*x%mod;y>>=1;
        }
        return ret;
    }
    #define maxn 524289
    poly pA,pB;
    int rev[maxn],_lstN,P[maxn],iP[maxn];
    il vd ntt(int*A,int N,int t){
        for(int i=0;i<N;++i)if(rev[i]>i)std::swap(A[i],A[rev[i]]);
        for(int o=1;o<N;o<<=1){
            int W=t?P[o]:iP[o];
            for(int*p=A;p!=A+N;p+=o<<1)
                for(int i=0,w=1;i<o;++i,w=1ll*w*W%mod){
                    int t=1ll*w*p[i+o]%mod;
                    p[i+o]=(p[i]-t+mod)%mod;p[i]=(p[i]+t)%mod;
                }
        }
        if(!t){
            int inv=pow(N,mod-2);
            for(int i=0;i<N;++i)A[i]=1ll*A[i]*inv%mod;
        }
    }
    int N,lg;
    il vd setN(int n){
        N=1,lg=0;
        while(N<n)N<<=1,++lg;
        if(N!=_lstN)for(int i=0;i<N;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<lg-1);
    }
    il vd ntt(poly&a,int t){
        static int A[maxn];
        for(int i=0;i<a.size();++i)A[i]=a[i];memset(A+a.size(),0,4*(N-a.size()));
        ntt(A,N,t);
        a.resize(N);
        for(int i=0;i<N;++i)a[i]=A[i];
        int s=a.size();while(s&&!a[s-1])--s;
        a.resize(s);
    }
    il poly mul(poly a,poly b,int newn=-1){
        if(newn==-1)newn=a.size()+b.size()-1;
        setN(a.size()+b.size()-1);
        ntt(a,1),ntt(b,1);
        for(int i=0;i<N;++i)a[i]=1ll*a[i]*b[i]%mod;
        ntt(a,0);a.resize(newn);
        return a;
    }
    il poly operator+(poly a,const poly&b){
        if(a.size()<b.size())a.resize(b.size());
        for(int i=0;i<a.size();++i)if(i<b.size())a[i]=(a[i]+b[i])%mod;
        return a;
    }
    il poly operator-(poly a,const poly&b){
        if(a.size()<b.size())a.resize(b.size());
        for(int i=0;i<a.size();++i)if(i<b.size())a[i]=(a[i]-b[i]+mod)%mod;
        return a;
    }
    il poly operator*(poly a,int b){
        for(auto&i:a)i=1ll*i*b%mod;
        return a;
    }
    il poly qiudao(poly a){
        for(int i=0;i<a.size()-1;++i)a[i]=1ll*a[i+1]*(i+1)%mod;
        a.erase(a.end()-1);
        return a;
    }
    il poly jifen(poly a){
        a.insert(a.begin(),0);
        for(int i=1;i<a.size();++i)a[i]=1ll*a[i]*pow(i,mod-2)%mod;
        return a;
    }
    il poly getinv(poly a){
        if(a.size()==1)return poly(1,pow(a[0],mod-2));
        int n=a.size(),m=a.size()+1>>1;
        poly _a(m);
        for(int i=0;i<m;++i)_a[i]=a[i];
        poly b=getinv(_a);
        setN(n+m*2-2);
        ntt(a,1);ntt(b,1);
        for(int i=0;i<N;++i)a[i]=1ll*a[i]*b[i]%mod*b[i]%mod;
        ntt(a,0),ntt(b,0);
        a.resize(n);
        return b*2-a;
    }
    il poly getln(poly a,int n=-1){
        if(n==-1)n=a.size();
        a.resize(n);
        return jifen(mul(qiudao(a),getinv(a),n));
    }
    il poly getexp(poly a){
        if(a.size()==1)return a[0]=1,a;
        int n=a.size(),m=a.size()+1>>1;
        poly _a(m);
        for(int i=0;i<m;++i)_a[i]=a[i];
        poly b=getexp(_a);
        return mul(b,poly(1,1)-getln(b,a.size())+a,a.size());
    }
    il poly operator^(poly a,int b){
        int n=a.size();
        a=getexp(getln(a)*b);a.resize(n);
        return a;
    }
    il poly sqrt(poly a){
        if(a.size()==1)return a;
        int n=a.size(),m=a.size()+1>>1;
        poly _a(m);
        for(int i=0;i<m;++i)_a[i]=a[i];
        poly b=sqrt(_a);b.resize(n);
        return (b+mul(a,getinv(b),n))*(mod+1>>1);
    }
    il vd poly_init(){
        int G=3,iG=pow(G,mod-2);
        for(int i=1;i<maxn;i<<=1)P[i]=pow(G,(mod-1)/(i<<1)),iP[i]=pow(iG,(mod-1)/(i<<1));
    }
    int fact[262147],ifact[262147];
    il poly get_stelin(int n){
        if(n==1){poly r(2);r[1]=1;return r;}
        int m=n>>1;
        poly g=get_stelin(m),_g(g),_h(m+1),h(m+1);
        for(int i=0,pm=1;i<=m;++i,pm=1ll*pm*m%mod)_g[i]=1ll*_g[i]*fact[i]%mod*pm%mod;
        for(int i=0;i<=m;++i)_h[i]=ifact[m-i];
        _h=mul(_h,_g);
        int invm=pow(m,mod-2);
        for(int i=0,pm=1;i<=m;++i,pm=1ll*pm*invm%mod)h[i]=1ll*_h[m+i]*pm%mod*ifact[i]%mod;
        g=mul(h,g,n+1);
        if(n&1){
            g.resize(n+1);
            for(int i=n;i;--i)g[i]=(g[i-1]+1ll*(n-1)*g[i])%mod;
            g[0]=1ll*g[0]*(n-1)%mod;
        }
        return g;
    }
    int main(){
    #ifdef XZZSB
        freopen("in.in","r",stdin);
        freopen("out.out","w",stdout);
    #endif
        poly_init();
        int n=gi();
        fact[0]=1;for(int i=1;i<=n;++i)fact[i]=1ll*fact[i-1]*i%mod;
        ifact[n]=pow(fact[n],mod-2);for(int i=n-1;~i;--i)ifact[i]=1ll*ifact[i+1]*(i+1)%mod;
        poly f=get_stelin(n);
        for(int i=0;i<=n;++i)printf("%d ",f[i]);
        return 0;
    }
    

    求一列

    https://www.luogu.org/problemnew/show/P5409

    这东西太神仙了

    考虑用两种方式拆开((x+1)^t)

    第一种就是快速幂的方法,((x+1)^t=exp(tln(x+1))=sum_{i=0}^{infty}frac{t^iln(x+1)^i}{i!})

    第二种是拆二项式,((x+1)^t=sum_{i=0}^{infty}x^iC_{t}^i)

    ((x+1)^t=sum_{i=0}^{infty}frac{x^i}{i!}t^{underline i})

    拆下降幂

    ((x+1)^t=sum_{i=0}^{infty}frac{x^i}{i!}sum_{j=0}^{infty}(-1)^{i-j}egin{bmatrix}i\jend{bmatrix}t^j)

    ((x+1)^t=sum_{j=0}^{infty}t^jsum_{i=0}^{infty}frac{x^i}{i!}(-1)^{i-j}egin{bmatrix}i\jend{bmatrix})

    将两种拆法放一起

    (sum_{i=0}^{infty}frac{t^iln(x+1)^i}{i!}=sum_{i=0}^{infty}t^isum_{j=0}^{infty}frac{x^j}{j!}(-1)^{j-i}egin{bmatrix}j\iend{bmatrix})

    (frac{t^kln(x+1)^k}{k!}=sum_{j=0}^{infty}frac{x^j}{j!}(-1)^{j-k}egin{bmatrix}j\iend{bmatrix})

    因为(t)是变量,所以任意取一个(i)两边都是相等的,(i)(k)就可以求得第(k)列的生成函数,(mod x^{n+1})可以求得0到n的(egin{bmatrix}i\kend{bmatrix})

    #pragma GCC optimize("O3,Ofast,no-stack-protector,unroll-loops,fast-math")
    #include<bits/stdc++.h>
    #define il inline
    #define vd void
    typedef long long ll;
    il ll gi(){
        ll x=0,f=1;
        char ch=getchar();
        while(!isdigit(ch))f^=ch=='-',ch=getchar();
        while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
        return f?x:-x;
    }
    //================多项式板子部分===================
    #define mod 167772161
    #define maxn 524289
    #define Gmod 3
    #define poly std::vector<int>
    il int pow(int x,int y){
        int ret=1;
        while(y){
            if(y&1)ret=1ll*ret*x%mod;
            x=1ll*x*x%mod;y>>=1;
        }
        return ret;
    }
    int rev[maxn],_lstN,P[maxn],iP[maxn];
    il vd ntt(int*A,int N,int t){
        for(int i=0;i<N;++i)if(rev[i]>i)std::swap(A[i],A[rev[i]]);
        for(int o=1;o<N;o<<=1){
            int W=t?P[o]:iP[o];
            for(int*p=A;p!=A+N;p+=o<<1)
                for(int i=0,w=1;i<o;++i,w=1ll*w*W%mod){
                    int t=1ll*w*p[i+o]%mod;
                    p[i+o]=(p[i]-t+mod)%mod;p[i]=(p[i]+t)%mod;
                }
        }
        if(!t){
            int inv=pow(N,mod-2);
            for(int i=0;i<N;++i)A[i]=1ll*A[i]*inv%mod;
        }
    }
    int N,lg;
    il vd setN(int n){
        N=1,lg=0;
        while(N<n)N<<=1,++lg;
        if(N!=_lstN){_lstN=N;for(int i=0;i<N;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<lg-1);}
    }
    il vd ntt(poly&a,int t){
        static int A[maxn];
        for(int i=0;i<a.size();++i)A[i]=a[i];memset(A+a.size(),0,4*(N-a.size()));
        ntt(A,N,t);
        a.resize(N);
        for(int i=0;i<N;++i)a[i]=A[i];
        int s=a.size();while(s&&!a[s-1])--s;
        a.resize(s);
    }
    il poly mul(poly a,poly b,int newn=-1){
        if(newn==-1)newn=a.size()+b.size()-1;
        setN(a.size()+b.size()-1);
        ntt(a,1),ntt(b,1);
        for(int i=0;i<N;++i)a[i]=1ll*a[i]*b[i]%mod;
        ntt(a,0);a.resize(newn);
        return a;
    }
    il poly operator+(poly a,const poly&b){
        if(a.size()<b.size())a.resize(b.size());
        for(int i=0;i<a.size();++i)if(i<b.size())a[i]=(a[i]+b[i])%mod;
        return a;
    }
    il poly operator-(poly a,const poly&b){
        if(a.size()<b.size())a.resize(b.size());
        for(int i=0;i<a.size();++i)if(i<b.size())a[i]=(a[i]-b[i]+mod)%mod;
        return a;
    }
    il poly operator*(poly a,int b){
        for(auto&i:a)i=1ll*i*b%mod;
        return a;
    }
    il poly qiudao(poly a){
        for(int i=0;i<a.size()-1;++i)a[i]=1ll*a[i+1]*(i+1)%mod;
        a.erase(a.end()-1);
        return a;
    }
    il poly jifen(poly a){
        a.insert(a.begin(),0);
        for(int i=1;i<a.size();++i)a[i]=1ll*a[i]*pow(i,mod-2)%mod;
        return a;
    }
    il poly getinv(poly a){
        if(a.size()==1)return poly(1,pow(a[0],mod-2));
        int n=a.size(),m=a.size()+1>>1;
        poly _a(m);
        for(int i=0;i<m;++i)_a[i]=a[i];
        poly b=getinv(_a);
        setN(n+m*2-2);
        ntt(a,1);ntt(b,1);
        for(int i=0;i<N;++i)a[i]=1ll*a[i]*b[i]%mod*b[i]%mod;
        ntt(a,0),ntt(b,0);
        a.resize(n);
        return b*2-a;
    }
    il poly getln(poly a,int n=-1){
        if(n==-1)n=a.size();
        a.resize(n);
        return jifen(mul(qiudao(a),getinv(a),n));
    }
    il poly getexp(poly a){
        if(a.size()==1)return a[0]=1,a;
        int n=a.size(),m=a.size()+1>>1;
        poly _a(m);
        for(int i=0;i<m;++i)_a[i]=a[i];
        poly b=getexp(_a);
        return mul(b,poly(1,1)-getln(b,a.size())+a,n);
    }
    il poly Pow(poly a,int b,int n){
        int cnt0=0;while(!a[cnt0]&&cnt0<a.size())++cnt0;
        a.erase(a.begin(),a.begin()+cnt0);
        if(1ll*b*cnt0>=n)return poly(n,0);
        int m=a.size(),a0=a[0],ia0=pow(a[0],mod-2);
        for(int i=0;i<m;++i)a[i]=1ll*ia0*a[i]%mod;
        a=getexp(getln(a)*b);a.resize(n);
        poly ret(n,0);int p=1ll*cnt0*b;
        for(int i=0;i<m;++i){
            ret[p++]=1ll*a[i]*a0%mod;
            if(p==n)break;
        }
        return ret;
    }
    il poly operator%(poly a,poly b){
        int n=a.size(),m=b.size();
        if(n<m)return a;
        std::reverse(a.begin(),a.end());
        std::reverse(b.begin(),b.end());
        b.resize(n);
        poly c=mul(a,getinv(b),n-m+1);
        std::reverse(a.begin(),a.end());
        b.resize(m);std::reverse(b.begin(),b.end());
        std::reverse(c.begin(),c.end());
        a=a-mul(b,c);
        int s=a.size();while(s&&!a[s-1])--s;
        a.resize(s);
        return a;
    }
    namespace sqrt_mod{
        struct numb{int real,imag;};
        int a,w;
        il numb operator*(const numb&a,const numb&b){return(numb){(1ll*a.real*b.real+1ll*a.imag*b.imag%mod*w)%mod,(1ll*a.real*b.imag+1ll*a.imag*b.real)%mod};}
        il int sqrt(int x){
            while(1){
                a=rand()%(mod-1)+1;w=(1ll*a*a-x+mod)%mod;
                if(pow(w,mod-1>>1)!=1)break;
            }
            numb s={a,1},ret={1,0};int y=mod+1>>1;
            while(y){
                if(y&1)ret=ret*s;
                s=s*s;y>>=1;
            }
            return std::min(mod-ret.real,ret.real);
        }
    }
    il poly sqrt(poly a){
        if(a.size()==1)return a[0]=sqrt_mod::sqrt(a[0]),a;
        int n=a.size(),m=a.size()+1>>1;
        poly _a(m);
        for(int i=0;i<m;++i)_a[i]=a[i];
        poly b=sqrt(_a);b.resize(n);
        return (b+mul(a,getinv(b),n))*(mod+1>>1);
    }
    il vd poly_init(){
        int G=Gmod,iG=pow(G,mod-2);
        for(int i=1;i<maxn;i<<=1)P[i]=pow(G,(mod-1)/(i<<1)),iP[i]=pow(iG,(mod-1)/(i<<1));
    }
    struct _poly_auto_init{_poly_auto_init(){poly_init();}}_auto_init;
    //End==============多项式板子部分===================
    int fact[131113],ifact[131113];
    int main(){
    #ifdef XZZSB
        freopen("in.in","r",stdin);
        freopen("out.out","w",stdout);
    #endif
        int n=gi(),k=gi();
        if(k>n){++n;while(n--)printf("0 ");return 0;}
        fact[0]=1;for(int i=1;i<=n;++i)fact[i]=1ll*fact[i-1]*i%mod;
        ifact[n]=pow(fact[n],mod-2);for(int i=n;i;--i)ifact[i-1]=1ll*ifact[i]*i%mod;
        poly f(2,1);
        f=Pow(getln(f,n-k+1),k,n+1);
        for(int i=0;i<k;++i)printf("0 ");
        for(int i=k;i<=n;++i)printf("%d ",1ll*f[i]*ifact[k]%mod*fact[i]%mod*(i-k&1?mod-1:1)%mod);
        return 0;
    }
    

    第二类斯特林数

    (egin{Bmatrix} n\mend{Bmatrix})表示将(n)个不同元素放进(m)个相同盒子的方案数。

    这个东西好像比第一类水一点,题也多一点?

    递推式:(egin{Bmatrix} n\mend{Bmatrix}=egin{Bmatrix} n-1\mend{Bmatrix}m+egin{Bmatrix} n-1\m-1end{Bmatrix})

    通项:容斥枚举放了(i)个空盒子。(sum_{i=0}^mfrac{(-1)^i}{i!}frac{(m-i)^n}{(m-i)!})

    性质:分解(m^n),枚举有一个盒子非空。(m^n=sum_{i=0}^mC_m^iegin{Bmatrix} n\iend{Bmatrix}i!)

    求一行

    直接用通项ntt

    求一列

    https://www.luogu.org/problemnew/show/P5396

    构造生成函数(F_m(x)=sum_{i=m}^{infty}egin{Bmatrix} i\mend{Bmatrix}x^i)

    根据递推式可得(F_m(x)=F_{m-1}(x)cdot x+F_m(x)cdot mx)

    移项得(F_m(x)=frac{F_{m-1}(x)}{1-mx}x)

    显然(F(0)=1),所以可以一路转移过去,最后写出来是这样的

    (F_m(x)=frac{x^m}{(1-x)(1-2x)cdots(1-mx)})

    可以分治ntt解决。

    具体数学上说(egin{bmatrix} n\mend{bmatrix}=egin{Bmatrix} -m\-nend{Bmatrix}),那么应该第2类的列和第1类的行有点关系,然而好像不可以用那种倍增方法优化,因为((1-x)cdots(1-nx))并不能修改x推到([1-(n+1)x]cdots(1-2nx))

    好像只能分治ntt了233

    一些题

    鸡你太美题

  • 相关阅读:
    div标签的闭合检查
    jquery easyui 显示和关闭数据加载的遮罩
    codeforces 446A DZY Loves Sequences
    android高速开发框架xUtils
    Android-spinner
    遗传算法优化策略
    面向对象的勾勾画画
    Android studio 解决setText中文乱码问题
    CAS—改动默认登录页
    android 使用post 提交
  • 原文地址:https://www.cnblogs.com/xzz_233/p/11008681.html
Copyright © 2011-2022 走看看