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 }