hdu 1588 求f(b) +f(k+b) +f(2k+b) +f((n-1)k +b) 之和 （矩阵快速幂）

g(i)=k*i+b; 0<=i<n
f(0)=0
f(1)=1
f(n)=f(n-1)+f(n-2) (n>=2)
求f(b) +f(k+b) +f(2*k+b) +f((n-1)*k +b) 之和

Sample Input
2 1 4 100 // k b n MOD
2 0 4 100

Sample Output
21
12

矩阵A      相当于

1 1          f（2）  f(1)

1 0          f（1）  f(0)

                     | 1       1| ^b          | f(b+1)    f(b)|  

         mat^b =|1       0 |          =  | f(b)    f(b-1)|  

求f（n） 就是求矩阵A的n次幂 再取第1行第2列的元素

要求的东西可化成  A^b*( I + A^k + (A^k)^2 + .... + (A^k)^(N-1) )  

矩阵ans1 = A^b

矩阵B = A^k

矩阵C =

B  I

O  I

C的n次幂后  再取右上的小矩阵  就是I+B+B^2....+B^(n-1)  赋给ans2

ans1 * ans  再取第1行第2列的元素 就是最终答案

# include <iostream>
# include <cstdio>
# include <cstring>
# include <algorithm>
# include <cmath>
# define LL long long
using namespace std ;

LL MOD ;

struct Matrix
{
    LL mat[6][6];
};

Matrix mul(Matrix a,Matrix b, int n)
{
    Matrix c;
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
        {
            c.mat[i][j]=0;
            for(int k=0;k<n;k++)
            {
                c.mat[i][j]=(c.mat[i][j] + a.mat[i][k]*b.mat[k][j])%MOD;
            }
        }
    return c;
}
Matrix pow_M(Matrix a,int k , int n)
{
    Matrix ans;
    memset(ans.mat,0,sizeof(ans.mat));
    for (int i=0;i<n;i++)
        ans.mat[i][i]=1;
    Matrix temp=a;
    while(k)
    {
        if(k&1)ans=mul(ans,temp,n);
        temp=mul(temp,temp,n);
        k>>=1;
    }
    return ans;
}



int main ()
{
    //freopen("in.txt","r",stdin) ;
    Matrix A ;
    A.mat[0][0] = A.mat[0][1] = A.mat[1][0] = 1 ;
    A.mat[1][1] = 0 ;
    int n , k , b ;
    while(cin>>k>>b>>n>>MOD)
    {
        int i ,j ;
        Matrix ans1 , ans2;
        ans1 = pow_M(A,b,2) ;
    
        Matrix B , C ;
        B = pow_M(A,k,2) ;
        memset(C.mat,0,sizeof(C.mat));
        for (i = 0 ; i < 2 ; i++) //扩展成4 * 4的矩阵C
        {
            for (j = 0 ; j < 2 ; j++)
            {
                C.mat[i][j] = B.mat[i][j] ;
            }
            C.mat[2+i][2+i] = 1 ;
            C.mat[i][2+i] = 1 ;
        }
        ans2 = pow_M(C,n,4) ; // 4*4
        
        ans2.mat[0][0] = ans2.mat[0][2] ;
        ans2.mat[0][1] = ans2.mat[0][3] ;
        ans2.mat[1][0] = ans2.mat[1][2] ;
        ans2.mat[1][1] = ans2.mat[1][3] ;
        ans1 = mul(ans1,ans2,2) ;

        cout<<ans1.mat[0][1]%MOD<<endl ;
    }

    return 0 ;
}