http://www.cnblogs.com/ACShiryu/archive/2011/08/09/poj3233.html
1 /* 这道题目借鉴他人的思路和代码，很有收获
2 先看下面这个快速幂求余的运算
3 递归用二分法，每个过程都求余。
4 long exp_mod(long a,long n,long b)
5 {
6 long t;
7 if(n==0) return 1%b;
8 if(n==1) return a%b;
9 t=exp_mod(a,n/2,b);
10 t=t*t%b;
11 if((n&1)==1) t=t*a%b;
12 return t;
13 }
14 本题用了两次二分的思想，一次用于求A^K,第二次用于求 (1+A^(k/2))*(A + A^2 + A^3 + … + A^(k/2))
15 */
16
17
18 #include<iostream>
19 #include<stdio.h>
20 using namespace std;
21 //=======================================================
22 int n,m,k;
23 struct Point//定义结构体更加方便,用数组时结构体有大作用的
24 {
25 int a[32][32];
26 void init()
27 {
28 memset(a,0,sizeof(a));
29 for(int i=1;i<=32;++i)
30 {
31 a[i][i]=1;
32 }
33 }
34 };
35 //=======================================================
36 Point mul(Point x,Point y)//乘法运算
37 {
38 Point z;
39 for(int i=1;i<=n;++i)
40 {
41 for(int j=1;j<=n;++j)
42 {
43 int sum=0;
44 for(int k=1;k<=n;++k)
45 {
46 sum+=((x.a[i][k]*y.a[k][j])%m);//每次都求余
47 z.a[i][j]=(sum%m);//这时也要求余。。。。。。
48 }
49 }
50 }
51 return z;
52 }
53
54 Point add(Point p,Point q)//矩阵加法运算
55 {
56 Point c;
57 for(int i=1;i<=n;++i)
58 {
59 for(int j=1;j<=n;++j)
60 {
61 c.a[i][j]=((p.a[i][j]+q.a[i][j])%m);//记得求余
62 }
63 }
64 return c;
65 }
66
67 void display(Point p)//打印函数
68 {
69 int i,j;
70 for( i=1;i<=n;++i)
71 {
72 for(j=1;j<=n-1;++j)
73 {
74 cout<<p.a[i][j]<<" ";//注意格式“ ”问题
75 }
76 if(j==n)cout<<p.a[i][j]<<endl;
77
78 }
79 }
80
81 Point maxtrimul(Point s,int k)//用递归定义实现快速幂的乘法
82 {
83 Point ans;//用来存放结果的
84 ans.init();
85 if(k==1)return s;
86 ans=maxtrimul(s,k>>1);//每次二分
87 ans=mul(ans,ans);
88 if(k&1) //奇数时再乘一次
89 {
90 ans=mul(s,ans);
91 }
92 return ans;
93
94 }
95
96 Point maxtrisum(Point s,int k)
97 {
98 if(k==1)return s;
99 Point ans;//用来保存答案
100 ans.init();
101 ans=add(ans,maxtrimul(s,k>>1));
102 ans=mul(ans,maxtrisum(s,k>>1));
103 if(k&1)//奇数时
104 {
105 ans=add(ans,maxtrimul(s,k));
106 }
107 return ans;
108
109 }
110
111 int main()
112 {
113 cin>>n>>k>>m;
114 Point s;
115 for(int i=1;i<=n;++i)
116 {
117 for(int j=1;j<=n;++j)
118 {
119 cin>>s.a[i][j];
120 }
121 }
122 s=maxtrisum(s,k);
123 display(s);
124
125 }
126