求出斐波那契中的 第 k*i+b 项的和。

思路分析：

定义斐波那契数列的矩阵

f(n)为斐波那契第n项

F(n) = f(n+1)

   f(n)

那么可以知道矩阵

A = 1 1

       1  0

使得 F（n） = A * F（n+1）

然后我们化简最后的答案

sum = F（b） +   F（K+b） +  F （2*k +b）....

sum = F（b） +  A^k F(b)    +   A^2k F(b).....

sum = （E+A^k + A^2k.....）*F(b)

那么我们把 矩阵  A^k  定义为矩阵 K。

再递推上面的求和公式。

E E   *   SUM = SUM + E

0 K  E           K

所以构造一个内嵌矩阵的矩阵。

然后求出和乘以F（b）即可。

#include <cstdio>
#include <iostream>
#include <cstring>
#include <iostream>
#define N 2
using namespace std;
typedef long long LL;
LL mod;
struct matrix
{
    LL data[N][N];
    friend matrix operator * (const matrix A,const matrix B)
    {
        matrix res;
        memset(res.data,0,sizeof res.data);
        for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
        for(int k=0;k<N;k++)
        res.data[i][j]+=(A.data[i][k]*B.data[k][j])%mod;
        return res;
    }
    friend matrix operator + (const matrix A,const matrix B)
    {
        matrix res;
        for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
        res.data[i][j]=(A.data[i][j]+B.data[i][j])%mod;
        return res;
    }
    friend matrix operator - (const matrix A,const matrix B)
    {
        matrix res;
        for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
        res.data[i][j]=((A.data[i][j]-B.data[i][j])+mod)%mod;
        return res;
    }
    void print()
    {
        for(int i=0;i<N;i++)
        {
            for(int j=0;j<N;j++)
            printf("%I64d ",data[i][j]);
            puts("");
        }
    }

}E,zero;
struct supermax
{
    matrix ret[N][N];
    friend supermax operator * (supermax A,supermax B)
    {
        supermax res;
        for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
        res.ret[i][j]=zero;

        for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
        for(int k=0;k<N;k++)
        {
            res.ret[i][j]=res.ret[i][j]+(A.ret[i][k]*B.ret[k][j]);
            for(int p=0;p<N;p++)
            for(int q=0;q<N;q++)
            res.ret[i][j].data[p][q]%=mod;
        }
        return res;
    }
};

matrix matmod(matrix origin,LL n)
{
    matrix res=E;

    while(n)
    {
        if(n&1)
        res=res*origin;
        n>>=1;
        origin=origin*origin;
    }
    return res;
}
supermax Do(supermax origin,LL n)
{
    supermax res;
    for(int i=0;i<N;i++)
    for(int j=0;j<N;j++)
    res.ret[i][j]=zero;
    for(int i=0;i<N;i++)
    res.ret[i][i]=E;

    while(n)
    {
        if(n&1)
        res=res*origin;
        n>>=1;
        origin=origin*origin;
    }
    return res;
}
int main()
{
    memset(zero.data,0,sizeof zero.data);
    memset(E.data,0,sizeof E.data);
    for(int i=0;i<N;i++)
    E.data[i][i]=1;

    LL k,b,n;
    while(scanf("%I64d%I64d%I64d%I64d",&k,&b,&n,&mod)!=EOF)
    {
        matrix fib;
        fib.data[0][0]=1;
        fib.data[0][1]=1;
        fib.data[1][0]=1;
        fib.data[1][1]=0;

        matrix K=matmod(fib,k);

        supermax o;
        o.ret[0][0]=E;
        o.ret[0][1]=E;
        o.ret[1][0]=zero;
        o.ret[1][1]=K;

        supermax final=Do(o,n);

        matrix tmp=(final.ret[0][0]*zero)+(final.ret[0][1]*E);

        matrix B=matmod(fib,b);
        matrix fibb,fib0;
        fib0.data[0][0]=1;
        fib0.data[1][0]=0;
        fib0.data[0][1]=fib0.data[1][1]=0;
        fibb=B*fib0;

        matrix ans = tmp*fibb;
        printf("%I64d\n",ans.data[1][0]%mod);
    }
    return 0;
}