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列的元素 就是最终答案

 1 # include <iostream>
 2 # include <cstdio>
 3 # include <cstring>
 4 # include <algorithm>
 5 # include <cmath>
 6 # define LL long long
 7 using namespace std ;
 8
 9 LL MOD ;
10
11 struct Matrix
12 {
13     LL mat[6][6];
14 };
15
16 Matrix mul(Matrix a,Matrix b, int n)
17 {
18     Matrix c;
19     for(int i=0;i<n;i++)
20         for(int j=0;j<n;j++)
21         {
22             c.mat[i][j]=0;
23             for(int k=0;k<n;k++)
24             {
25                 c.mat[i][j]=(c.mat[i][j] + a.mat[i][k]*b.mat[k][j])%MOD;
26             }
27         }
28     return c;
29 }
30 Matrix pow_M(Matrix a,int k , int n)
31 {
32     Matrix ans;
33     memset(ans.mat,0,sizeof(ans.mat));
34     for (int i=0;i<n;i++)
35         ans.mat[i][i]=1;
36     Matrix temp=a;
37     while(k)
38     {
39         if(k&1)ans=mul(ans,temp,n);
40         temp=mul(temp,temp,n);
41         k>>=1;
42     }
43     return ans;
44 }
45
46
47
48 int main ()
49 {
50     //freopen("in.txt","r",stdin) ;
51     Matrix A ;
52     A.mat[0][0] = A.mat[0][1] = A.mat[1][0] = 1 ;
53     A.mat[1][1] = 0 ;
54     int n , k , b ;
55     while(cin>>k>>b>>n>>MOD)
56     {
57         int i ,j ;
58         Matrix ans1 , ans2;
59         ans1 = pow_M(A,b,2) ;
60
61         Matrix B , C ;
62         B = pow_M(A,k,2) ;
63         memset(C.mat,0,sizeof(C.mat));
64         for (i = 0 ; i < 2 ; i++) //扩展成4 * 4的矩阵C
65         {
66             for (j = 0 ; j < 2 ; j++)
67             {
68                 C.mat[i][j] = B.mat[i][j] ;
69             }
70             C.mat[2+i][2+i] = 1 ;
71             C.mat[i][2+i] = 1 ;
72         }
73         ans2 = pow_M(C,n,4) ; // 4*4
74
75         ans2.mat[0][0] = ans2.mat[0][2] ;
76         ans2.mat[0][1] = ans2.mat[0][3] ;
77         ans2.mat[1][0] = ans2.mat[1][2] ;
78         ans2.mat[1][1] = ans2.mat[1][3] ;
79         ans1 = mul(ans1,ans2,2) ;
80
81         cout<<ans1.mat[0][1]%MOD<<endl ;
82     }
83
84     return 0 ;
85 }