F(g(i)) = F(b) + F(b+k)+F(b+2k)+....+F(b+nk)

           = F(b) + (A^k)F(b) + (A^2k)F(b)+….+(A^nk)F(b)

           = F(b) [ E +A^k + A^2k + ….+ A^（nk）]   (E為單位矩陣)

令 K = A^k 

   原式=E +A^k + A^2k + ….+ A^nk 變成 K^0+K^1+K^2+…+K^n

構造矩陣

              （圖源網路）

#include<iostream>
#include<string>
#include<cstdio>
#include<cstring>
#include<queue>
#include<map>
#include<cmath>
#include<stack>
#include<set>
#include<vector>
#include<algorithm>
#define ll long long
#define inf 1<<30
using namespace std;
int MAXN;
int mod;
struct Matrix
{
	ll mat[15][15];
	Matrix() {}
	Matrix operator*(Matrix const &b)const
	{
		Matrix res;
		memset(res.mat, 0, sizeof(res.mat));
		for (int i = 1 ;i <=MAXN; i++)
			for (int j = 1; j <= MAXN; j++)
				for (int k = 1; k <=MAXN; k++)
					res.mat[i][j] = (res.mat[i][j]+this->mat[i][k] * b.mat[k][j])%mod;
		return res;
	}
};
Matrix base,ans1,ans2,ans;
Matrix pow_mod(Matrix base, int n)
{
	Matrix res;
	memset(res.mat, 0, sizeof(res.mat));
	for (int i = 1; i <=MAXN; i++)
		res.mat[i][i] = 1;
	while (n > 0)
	{
		if (n & 1) res = res*base;
		base = base*base;
		n >>= 1;
	}
	return res;
}

int main()
{
    int n,k,m,t,b;
    while(~scanf("%d%d%d%d",&k,&b,&n,&m))
    {
        mod=m;
        memset(base.mat,0,sizeof(base.mat));
        memset(ans1.mat,0,sizeof(ans1.mat));
        base.mat[1][1]=1;
        base.mat[1][2]=1;
        base.mat[2][1]=1;
        base.mat[2][2]=0;
        MAXN=2;
        ans1=pow_mod(base,k);//K
        ans2=pow_mod(base,b);
        ans1.mat[1][3]=1;
        ans1.mat[2][4]=1;
        ans1.mat[3][3]=1;
        ans1.mat[4][4]=1;
        MAXN=4;
        ll sum=0;
        Matrix ans3=pow_mod(ans1,n);
        for(int i=1;i<=2;i++)
            sum=(sum+(ans2.mat[2][i]*ans3.mat[i][3])%mod)%mod;
        printf("%lld\n",sum);

    }
}