http://acm.hdu.edu.cn/showproblem.php?pid=1588
/*
把斐波那契数列转化为矩阵:
A={1 1}
{1,0};
{f[n+1],f[n]}
{f[n],f[n-1]} = A^n ;最后输出右上角那项
或者用
{f[n+2],f[n+1]}
{f[n+1], f[n] } = A^(n+1); 最后输出右下角那项
我们用第一个公式
所求即为A^b + A^(k+b) + A^(2*k+b) + ... + A^((n-1)*k+b)
=A^b * ( A^0 + A^k + A^(2*k) + ... + A^((n-1)*k) )
=A^b * ( (A^k)^0 + (A^k)^1 + (A^k)^2 + ...+ (A^k)^(n-1) );
B=A^k;
上式=A^b * ( B^0 + B^1 + B^2 + ... + B^(n-1) );
B^0 + B^1 + B^2 + ... + B^(n-1)用上篇介绍到的递归二分 方法求解
最后输出矩阵的第二项(右上角)即可；
对于求解 B^0 + B^1 + B^2 + ... + B^(n-1)
我们也可以构造矩阵的方法。
我们来设置这样一个矩阵
B I
O I
其中O是零矩阵，I是单位矩阵
将它乘方，得到
B^2 I+B
O I
乘三方，得到
B^3 I+B+B^2
O I
乘四方，得到
B^4 I+B+B^2+B^3
O I
用快速幂求出n方，让第二项再与A^b相乘即可。
*/
View Code
1 # include<stdio.h>
2 # include<string.h>
3 int m;
4 struct matrix{
5 __int64 a[3][3];
6 };
7 matrix ans2;
8 matrix power(matrix un,matrix un1)
9 {
10 matrix ans;
11 int i,j,h;
12 for(i=1;i<=2;i++)
13 for(j=1;j<=2;j++)
14 {
15 ans.a[i][j]=0;
16 for(h=1;h<=2;h++)
17 {
18 ans.a[i][j]+=un.a[i][h]*un1.a[h][j];
19 ans.a[i][j]%=m;
20 }
21 }
22 return ans;
23 }
24 matrix mod(matrix un,int k)//快速幂
25 {
26 matrix ans;
27 ans.a[1][1]=1;
28 ans.a[1][2]=0;
29 ans.a[2][1]=0;
30 ans.a[2][2]=1;
31 while(k)
32 {
33 if(k%2) ans=power(ans,un);
34 un=power(un,un);
35 k/=2;
36 }
37 return ans;
38 }
39 matrix dfs(int n)//递归二分
40 {
41 matrix ans1,ans;
42 int i,j;
43 if(n==0)
44 {
45 for(i=1;i<=2;i++)
46 for(j=1;j<=2;j++)
47 ans.a[i][j]=0;
48 return ans;
49 }
50 if(n==1) return ans2;
51 ans1=mod(ans2,n/2);
52 ans1.a[1][1]++;
53 ans1.a[2][2]++;//(1 + A^(n/2) )
54 ans=dfs(n/2);//（ A^1 + A^2 + ... + A^(n/2)）;
55 ans=power(ans1,ans);
56 if(n%2)//判奇偶
57 {
58 ans1=mod(ans2,n);
59 for(i=1;i<=2;i++)
60 for(j=1;j<=2;j++)
61 {
62 ans.a[i][j]+=ans1.a[i][j];
63 ans.a[i][j]%=m;
64 }
65 }
66 return ans;
67 }
68 int main()
69 {
70 int i,j,k,b,n;
71 matrix ans,un,ans1,mm;
72 while(scanf("%d%d%d%d",&k,&b,&n,&m)!=EOF)
73 {
74 un.a[1][1]=1;
75 un.a[1][2]=1;
76 un.a[2][1]=1;
77 un.a[2][2]=0;
78 ans=mod(un,b+1);//把f[n]放在右下角
79 //ans=mod(un,b+1);把f[n]放在右上角或者左下角
80
81 ans1.a[1][1]=1;
82 ans1.a[1][2]=0;
83 ans1.a[2][1]=0;
84 ans1.a[2][2]=1;//单位矩阵
85
86 ans2=mod(un,k);
87
88 n--;
89 mm=dfs(n);
90 for(i=1;i<=2;i++)
91 for(j=1;j<=2;j++)
92 {
93 mm.a[i][j]+=ans1.a[i][j];
94 mm.a[i][j]%=m;
95 }
96 mm=power(ans,mm);
97 printf("%d\n",mm.a[2][2]%m);//f[n]在右下角
98 //f[n]在右上角或者左下角时输出mm.a[1][2]%m
99 }
100 return 0;
101 }
