思路:显然的有f(n)=A^n,那么把g(i)代入后有sum = A^b+A^(k+b)+A^(2k+b)+...A^((n-1)k+b)
提取A^b之后sum = A^b*(A^k^0+A^k^1+A^k^2+.....A^k^(n-1)),将A^k设为B的话那么就是一个简单的二分求和过程,注意最后要加上B^0
坑点:虽然最后结果不会爆long long,可是中途乘法过程中是会爆的
#include<iostream>
#include<cstdio>
#include<cstring>
using namespace std;
#define LL long long
int n,mod,k,b;
struct Mat
{
LL a[31][31]; //矩阵大小
}pi;
Mat mul(Mat a,Mat b)
{
Mat t;
memset(t.a,0,sizeof(t.a));
for(int i = 0;i<2;i++)
{
for(int k = 0;k<2;k++)
{
if(a.a[i][k])
for(int j = 0;j<2;j++)
{
t.a[i][j]+=a.a[i][k]*b.a[k][j];
if(t.a[i][j]>=mod)
t.a[i][j]%=mod;
}
}
}
return t;
}
Mat expo(Mat p,int k)
{
if(k==1)return p;
Mat e;
memset(e.a,0,sizeof(e.a));
for(int i = 0;i<2;i++) //初始化单位矩阵
e.a[i][i]=1;
if(k==0)return e;
while(k)
{
if(k&1)
e = mul(p,e);
p = mul(p,p);
k>>=1;
}
return e;
}
Mat add(Mat a,Mat b)
{
Mat t;
memset(t.a,0,sizeof(t.a));
for(int i = 0;i<2;i++)
for(int j = 0;j<2;j++)
t.a[i][j]=(a.a[i][j]+b.a[i][j])%mod;
return t;
}
Mat sum(int k)
{
if(k==1)
return pi;
if(k&1)
return add(sum(k-1),expo(pi,k));
else
{
Mat tmp = expo(pi,k/2);
Mat temp = sum(k/2);
return add(temp,mul(tmp,temp));
}
}
int main()
{
while(scanf("%d%d%d%d",&k,&b,&n,&mod)!=EOF)
{
memset(pi.a,0,sizeof(pi.a));
Mat p;
p.a[0][0]=p.a[0][1]=p.a[1][0]=1;
p.a[1][1]=0;
pi = expo(p,k); //A^k
Mat res = expo(p,b); //A^b
Mat tmp = sum(n-1);
Mat e;
memset(e.a,0,sizeof(e.a));
for(int i = 0;i<2;i++) //初始化单位矩阵
e.a[i][i]=1;
tmp = add(tmp,e);
tmp = mul(tmp,res);
printf("%d\n",tmp.a[0][1]);
}
}