注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

fudq's AC Road

何以解忧,唯有AC!

 
 
 

日志

 
 

hdu 4767 Bell  

2013-10-01 13:58:49|  分类: ACM-hdu |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
http://acm.hdu.edu.cn/showproblem.php?pid=4767
题意:求第n个Bell数,n最大2^31。结果对95041567取余。
题解:
看维基百科(http://zh.wikipedia.org/wiki/%E8%B4%9D%E5%B0%94%E6%95%B0)了解关于Bell数的几条性质:
每个Bell数都是第二类斯特林数之和
如果p是任意质数,则有 B_{p+n}\equiv B_n+B_{n+1}\ (\operatorname{mod}\ p).
还有可以发现的是95041567=31*37*41*43*47。
利用以上信息我们便可解此题。
如果能够得到B(n)%31,B(n)%37,B(n)%41,B(n)%43,B(n)%47的结果,我们便可以用中国剩余定理求出B(n)%95041567的结果
现在问题转换成求B(n)%w[i],由于w[i]是质数,便可利用上面Bell数的性质。
举个例子,如果w[i]=5,可以由B0 B1 B2 B3 B4得到B5 B6 B7 B8,又由B4 B5 B6 B7 B8得到B9 B10 B11 B12,依次类推。
这时可以利用矩阵二次幂快速求解,构造一个如下所示的01矩阵:
0 1 0 0 0
0 1 1 0 0
0 0 1 1 0
0 0 0 1 1
1 0 0 0 1
便可以由B0 B1 B2 B3 B4得到B4 B5 B6 B7 B8。
对于一个n,需要将构造矩阵乘上n/(w[i]-1)次,第n%(w[i]-1)列的结果便是B(n)%w[i]。

/*
* pro.cpp
*
* Created on: 2013-09-30
* Author: fudq
*/
#include <functional>
#include <algorithm>
#include <iostream>
//#include <fstream>
#include <sstream>
#include <iomanip>
#include <numeric>
#include <cstring>
#include <cassert>
#include <cstdio>
#include <string>
#include <vector>
#include <bitset>
#include <queue>
#include <stack>
#include <cmath>
#include <ctime>
#include <list>
#include <set>
#include <map>
using namespace std;
//#pragma comment(linker,"/STACK:102400000,102400000")

#define FOR(i,a) for((i)=0;i<(a);(i)++)
#define MEM(a) (memset((a),0,sizeof(a)))
#define LL long long

const int N=55;
const int M=32010;
const int MOD=95041567;
const int INF=0x7fffffff;
const int dir[4][2]={{-1,0},{1,0},{0,-1},{0,1}};
//const int dir[8][2]={{-1,0},{1,0},{0,-1},{0,1},{-1,1},{1,-1},{-1,-1},{1,1}};
const double eps=1e-16;
const double PI=acos(-1.0);

inline int sign(double x){return (x>eps)-(x<-eps);}
template<class T> T gcd(T a,T b){return b?gcd(b,a%b):a;}
template<class T> T lcm(T a,T b){return a/gcd(a,b)*b;}
template<class T> inline T lcm(T a,T b,T d){return a/d*b;}
template<class T> inline T Min(T a,T b){return a<b?a:b;}
template<class T> inline T Max(T a,T b){return a>b?a:b;}
template<class T> inline T Min(T a,T b,T c){return min(min(a, b),c);}
template<class T> inline T Max(T a,T b,T c){return max(max(a, b),c);}
template<class T> inline T Min(T a,T b,T c,T d){return min(min(a, b),min(c,d));}
template<class T> inline T Max(T a,T b,T c,T d){return max(max(a, b),max(c,d));}
/*************************/

LL bell[N][6],num[N][N][6];
LL a[6],w[6]={31,37,41,43,47};

void init() //打印前50第二类斯特林数,每行相加即为贝尔数
{
for(int i=0;i<5;i++)
{
num[0][0][i]=1;
bell[0][i]=bell[1][i]=1;
}
for(int i=1;i<=50;i++)
{
for(int j=0;j<5;j++)
num[i][1][j]=1;
for(int j=1;j<=i;j++)
for(int k=0;k<5;k++)
num[i][j][k]=(num[i-1][j-1][k]+j*num[i-1][j][k])%w[k];
for(int j=0;j<5;j++)
{
bell[i][j]=0;
for(int k=1;k<=i;k++)
bell[i][j]=(bell[i][j]+num[i][k][j])%w[j];
}
}
}

LL Extended_Euclid(LL a,LL b,LL &x,LL &y) //扩展欧几里得算法
{
LL d;
if(b==0)
{
x=1;y=0;
return a;
}
d=Extended_Euclid(b,a%b,y,x);
y-=a/b*x;
return d;
}

LL Chinese_Remainder(LL a[],LL w[],LL len) //中国剩余定理 a[]存放余数 w[]存放两两互质的数
{
LL d,x,y,m,s=1,ret=0;
for(int i=0;i<len;i++)
s*=w[i];
for(int i=0;i<len;i++)
{
m=s/w[i];
d=Extended_Euclid(w[i],m,x,y);
ret=(ret+y*m*a[i])%s;
}
return (s+ret%s)%s;
}

struct Matrix{
int v[N][N];
};

Matrix Mat_mul(Matrix m1,Matrix m2,LL pri) //矩阵相乘
{
Matrix c;
MEM(c.v);
for(int i=0;i<N;i++)
for(int j=0;j<N;j++)
for(int k=0;k<N;k++)
c.v[i][j]=(c.v[i][j]+m1.v[i][k]*m2.v[k][j])%pri;
return c;
}

Matrix Mpow(Matrix A,LL n,LL pri) //矩阵快速幂
{
Matrix c,x=A;
MEM(c.v);
for(int i=0;i<N;i++)
c.v[i][i]=1;
while(n >= 1)
{
if(n & 1)
c=Mat_mul(c,x,pri);
n>>=1;
x=Mat_mul(x,x,pri);
}
return c;
}

LL work(LL n,LL t) //计算B(n)%w[t]
{
LL res=0;
Matrix tmp,ans;
MEM(tmp.v); //构造定义矩阵tmp
for(int i=1;i<w[t];i++)
tmp.v[i][i]=tmp.v[i-1][i]=1;
tmp.v[w[t]-1][0]=1;
LL p1,p2;
p1=n/(w[t]-1);p2=n%(w[t]-1);
ans=Mpow(tmp,p1,w[t]);
for(int i=0;i<w[t];i++)
res=(res+bell[i][t]*ans.v[i][p2])%w[t];
return res;
}

void solve(LL n)
{
if(n < 50)
{
for(int i=0;i<5;i++)
a[i]=bell[n][i];
printf("%I64d\n",Chinese_Remainder(a,w,5));
return ;
}
for(int i=0;i<5;i++) //计算对每个数的余数
a[i]=work(n,i);
printf("%I64d\n",Chinese_Remainder(a,w,5));
}

int main()
{
#ifndef ONLINE_JUDGE
freopen("testin.txt", "r", stdin);
//freopen("textout.txt", "w", stdout);
#endif
init();
int T;
LL n;
scanf("%d",&T);
while(T--)
{
scanf("%I64d",&n);
solve(n);
}
return 0;
}



  评论这张
 
阅读(206)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017