链接:http://acm.hdu.edu.cn/showproblem.php?pid=1588
题意很好理解,不多说。现在开始做矩阵专题,都是入门类型的,但是对于矩阵的构造还是很难想到,所以这个题目也想了两天。主要想说说我在做这题时构造矩阵的想法。
首先,斐波拉契数列是可以写成矩阵形式的,这应该都懂,所以无形之中给了矩阵构造的提示。下边就是我做完后整理的思路。原谅菜鸟的思路
代码:
#include <string.h>
#include <stdio.h>
//#define mod 2008
#define MAX_dimension 2
typedef int MATRIX_TYPE;
typedef __int64 MAX_INT_TYPE; //矩阵乘法运算时和加法运算时,为避免溢出,在取模之前用这个类型暂存。
typedef MATRIX_TYPE Matrix[MAX_dimension][MAX_dimension]; //矩阵类型的定义
int ndim=MAX_dimension;
int mod;
Matrix A={ {1,1},
{1,0}};
void m_zero(Matrix x)
{
memset(x, 0, sizeof(MATRIX_TYPE)*MAX_dimension*MAX_dimension);
}
void m_one(Matrix x)
{
int i;
m_zero(x);
for(i=0;i<ndim;++i)
x[i][i]=1;
}
void m_copy(Matrix src,Matrix dest)
{
memcpy(dest,src, sizeof(MATRIX_TYPE)*MAX_dimension*MAX_dimension);
}
//z=x+y;
void m_add(Matrix x,Matrix y,Matrix z)
{
int i,j;
for(i=0;i<ndim;i++)
for(j=0;j<ndim;j++)
if(mod<=1)z[i][j]=x[i][j]+y[i][j];
else z[i][j]=(x[i][j]+(MAX_INT_TYPE)y[i][j])%mod;
}
//c=a*b
void m_multiple(Matrix a,Matrix b,Matrix c)
{
int i,j,k;
MAX_INT_TYPE t;
for(i=0;i<ndim;i++)
for(j=0;j<ndim;j++)
{
c[i][j]=0;
if(mod<=1)
for(k=0;k<ndim;k++) c[i][j]+=a[i][k]*b[k][j];
else
for(k=0;k<ndim;k++)
{
c[i][j]+=(a[i][k]*(MAX_INT_TYPE)b[k][j])%mod;
c[i][j]%=mod;
}
}
}
//y=x^n
void m_pow(Matrix x, unsigned int n, Matrix y)
{
Matrix temp,temp_x;
m_one(y);
m_copy(x,temp_x);
while(n)
{
if (n & 1)
{
m_multiple(temp_x,y,temp);
m_copy(temp,y);
}
if ((n >>= 1))
{
m_multiple(temp_x,temp_x,temp);
m_copy(temp,temp_x);
}
}
}
//根据已经算出来的x^1 x^2 x^4 ...去求矩阵的幂y=x^n
void m_pow_with_saved(Matrix x_p[],unsigned int n, Matrix y)
{
Matrix temp;
m_one(y);
for(int k=1;n;++k,n>>=1)
{
if (n & 1)
{
m_multiple(x_p[k],y,temp);
m_copy(temp,y);
}
}
}
//计算根据已经算出来的x^1 x^2 x^4 ...去求矩阵的幂和y=sum(x^n)
void m_pow_sum1(Matrix x_p[],unsigned int n, Matrix y)
{
m_zero(y);
if(n==0)return;
if(n==1) m_copy(x_p[1],y);
else
{
Matrix temp;
m_pow_sum1(x_p,n>>1,temp);//计算 x^1+...+^(n/2)
m_add(temp,y,y);
//开始计算 x^(1+n/2)+...+x^n
Matrix temp2;
m_pow_with_saved(x_p,n>>1,temp2);
Matrix temp3;
m_multiple(temp,temp2,temp3);
m_add(temp3,y,y);
if(n&1)
{
m_multiple(temp2,temp2,temp3);//计算x^(n-1)
m_multiple(temp3,x_p[1],temp2);//计算 x^n
m_add(temp2,y,y);//add x^n
}
}
}
//计算x^1 x^2 x^4 ...,然后调用求幂和的地方
void m_pow_sum(Matrix x, unsigned int n, Matrix y)
{
//calculate x^1 x^2 x^4 ... x^logn
unsigned int i=0,logn,n2=n;
for(logn=0,n2=n;n2;logn++,n2 >>= 1);
Matrix *pow_arry=new Matrix[logn==0?2:(logn+1)];
m_one(pow_arry[0]);
m_copy(x,pow_arry[1]);
for(i=1;i<logn;i++)
{
m_multiple(pow_arry[i],pow_arry[i],pow_arry[i+1]);
}
m_pow_sum1(pow_arry,n,y);
delete []pow_arry;
}
int main()
{
//freopen("in.txt","r",stdin);
int k,b,n;
Matrix ans,res,sum;
while(scanf("%d%d%d%d",&k,&b,&n,&mod)!=EOF)
{
if(b)
m_pow(A,b,ans);//计算A^b
else
m_one(ans);
m_pow(A,k,res);//计算A^k
m_zero(sum);
m_pow_sum(res,n-1,sum);//计算A^(k*i)之和
m_multiple(ans,sum,res);//计算A^b*(sum(A^(k*i))),1<=i<=n
m_add(res,ans,sum);//加上i等于0的值
printf("%d\n",sum[0][1]%mod);
}
return 0;
}