HDU4767 Bell Bell数的性质，矩阵快速幂

 

Description

 

What? MMM is learning Combinatorics!?
Looks like she is playing with the bell sequence now:
bell[n] = number of ways to partition of the set {1, 2, 3, ..., n}

e.g. n = 3:
{1, 2, 3}
{1} {2 3}
{1, 2} {3}
{1, 3} {2}
{1} {2} {3}
so, bell[3] = 5

MMM wonders how to compute the bell sequence efficiently.

Input

T -- number of cases
for each case:
　　n (1 <= n < 2^31)

Output

for each case:
　　bell[n] % 95041567

Sample Input

6
1
2
3
4
5
6

Sample Output

1
2
5
15
52
203

题意：输出$Bell$数对$95041567$取模

本身并不是对质数取模，一开始可以对模数分解

这样子所有的质数模就不超过50

考虑公式$ B_{n+p}equiv B_{n}+B_{n+1}mod p$

我们可以尝试构造一个矩阵

使用矩阵快速幂加速

最后用中国剩余定理合并即可

#include<bits/stdc++.h>
using namespace std;

template<typename T>T read(){
    T x=0,f=1;char c=getchar();
    
    for(;!isdigit(c);c=getchar())if(c=='-')f=-1;
    for(;isdigit(c);c=getchar())x=(x<<3)+(x<<1)+(c^48);

    return x*f;
}

const int maxn=100010;
typedef long long ll;

const int mod[5]={47,43,41,37,31};
const int Mod=95041567;

int n,Bell[60],C[60],kase;

int A[50];

struct matrix{
    int m[50][50];
    matrix(){
        memset(m,0,sizeof m);
    };
    matrix(int v){
        memset(m,0,sizeof m);
        for(int i=0;i<50;++i)m[i][i]=v;
    }
};

matrix mul(matrix a,matrix b){
    matrix ret;
    for(int i=0;i<50;++i){
        for(int j=0;j<50;++j){
            for(int k=0;k<50;++k){
                (ret.m[i][j]+=1ll*a.m[i][k]*b.m[k][j]%mod[kase]);
                (ret.m[i][j])%=mod[kase];
            }
        }
    }
    return ret;
}

matrix fpm(matrix a,ll b){
    matrix ret(1);
    for(;b;b>>=1,a=mul(a,a))
        if(b&1)ret=mul(ret,a);
    return ret;
}

int fpm(int a,int i){
    int ret=1,b=mod[i]-2;
    for(;b;b>>=1,a=a*a%mod[i])
        if(b&1)ret=ret*a%mod[i];
    return ret;
}

int CRT(){
    int ret=0;
    for(int i=0;i<5;++i){
        int mm=Mod/mod[i];
        int inv=fpm(mm%mod[i],i);
        ret=(ret+1ll*A[i]*mm*inv)%Mod;
    }
    return ret;
}

void solve(){
    scanf("%d",&n);
    if(n<48){
        printf("%d
",Bell[n]);
        return;
    }
    for(kase=0;kase<5;++kase){
        matrix tmp,vec;
        
        for(int i=0;i<mod[kase];++i)
            vec.m[i][0]=Bell[i]%mod[kase];
        for(int i=0;i<mod[kase]-1;++i)
            tmp.m[i][i]=tmp.m[i][i+1]=1;
        tmp.m[mod[kase]-1][mod[kase]-1]=1;
        tmp.m[mod[kase]-1][0]=1;
        tmp.m[mod[kase]-1][1]=1;
        
        tmp=fpm(tmp,n/mod[kase]);
        vec=mul(tmp,vec);

        A[kase]=vec.m[n%mod[kase]][0];
    }

    printf("%d
",CRT());
}

void init(){
    Bell[0]=Bell[1]=1;
    C[0]=1;
    for(int i=2;i<50;++i){
        C[i-1]=Bell[i-1];
        for(int j=i-2;~j;--j)
            C[j]=(C[j]+C[j+1])%Mod;
        Bell[i]=C[0];
    }

}

int main(){
    #ifndef ONLINE_JUDGE
    freopen("chk.in","r",stdin);
    freopen("chk.out","w",stdout);
    #endif
    
    init();
    for(int T=read<int>()+1;--T;){
        solve();
    }
    return 0;
}