标签：csdn csharp stream include namespace src int type img
题意：
给你两个公式：
对于任意质数p:
让你求Bell_n模95041567，n<=2147483647；
其中95041567=31*37*41*43*47;
考虑这些质数都很小，那么我们可以考虑对每一个质数分别处理，然后用CRT合并；
然后对于每个质数，我们发现第二个公式是一个类似斐波那契的递推，我们考虑用一个50*50的矩阵乘法优化递推；
下面给出CRT合并的公式(蒯自ACdreamer)：
设正整数
有整数解。并且在模
其中
//MADE BY QT666
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
typedef long long ll;
const int N=100050;
const int Mod=95041567;
int m[]={31,37,41,43,47};
ll c[55][55],bell[55],ans[3][55][55],s[55][55],a[N];
void work(int a,int b,int mo){
for(int i=1;i<=50;i++){
for(int j=1;j<=50;j++){
for(int k=1;k<=50;k++){
(s[i][j]+=(ans[a][i][k]*ans[b][k][j])%mo)%=mo;
}
}
}
for(int i=1;i<=50;i++){
for(int j=1;j<=50;j++){
ans[a][i][j]=s[i][j],s[i][j]=0;
}
}
}
ll calc(int k,int y){
memset(ans,0,sizeof(ans));
memset(bell,0,sizeof(bell));
memset(c,0,sizeof(c));
for(int i=0;i<=50;++i) c[i][0]=1;
for(int i=1;i<=50;++i)
for(int j=1;j<=i;++j){
c[i][j]=(c[i-1][j-1]+c[i-1][j])%m[k];
}
bell[0]=1;
for(int i=1;i<=50;i++){
for(int j=0;j<i;j++) (bell[i]+=(bell[j]*c[i-1][j])%m[k])%=m[k];
}
for(int i=1;i<=50;i++) ans[1][1][i]=bell[i];
ans[2][50-m[k]+1][50]=1,ans[2][50-m[k]+2][50]=1;
for(int i=1;i<=49;i++) ans[2][i+1][i]=1;
if(y<=50) return ans[1][1][y];
else{
y-=50;
while(y){
if(y&1) work(1,2,m[k]);
work(2,2,m[k]);y>>=1;
}
return ans[1][1][50];
}
}
void exgcd(ll a,ll b,ll &x,ll &y){
if(!b) {x=1,y=0;return;}
exgcd(b,a%b,y,x),y-=x*(a/b);
}
ll qpow(ll x,ll y,ll mo){
ll ret=1;
while(y){
if(y&1) (ret*=x)%=mo;
(x*=x)%=mo;y>>=1;
}
return ret;
}
void Crt(int n){
ll Ans=0;
for(int i=0;i<5;i++){
ll M=Mod/m[i];
(Ans+=qpow(M,m[i]-2,m[i])*M%Mod*calc(i,n)%Mod)%=Mod;
}
printf("%lld\n",Ans);
}
int main(){
freopen("bell.in","r",stdin);
freopen("bell.out","w",stdout);
int T;scanf("%d",&T);
while(T--){
int n;scanf("%d",&n);Crt(n);
}
return 0;
}
标签：csdn csharp stream include namespace src int type img
原文地址：http://www.cnblogs.com/qt666/p/7632004.html
