sum=f(b)+f(k*1+b)+f(k*2+b)+...+f(k*n+b)=(P^(b-1)+P^(k*1+b-1)+...+P^(k*n+b-1))(f(1) f(0));把P^(b-1)提出去,
#include <cstdio>
#include <iostream>using namespace std;
int mod;
typedef struct
{
long long m[2][2];
}matrix;
typedef struct
{
long long m[4][4];
}matrix1;
matrix I={1,0,0,1};
matrix P={1,1,1,0};
matrix O={0,0,0,0};
matrix AA={0,1,1,-1}; //AA为P的逆矩阵,因为考虑到b=0时,P^(b-1)就变为P的-1次方,正好是P的逆矩阵
matrix1 I1={1,0,0,0,
0,1,0,0,
0,0,1,0,
0,0,0,1};
matrix1 P1={0,0,1,0,
0,0,0,1,
0,0,1,0,
0,0,0,1};
matrix mul(matrix a,matrix b)
{
int i,j,k;
matrix c;
for(i=0;i<2;i++)
for(j=0;j<2;j++)
{
c.m[i][j]=0;
for(k=0;k<2;k++)
c.m[i][j]+=((a.m[i][k]%mod)*(b.m[k][j]%mod))%mod;
c.m[i][j]%=mod;
}
return c;
}
matrix quick_mod(int n)
{
matrix a=P,b=I;
while(n>0)
{
if(n&1)
b=mul(b,a);
n=n>>1;
a=mul(a,a);
}
return b;
}
matrix1 mul1(matrix1 a,matrix1 b)
{
int i,j,k;
matrix1 c;
for(i=0;i<4;i++)
for(j=0;j<4;j++)
{
c.m[i][j]=0;
for(k=0;k<4;k++)
c.m[i][j]+=((a.m[i][k]%mod)*(b.m[k][j]%mod))%mod;
c.m[i][j]%=mod;
}
return c;
}
matrix1 quick_mod1(int n)
{
matrix1 a=P1,b=I1;
while(n>0)
{
if(n&1)
b=mul1(b,a);
n=n>>1;
a=mul1(a,a);
}
return b;
}
int main()
{
int k,b,n;
while(scanf("%d%d%d%d",&k,&b,&n,&mod)==4)
{
matrix B=quick_mod(k);
P1.m[0][0]=B.m[0][0];
P1.m[0][1]=B.m[0][1];
P1.m[1][0]=B.m[1][0];
P1.m[1][1]=B.m[1][1];
matrix1 R=quick_mod1(n);
matrix Q,A;
Q.m[0][0]=R.m[0][2];
Q.m[0][1]=R.m[0][3];
Q.m[1][0]=R.m[1][2];
Q.m[1][1]=R.m[1][3];
if(b==0) A=AA;
else A=quick_mod(b-1);
Q=mul(Q,A);
printf("%d\n",(Q.m[0][0]+mod)%mod);
}
return 0;
}