hdu1588 矩阵快速幂

//看了很多的博客 后来队友指点才懂
//sum=f(g(0))+f(g(1))+....
//sum=A^(b-1)*|...|....
//要将b-1换，防止出现b=0时有负一，用A^b代替，取下面的即可
//这样问题成了 sum=A^b(A+A^(2k)+A^(3k)+...+A^(k(n-1)));
//令B=A^k次，就简单了。
/*
主要要求1+A+A^2+A^3+...+A^(n-1)次方
|  A   A  |        | A^2  A^2+A |                        |  A^(k-1)   A^(k-1)+A^(k-2)+A^(k-3)... |
令B=  |         |    B^2=|            |  这样可以得到  B^(n-1)=|                                       |
|  0   1  |        | 0      1   |                        |     0                1                |
*/
#include<stdio.h>
#include<string.h>
#define maxn 30
#define ll __int64
ll n,mod;
struct Mat
{
ll mat[maxn][maxn];
};
Mat cal1(Mat a,Mat b,int nn)//矩阵乘法
{
Mat c;
memset(c.mat,0,sizeof(c.mat));
int i,j,k;
for(i=0;i<nn;i++)
for(j=0;j<nn;j++)
for(k=0;k<nn;k++)
{
c.mat[i][j]+=((a.mat[i][k]*b.mat[k][j])%mod);
c.mat[i][j]%=mod;
}
return c;
}
Mat cal2(Mat a,ll k,int nn)//矩阵幂
{
int i,j;
Mat c;
for(i=0;i<nn;i++)
for(j=0;j<nn;j++)
if(i==j)c.mat[i][j]=1;
else c.mat[i][j]=0;
while(k)
{
if(k&1)
c=cal1(c,a,nn);
k=k>>1;
a=cal1(a,a,nn);
}
return c;
}
Mat A,B,S;
void initA()
{
A.mat[0][0]=1;
A.mat[0][1]=1;
A.mat[1][0]=1;
A.mat[1][1]=0;
}
void initB()
{
int i,j;
for(i=0;i<2;i++)
for(j=0;j<2;j++)
{
B.mat[i][j]=A.mat[i][j];
B.mat[i][j+2]=A.mat[i][j];
}
for(i=2;i<4;i++)
for(j=0;j<2;j++)
B.mat[i][j]=0;
for(i=2;i<4;i++)
for(j=2;j<4;j++)
if(i==j)
B.mat[i][j]=1;
else B.mat[i][j]=0;
}
Mat getmat()
{
int i,j;
Mat c;
for(i=0;i<2;i++)
for(j=2;j<4;j++)
c.mat[i][j-2]=B.mat[i][j];
for(i=0;i<2;i++)
for(j=0;j<2;j++)
if(i==j)
c.mat[i][j]+=1;
return c;
}
int main()
{
ll i,j,k,b;
while(scanf("%I64d %I64d %I64d %I64d",&k,&b,&n,&mod)!=EOF)
{
initA();
S=cal2(A,b,2);
A=cal2(A,k,2);
initB();
B=cal2(B,n-1,4);
Mat temp=getmat();
S=cal1(S,temp,2);
/*for(i=0;i<2;i++)
{
for(j=0;j<2;j++)
printf("%I64d ",S.mat[i][j]);
printf("\n");
}*/
printf("%d\n",S.mat[1][0]);

}
}