2013年长春网络赛的题。
比赛时很快就找到了bell数的性质,看到数据范围也马上想到要用矩阵做,只可惜当时没接触过中国剩余定理……
以下题解转自:http://fudq.blog.163.com/blog/static/19135023820139114642927
对于贝尔数,如果p是任意质数,则有B[p+n]=B[n]+B(n+1)%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]。
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
const int N=55;
int stir2[5][N][N]; //第二类斯特林数
int bell[5][N]; //贝尔数
int data[5];
int mod[] = {31, 37, 41, 43, 47};
int modd; //矩阵计算时的模值
class Matrix
{
public:
int a[N][N];
int n; //矩阵大小
void init (int x)
{
memset(a,0,sizeof(a));
if (x) for (int i=0;i<N;i++)
a[i][i]=1;
}
Matrix operator * (Matrix b)
{
Matrix p;
p.n=b.n; //
p.init(0);
for (int i=0;i<n;i++) for (int j=0;j<n;j++)
for (int k=0;k<n;k++)
(p.a[i][j]+=a[i][k]*b.a[k][j]%modd)%=modd;
return p;
}
Matrix power (int t)
{
Matrix ans,p=*this;
ans.n=p.n; //
ans.init(1);
while (t)
{
if (t&1) ans=ans*p;
p=p*p;
t>>=1;
}
return ans;
}
}a;
void Init () //预处理50以下的贝尔数
{
int i,j,k;
for (i=0;i<5;i++)
{
stir2[i][0][0]=1;
bell[i][0]=bell[i][1]=1;
}
for (i=1;i<=50;i++)
{
for (j=0;j<5;j++)
stir2[j][i][1]=stir2[j][i][i]=1;
for (j=1;j<=i;j++) for (k=0;k<5;k++)
stir2[k][i][j]=(stir2[k][i-1][j-1]+j*stir2[k][i-1][j])%mod[k];
}
for (i=1;i<=50;i++)
for (k=0;k<5;k++)
{
bell[k][i]=0;
for (j=0;j<=i;j++)
bell[k][i]=(bell[k][i]+stir2[k][i][j])%mod[k];
}
}
int Extended_Euclid (int a,int b,int &x,int &y)
{//扩展欧几里得算法,求ax+by=gcd(a,b)的一组解(x,y),d=gcd(a,b)
int d;
if (b==0)
{
x=1;y=0;
return a;
}
d=Extended_Euclid(b,a%b,y,x);
y-=a/b*x;
return d;
}
int Chinese_Remainder (int a[],int w[],int len)
{//中国剩余定理 a[]存放余数 w[]存放两两互质的数
int x,y,m;
int lcm=1,i,ans=0;
for (i=0;i<len;i++)
lcm*=w[i];
for (i=0;i<len;i++)
{
m=lcm/w[i];
Extended_Euclid(w[i],m,x,y);
ans=(ans+y*m*a[i])%lcm;
}
return (lcm+ans%lcm)%lcm;
}
int Cal (int n,int t) //构造矩阵a,计算bell[n]%w[t]
{
int ans=0,i;
a.init(0);
a.n=mod[t];
modd=mod[t];
for (i=1;i<mod[t];i++)
a.a[i][i]=a.a[i-1][i]=1;
a.a[ mod[t]-1 ][0]=1;
int p1=n/(mod[t]-1);
int p2=n%(mod[t]-1);
a=a.power(p1);
for (i=0;i<mod[t];i++)
ans=(ans+bell[t][i]*a.a[i][p2])%mod[t];
return ans;
}
int main ()
{
#ifdef ONLINE_JUDGE
#else
freopen("read.txt","r",stdin);
#endif
int T,n,i;
scanf("%d",&T);
Init ();
while (T--)
{
scanf("%d",&n);
if (n<50)
{
for (i=0;i<5;i++)
data[i]=bell[i][n];
printf("%d\n",Chinese_Remainder(data,mod,5));
continue;
}
for (i=0;i<5;i++)
data[i]=Cal(n,i);
printf("%d\n",Chinese_Remainder(data,mod,5));
}
return 0;
}