zoukankan      html  css  js  c++  java
  • luoguP5824 十二重计数法 组合+生成函数+第二类斯特林数

    12种组合计数问题合在一起.   

    code:

    #include <cmath>
    #include <cstring>
    #include <algorithm>
    #include <cstdio>
    #include <string>
    #define ll long long
    #define ull unsigned long long
    using namespace std;
    namespace IO
    {
        char buf[100000],*p1,*p2;
        #define nc() (p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++)
        int rd()
        {
            int x=0; char s=nc();
            while(s<'0') s=nc();
            while(s>='0') x=(((x<<2)+x)<<1)+s-'0',s=nc();
            return x;
        }    
        void print(int x) {if(x>=10) print(x/10);putchar(x%10+'0');}
        void setIO(string s)
        {
            string in=s+".in";
            string out=s+".out";
            freopen(in.c_str(),"r",stdin);
            // freopen(out.c_str(),"w",stdout);
        }
    }; 
    const int G=3;
    const int N=2000005;      
    const int mod=998244353;                   
    int A[N],B[N],w[2][N],mem[N*100],*ptr=mem,tmpa[N],tmpb[N],aa[N],bb[N];     
    inline int qpow(int x,int y)
    {
        int tmp=1;
        for(;y;y>>=1,x=(ll)x*x%mod)     if(y&1) tmp=(ll)tmp*x%mod;
        return tmp;
    } 
    inline int INV(int a) { return qpow(a,mod-2); }   
    inline void ntt_init(int len)
    {
        int i,j,k,mid,x,y; 
        w[1][0]=w[0][0]=1,x=qpow(3,(mod-1)/len),y=qpow(x,mod-2);
        for (i=1;i<len;++i) w[0][i]=(ll)w[0][i-1]*x%mod,w[1][i]=(ll)w[1][i-1]*y%mod;    
    }
    void NTT(int *a,int len,int flag)
    {
        int i,j,k,mid,x,y;          
        for(i=k=0;i<len;++i)
        {
            if(i>k)    swap(a[i],a[k]);
            for(j=len>>1;(k^=j)<j;j>>=1);
        }
        for(mid=1;mid<len;mid<<=1)       
            for(i=0;i<len;i+=mid<<1)
                for(j=0;j<mid;++j)     
                {
                    x=a[i+j], y=(ll)w[flag==-1][len/(mid<<1)*j]*a[i+j+mid]%mod;
                    a[i+j]=(x+y)%mod;
                    a[i+j+mid]=(x-y+mod)%mod;
                }
        if(flag==-1)
        {
            int rev=INV(len);
            for(i=0;i<len;++i)    a[i]=(ll)a[i]*rev%mod;
        }
    }         
    inline void getinv(int *a,int *b,int len,int la)      
    {
        if(len==1) { b[0]=INV(a[0]);   return; }
        getinv(a,b,len>>1,la);
        int l=len<<1,i;
        memset(A,0,l*sizeof(A[0]));         
        memset(B,0,l*sizeof(A[0]));
        memcpy(A,a,min(la,len)*sizeof(a[0]));                                                 
        memcpy(B,b,len*sizeof(b[0]));        
        ntt_init(l);
        NTT(A,l,1),NTT(B,l,1); 
        for(i=0;i<l;++i)  A[i]=((ll)2-(ll)A[i]*B[i]%mod+mod)*B[i]%mod;
        NTT(A,l,-1);                            
        memcpy(b,A,len<<2);     
    }    
    void get_dao(int *a,int *b,int len)
    {
        for(int i=1;i<len;++i) b[i-1]=(ll)i*a[i]%mod;
        b[len-1]=0;
    }                
    void get_jifen(int *a,int *b,int len)
    {
        for(int i=1;i<len;++i) b[i]=(ll)INV(i)*a[i-1]%mod;
        b[0]=0;
    }
    void get_ln(int *a,int *b,int len,int la)
    {
        int l=len<<1,i;
        memset(tmpa,0,l<<2);
        memset(tmpb,0,l<<2);
        get_dao(a,tmpa,min(len,la));
        getinv(a,tmpb,len,la);
        ntt_init(l);      
        NTT(tmpa,l,1),NTT(tmpb,l,1);
        for(i=0;i<l;++i) tmpa[i]=(ll)tmpa[i]*tmpb[i]%mod;
        NTT(tmpa,l,-1);
        get_jifen(tmpa,b,len);  
    }   
    void get_exp(int *a,int *b,int len,int la)
    {
        if(len==1) { b[0]=1; return; }                       
        int l=len<<1,i;
        get_exp(a,b,len>>1,la);           
        for(i=0;i<l;++i)  aa[i]=bb[i]=0;
        for(i=0;i<(len>>1);++i) aa[i]=b[i];        
        get_ln(b,bb,len,len>>1);                                           
        for(i=0;i<len;++i) bb[i]=(ll)(mod-bb[i]+(i>=la?0:a[i]))%mod;                           
        bb[0]=(bb[0]+1)%mod;
        ntt_init(l);
        NTT(aa,l,1),NTT(bb,l,1);
        for(i=0;i<l;++i) aa[i]=(ll)aa[i]*bb[i]%mod; 
        NTT(aa,l,-1);  
        for(i=0;i<len;++i)  b[i]=aa[i]; 
    }
    struct poly
    {
        int len,*a;
        poly(){}  
        poly(int l) {len=l,a=ptr,ptr+=l; }       
        inline void rev() { reverse(a,a+len); }  
        inline void fix(int l) {len=l,a=ptr,ptr+=l;}
        inline void get_mod(int l) { for(int i=l;i<len;++i) a[i]=0;  len=l;  }
        inline poly dao()
        {   
            poly re(len-1);
            for(int i=1;i<len;++i)  re.a[i-1]=(ll)i*a[i]%mod;    
            return re;
        }
        inline poly jifen()
        {
            poly c;
            c.fix(len+1);   
            c.a[0]=0;
            for(int i=1;i<=len;++i) c.a[i]=(ll)a[i-1]*INV(i)%mod;           
            return c;
        }      
        inline poly Inv(int l)
        {               
            int lim=1;  
            while(lim<=l) lim<<=1;  
            poly b(lim);
            getinv(a,b.a,lim,len);                             
            b.get_mod(l); 
            return b;                   
        }              
        inline poly ln(int l)
        {
            int lim=1;
            while(lim<=l) lim<<=1;                 
            poly b(lim);
            get_ln(a,b.a,lim,len);
            return b;
        }                     
        inline poly exp(int l)
        {
            int lim=1;
            while(lim<=l) lim<<=1;
            poly b(lim);
            get_exp(a,b.a,lim,len);    
            return b;
        }  
        inline poly q_pow(int k,int l)
        {          
            int lim=1; 
            while(lim<=l) lim<<=1; 
            poly b(lim),c(lim); 
            get_ln(a,b.a,lim,len);   
            for(int i=0;i<b.len;++i) b.a[i]=(ll)b.a[i]*k%mod;    
            get_exp(b.a,c.a,lim,b.len);  
            c.get_mod(l);
            return c; 
        }                   
        inline poly operator*(const poly &b) const
        {
            poly c(len+b.len-1);
            if(c.len<=500)
            {    
                for(int i=0;i<len;++i)
                    if(a[i])   for(int j=0;j<b.len;++j)  c.a[i+j]=(c.a[i+j]+(ll)(a[i])*b.a[j])%mod; 
                return c;
            }
            int n=1;
            while(n<(len+b.len)) n<<=1;
            memset(A,0,n<<2);
            memset(B,0,n<<2);
            memcpy(A,a,len<<2);                        
            memcpy(B,b.a,b.len<<2);             
            ntt_init(n);   
            NTT(A,n,1), NTT(B,n,1);
            for(int i=0;i<n;++i) A[i]=(ll)A[i]*B[i]%mod;
            NTT(A,n,-1);
            memcpy(c.a,A,c.len<<2);
            return c;  
        }
        poly operator+(const poly &b) const
        {
            poly c(max(len,b.len));
            for(int i=0;i<c.len;++i)  c.a[i]=((i<len?a[i]:0)+(i<b.len?b.a[i]:0))%mod;
            return c;
        }
        poly operator-(const poly &b) const
        {
            poly c(len);  
            for(int i=0;i<len;++i)
            {
                if(i>=b.len)   c.a[i]=a[i];
                else c.a[i]=(a[i]-b.a[i]+mod)%mod;
            }
            return c;
        }
        poly operator/(poly u)
        {
            int n=len,m=u.len,l=1;
            while(l<(n-m+1)) l<<=1;                      
            rev(),u.rev();       
            poly v=u.Inv(l);
            v.get_mod(n-m+1);   
            poly re=(*this)*v;
            rev(),u.rev();
            re.get_mod(n-m+1);    
            re.rev();
            return re;
        } 
        poly operator%(poly u)
        { 
            poly re=(*this)-u*(*this/u);   
            re.get_mod(u.len-1);  
            return re;
        }                
    }Stirling,Ferrers;     
    #define MAX 500002   
    int n,m;            
    int fac[N],inv[N],in[N];     
    int C(int x,int y) 
    {
        if(x<0||y<0||x<y) return 0;  
        return (ll)fac[x]*inv[y]%mod*inv[x-y]%mod;     
    }
    void init()  
    { 
        int i,j;  
        fac[0]=inv[0]=in[0]=1;   
        for(i=1;i<MAX;++i) fac[i]=(ll)fac[i-1]*i%mod,inv[i]=INV(fac[i]),in[i]=INV(i);          
    }
    void case_1() { printf("%d
    ",qpow(m,n)); }  
    void case_2() { printf("%d
    ",(ll)C(m,n)*fac[n]%mod); }  
    void get_Stirling() 
    {
        int i,j;   
        poly a(n+1),b(n+1);   
        for(i=0;i<a.len;++i) 
        {
            a.a[i]=(ll)((i&1)?mod-1:1)*inv[i]%mod;   
            b.a[i]=(ll)qpow(i,n)*inv[i]%mod;   
        }
        Stirling=a*b;           
    }
    void case_3() { printf("%d
    ",m<=n?(ll)Stirling.a[m]*fac[m]%mod:0); }
    void case_4() 
    { 
        int ans=0,i,j;    
        for(i=1;i<=min(n,m);++i) (ans+=Stirling.a[i])%=mod;   
        printf("%d
    ",ans);   
    }
    void case_5() { printf("%d
    ",m<n?0:1);  } 
    void case_6() { printf("%d
    ",m<=n?Stirling.a[m]:0); }      
    void case_7() { printf("%d
    ",C(m+n-1,n)); }  
    void case_8() { printf("%d
    ",C(m,n)); }
    void case_9() { printf("%d
    ",m>n?0:C(n-1,n-m)); }   
    void get_Ferrers() 
    {   
        int i,j;  
        Ferrers.fix(n+1);     
        for(i=1;i<=m;++i)  
        {
            for(j=i;j<=n;j+=i) 
            {
                int t=j/i;
                (Ferrers.a[j]+=in[t])%=mod;    
            }
        }    
        Ferrers=Ferrers.exp(n+1);             
    }
    void case_10() { printf("%d
    ",Ferrers.a[n]);} 
    void case_11() { printf("%d
    ",n<=m); }    
    void case_12() { printf("%d
    ",m>n?0:Ferrers.a[n-m]); }         
    int main()  
    {
        // IO::setIO("input");   
        int i,j;  
        scanf("%d%d",&n,&m);  
        init();  
        get_Stirling(); 
        case_1(); 
        case_2();       
        case_3(); 
        case_4();     
        case_5(); 
        case_6();               
        case_7();    
        case_8();   
        case_9();    
        get_Ferrers();    
        case_10();    
        case_11();   
        case_12();     
        return 0;
    }   
    

      

  • 相关阅读:
    Java实现web页面内容抓取
    Java写入文件的几种方法及性能对比
    Java实现导出excel
    win10系统安装VMware虚拟机软件以及linux系统
    oracle11g安装教程
    oracle怎么建立本地连接
    工厂模式
    计算一个字符串中每个字符出现的次数
    MySql多表查询
    如何查看MySql的sql语句性能
  • 原文地址:https://www.cnblogs.com/guangheli/p/12259496.html
Copyright © 2011-2022 走看看