矩阵运用

十个利用矩阵乘法解决的经典题目：

http://www.matrix67.com/blog/archives/276

1.hodj 1575 TrA

求矩阵的k次方

 

代码 
 1 #include <iostream>
 2 #include <cstdio>
 3  using namespace std;
 4 
 5  struct Matrix
 6 {
 7     int matrix[12][12];
 8 }unit, init, res;
 9 
10  int n;
11 Matrix Mul(Matrix a, Matrix b)
12 {
13     Matrix c;
14     int i, j, k;
15     for(i=0; i<n; ++i){
16         for(j=0; j<n; ++j){
17             int temp = 0;
18             for(k=0; k<n; ++k){
19                 temp += (a.matrix[i][k] * b.matrix[k][j]) % 9973;
20                 temp %= 9973;
21             }
22             c.matrix[i][j] = temp;
23         }
24     }
25     return c;
26 }
27 Matrix Cal(int k)
28 {
29     Matrix p, q;
30     p = unit;
31     q = init;
32     while(k != 1){
33         if(k&1){
34             p = Mul(p, q);
35             k --;
36         }else {
37             q = Mul(q, q);
38             k >>= 1;
39         }
40     }
41     p = Mul(p, q);
42     return p;
43 }
44  int main()
45 {
46  //    freopen("c:/aaa.txt", "r", stdin);
47      int T, i, j, k;
48     scanf("%d", &T);
49     while(T--){
50         scanf("%d %d", &n, &k);
51         for(i=0; i<n; ++i){
52             for(j=0; j<n; ++j){
53                 scanf("%d", &init.matrix[i][j]);
54                 unit.matrix[i][j] = (i == j);
55             }
56         }
57         res = Cal(k);
58         int ans = res.matrix[0][0];
59         for(i=1; i<n; ++i){
60             ans = (ans + res.matrix[i][i]) % 9973;
61         }
62         printf("%d\n", ans);
63     }
64     return 0;
65 }
66  

 

2.hdoj 1757 A Simple Math Problem

递推的矩阵求法。注意初始化，如 f(n) = 4f(n-1) – 3f(n-2) + 2f(n-4)

第k项为：

 

代码 
 1 #include <iostream>
 2 #include <cstdio>
 3  using namespace std;
 4 
 5  struct Matrix
 6 {
 7     int matrix[12][12];
 8 }unit, init, res;
 9 
10  int m;
11 Matrix Mul(Matrix a, Matrix b)
12 {
13     Matrix c;
14     int i, j, k;
15     for(i=0; i<10; ++i){
16         for(j=0; j<10; ++j){
17             int temp = 0;
18             for(k=0; k<10; ++k){
19                 temp += (a.matrix[i][k] * b.matrix[k][j]) % m;
20                 temp %= m;
21             }
22             c.matrix[i][j] = temp;
23         }
24     }
25     return c;
26 }
27 Matrix Cal(int k)
28 {
29     Matrix p, q;
30     p = unit;
31     q = init;
32     while(k != 1){
33         if(k&1){
34             p = Mul(p, q);
35             k --;
36         }else {
37             q = Mul(q, q);
38             k >>= 1;
39         }
40     }
41     p = Mul(p, q);
42     return p;
43 }
44  int main()
45 {
46  //    freopen("c:/aaa.txt", "r", stdin);
47      int i, j, k;   
48     while(scanf("%d %d", &k, &m) == 2){
49         if(k < 10){
50             printf("%d\n", k%m);
51             continue;
52         }
53         for(i=0; i<10; ++i){
54             for(j=0; j<10; ++j){
55                 init.matrix[i][j] = 0;
56                 unit.matrix[i][j] = (i == j);
57             }
58         }
59         for(i=0; i<9; ++i) init.matrix[i][1+i] = 1;
60         for(i=0; i<10; ++i) scanf("%d", &init.matrix[9][9-i]);
61         k -= 9;
62         res = Cal(k);
63         int temp = 0;
64         for(i=0; i<10; ++i){
65             temp += (res.matrix[9][i] * i) % m;
66             temp %= m;
67         }
68         printf("%d\n", temp);
69     }
70     return 0;
71 }

3.hdoj 2254 奥运

即求求A + A^2 + A^3 + … + A^k。二分。

 

代码 
  1 #include <iostream>
  2 #include <cstdio>
  3 #include <set>
  4 #include <iterator>
  5 #include <algorithm>
  6 using namespace std;
  7 
  8 typedef __int64 lld;
  9 const int maxn = 30;
 10 struct Matrix
 11 {
 12     int matrix[maxn][maxn];
 13 }unit, init;
 14 
 15 struct Node
 16 {
 17     lld a, b;
 18 }ns[10010];
 19 
 20 int n, idx;
 21 lld num[35];
 22 Matrix Add(Matrix A, Matrix B)
 23 {
 24     int i, j;
 25     Matrix C;
 26     for(i=0; i<idx; ++i){
 27         for(j=0; j<idx; ++j){
 28             C.matrix[i][j] = (A.matrix[i][j] + B.matrix[i][j]) % 2008;
 29         }
 30     }
 31     return C;
 32 }
 33 Matrix Mul(Matrix A, Matrix B)
 34 {
 35     int i, j, k, temp;
 36     Matrix C;
 37     for(i=0; i<idx; ++i){
 38         for(j=0; j<idx; ++j){
 39             temp = 0;
 40             for(k=0; k<idx; ++k){
 41                 temp += (A.matrix[i][k] * B.matrix[k][j]) % 2008;
 42                 temp %= 2008;
 43             }
 44             C.matrix[i][j] = temp;
 45         }
 46     }
 47     return C;
 48 }
 49 Matrix Pow(Matrix A, int k)
 50 {
 51     Matrix p, q;
 52     p = unit;
 53     q = init;
 54     while(k != 1){
 55         if(k & 1){
 56             p = Mul(p, q);
 57             k --;
 58         }else{
 59             q = Mul(q, q);
 60             k >>= 1;
 61         }
 62     }
 63     p = Mul(p, q);
 64     return p;
 65 }
 66 Matrix Sum(Matrix A, int k)
 67 {
 68     Matrix B, C;
 69     if(k == 1) return A;
 70     B = Sum(A, k>>1);
 71     if(k & 1){
 72         C = Pow(A, (k>>1)+1);   
 73         return Add(Add(B, C), Mul(C, B));
 74     }else {
 75         C = Pow(A, k>>1);
 76         return Add(B, Mul(C, B));
 77     }
 78 }
 79    
 80 
 81 int find(lld x)
 82 {
 83     int l, r, mid;
 84     l = 0;
 85     r = idx - 1;
 86     while(l <= r){
 87         mid = (l + r) >> 1;
 88         if(num[mid] == x) return mid;
 89         else if(num[mid] < x){
 90             l = mid + 1;
 91         }else {
 92             r = mid - 1;
 93         }
 94     }
 95     return -1;
 96 }
 97 
 98 void solve()
 99 {
100     int t1, t2, i, j, k;
101     lld v1, v2;
102     scanf("%d", &n);
103     for(i=0; i<n; ++i){
104         scanf("%I64d %I64d %d %d", &v1, &v2, &t1, &t2);
105         if(0 == t1 && 0 == t2){
106             printf("0\n");
107             continue;
108         }
109         int finda = find(v1);
110         int findb = find(v2);
111         if(finda == -1 || findb == -1){
112             printf("0\n");
113             continue;
114         }
115         if(t1 > t2){
116             int sw;
117             sw = t1;
118             t1 = t2;
119             t2 = sw;
120         }
121         int temp;
122         Matrix p, q;
123         if(t1 == 0 || t1 == 1){
124             q = Sum(init, t2);
125             temp = q.matrix[finda][findb];
126         }else{
127             p = Sum(init, t1 - 1);
128             q = Sum(init, t2);
129             temp = (q.matrix[finda][findb] - p.matrix[finda][findb]) % 2008;
130         }
131         while(temp < 0) temp += 2008;
132         temp %= 2008;
133         printf("%d\n", temp);
134     }
135 }
136 void Init()
137 {
138     int i, j;
139     set<lld> ss;
140     set<lld>::iterator ite;
141     ss.clear();
142     for(i=0; i<n; ++i){
143         scanf("%I64d %I64d", &ns[i].a, &ns[i].b);
144         ss.insert(ns[i].a);
145         ss.insert(ns[i].b);
146     }
147     idx = 0;
148     for(ite=ss.begin(); ite!=ss.end(); ++ite)
149         num[idx++] = *ite;
150 
151     for(i=0; i<idx; ++i)
152         for(j=0; j<idx; ++j){
153             unit.matrix[i][j] = (i == j);
154             init.matrix[i][j] = 0;
155         }
156    
157            
158     for(i=0; i<n; ++i){
159         int finda, findb;
160         finda = find(ns[i].a);
161         findb = find(ns[i].b);
162         init.matrix[finda][findb] ++;
163         init.matrix[finda][findb] %= 2008;
164     }
165 }
166    
167 
168 int main()
169 {
170  //   freopen("c:/aaa.txt", "r", stdin);   
171     while(scanf("%d", &n) == 1){
172         Init();
173         solve();
174     }
175     return 0;
176 }

4.hdoj 2971 Tower

An = 2 * a2 * An-1 – An-2    ;   Sn = A1^2 + A2^2 + A3^2 + …+An^2

Sn = An^2 + Sn-1  ==> Sn – Sn-1 = An^2 ;

所以A的平方项都可以用S代替，Sn = An^2 + Sn-1  ==> Sn = 4*A2^2*Sn-1 + (2-8*A2^2) * Sn-2 + 4*A2^2*Sn-3 – Sn-4;

 

代码 
 1 #include <iostream>
 2 #include <cstdio>
 3 #include <cstring>
 4 using namespace std;
 5 
 6 typedef __int64 lld;
 7 const int maxn = 4;
 8 struct Matrix
 9 {
10     lld matrix[maxn][maxn];
11 }mat, unit;
12 lld t, n, mod;
13 lld s1, s2, s3, s4;
14 Matrix Mul(Matrix &a, Matrix &b)
15 {
16     Matrix c;
17     int i, j, k;
18     for(i=0; i<4; ++i){
19         for(j=0; j<4; ++j){
20             c.matrix[i][j] = 0;
21             for(k=0; k<4; ++k){
22                 if(a.matrix[i][k] && b.matrix[k][j])
23                     c.matrix[i][j] = (c.matrix[i][j] + a.matrix[i][k] * b.matrix[k][j]) % mod;
24             }
25         }
26     }
27     return c;
28 }
29 
30 Matrix Cal(int k)
31 {
32     Matrix p, q;
33     p = unit;
34     q = mat;
35     while(k){
36         if(k & 1) p = Mul(p, q);
37         q = Mul(q, q);
38         k >>= 1;
39     }
40     return p;
41 }
42 lld solve()
43 {
44     int i, j;
45     if(n == 1) return s1 % mod;
46     if(n == 2) return s2 % mod;
47     if(n == 3) return s3 % mod;
48     if(n == 4) return s4 % mod;
49     Matrix ans = Cal(n - 4);
50     lld temp;
51     temp = (ans.matrix[3][0] * s1 % mod + ans.matrix[3][1] * s2 % mod +
52             ans.matrix[3][2] * s3 % mod + ans.matrix[3][3] * s4 % mod) % mod;
53     if(temp < 0) temp = (temp + mod) % mod;
54     return temp;
55 }
56 void init()
57 {
58     int i, j;
59     lld a1, a2, a3, a4;
60     a1 = 1;
61     a2 = t % mod;
62     a3 = (2 * t * a2 - a1) % mod;
63     a4 = (2 * t * a3 - a2) % mod;
64     s1 = 1;
65     s2 = (s1 + a2 * a2) % mod;
66     s3 = (s2 + a3 * a3) % mod;
67     s4 = (s3 + a4 * a4) % mod;
68     for(i=0; i<4; ++i){
69         for(j=0; j<4; ++j){
70             mat.matrix[i][j] = 0;
71             unit.matrix[i][j] = (i == j);
72         }
73     }
74     mat.matrix[0][1] = mat.matrix[1][2] = mat.matrix[2][3] = 1;
75     mat.matrix[3][0] = -1;
76     mat.matrix[3][1] = (4 * t * t) % mod;
77     mat.matrix[3][2] = (2 - 8 * t * t) % mod;
78     mat.matrix[3][3] = (4 * t * t) % mod;
79 }
80 
81 int main()
82 {
83  //   freopen("c:/aaa.txt", "r", stdin);
84     int T;
85     scanf("%d", &T);
86     while(T--){
87         scanf("%I64d %I64d %I64d", &t, &n, &mod);
88         init();
89         printf("%I64d\n", solve());
90     }
91     return 0;
92 }

5.hdoj 2294 Pendant

dp[i][j] = dp[i-1][j] * j + dp[i-1][j-1] * (k-j+1);

令f[i] = dp[i][j]; f[i] = f[i-1] * A; A为转移矩阵

则 S = f[1] + f[2] + f[3] + ... + f[n] = f[1] (A^0 + A^1 +　A^2 + ... + A^n-1);

而f[1] = A * f^0 = A * E;

代入得 S = A (A^0 + A^1 +　A^2 + ... + A^n-1);

代码 
 1 #include <iostream>
 2 #include <cstring>
 3  using namespace std;
 4 
 5 typedef __int64 lld;
 6  const int maxn = 62;
 7 const int mod = 1234567891;
 8 struct Matrix
 9 {
10     int matrix[maxn][maxn];
11 }mat, unit, mat2, unit2;
12 
13 int n, k;
14 
15 void init()
16 {
17     int i, j;
18     memset(mat.matrix, 0, sizeof(mat.matrix));
19     memset(unit.matrix, 0, sizeof(unit.matrix));
20     memset(mat2.matrix, 0, sizeof(mat2.matrix));
21     memset(unit2.matrix, 0, sizeof(unit2.matrix));
22     for(i=1; i<=k; ++i){ //构造A矩阵
23         mat.matrix[i][i] = i;
24         mat.matrix[i-1][i] = k + 1 - i;
25     }
26     mat2 = mat;
27     for(i=0; i<=k; ++i)
28         mat2.matrix[i][i+k+1] = 1;
29     for(i=k+1; i<=2*k+1; ++i)
30         mat2.matrix[i][i] = 1;
31     for(i=0; i<=2*k+1; ++i)
32         unit2.matrix[i][i] = unit.matrix[i][i] = 1;
33 }
34 
35 
36 Matrix Mul(Matrix a, Matrix b)
37 {
38     Matrix c;
39     int i, j, g;
40     for(i=0; i<=2*k+1; ++i){
41         for(j=0; j<=2*k+1; ++j){
42             c.matrix[i][j] = 0;
43             for(g=0; g<=2*k+1; ++g){
44                 c.matrix[i][j] = ( (lld)c.matrix[i][j] + (lld)a.matrix[i][g] * b.matrix[g][j] ) % mod;
45             }
46         }
47     }
48     return c;
49 }
50 
51 
52 
53 int solve()
54 {
55     Matrix p, q;
56     int i, j;
57     p = unit2;
58     q = mat2;
59     while(n){
60         if(n & 1){
61             p = Mul(p, q);
62         }
63         q = Mul(q, q);
64         n >>= 1;
65     }
66     mat2 = p;
67     for(i=0; i<=k; ++i){
68         for(j=0; j<=k; ++j){
69             mat2.matrix[i][j] = mat2.matrix[i][j+k+1];
70             if(i == j) mat2.matrix[i][j] =((lld)mat2.matrix[i][i] + 1) % mod;
71         }
72     }
73     return (int)( (lld)k * mat2.matrix[1][k] % mod);
74 }
75 
76 
77 int main()
78 {
79     int T;
80     scanf("%d", &T);
81     while(T--){
82         scanf("%d %d", &n, &k);
83         init();
84         printf("%d\n", solve());
85     }
86     return 0;
87 }

 6.hdoj 3117 Fibonacci Numbers

首部和尾部分开求

求尾部时矩阵运算，然后取模。

求首部可以用公式，并且由取数字首部的特征可以想到用log运算。

令f[n] = a * 10^b ; (a < 1)

log(f[n]) = log(a) + b;  log(a) < 1;

a = 10 ^ (log(a)); 然后将小数a不断乘10直到四位数

代码 
 1 #include <iostream>
 2 #include <cstring>
 3 #include <cmath>
 4  using namespace std;
 5 
 6 const int maxn = 2;
 7 struct Matrix
 8 {
 9     int matrix[maxn][maxn];
10 }mat, unit;
11 
12 int n;
13 int num[40];
14 
15 
16 Matrix Mul(Matrix a, Matrix b)
17 {
18     Matrix c;
19     int i, j, g;
20     for(i=0; i<2; ++i){
21         for(j=0; j<2; ++j){
22             c.matrix[i][j] = 0;
23             for(g=0; g<2; ++g){
24                 c.matrix[i][j] = (c.matrix[i][j] + a.matrix[i][g] * b.matrix[g][j] ) % 10000;
25             }
26         }
27     }
28     return c;
29 }
30 
31 
32 
33 int getrear()
34 {
35     if(n == 40) return 4155;
36     int k = n - 40;
37     Matrix p, q;
38     p = unit;
39     q = mat;
40     while(k){
41         if(k & 1) p = Mul(p, q);
42         k >>= 1;
43         q = Mul(q, q);
44     }
45     return (p.matrix[1][0] * 5986 + p.matrix[1][1] * 4155) % 10000;    
46 }
47 
48 
49 
50 int gethead()
51 {
52     double t;
53     t = log10(1.0 / sqrt(5.0)) + (double)n * log10( (1+sqrt(5.0))/2.0 );
54     t = t - (int)t;
55     t = pow(10.0, t);
56     while(t < 1000) t *= 10;
57     return (int)t;
58 }
59 
60 
61 void init()
62 {
63     int i, j;
64     num[0] = 0;
65     num[1] = 1;
66     for(i=2; i<=39; ++i) num[i] = num[i-1] + num[i-2];
67     mat.matrix[0][0] = 0;
68     mat.matrix[0][1] = mat.matrix[1][0] = mat.matrix[1][1] = 1;
69     memset(unit.matrix, 0, sizeof(unit.matrix));
70     unit.matrix[0][0] = unit.matrix[1][1] = 1;
71 }
72 
73 
74 
75 int main()
76 {
77     init();
78     while(scanf("%d", &n) != EOF){
79         if(n < 40) printf("%d\n", num[n]);
80         else {
81             int head, rear;
82             head = gethead();
83             rear = getrear();
84             printf("%d...%0.4d\n", head, rear);
85         }
86     }
87     return 0;
88 }

 7.hdoj 2842 Chinese Rings

f[1] = 1; f[2] = 2;

f[n] = 2 * f[n-2] + 1 + f[n-1]; 取第n个的时候先取前n-2个，然后去第n个，然后放上前n-2个，然后取下前n - 1个；

 

代码
 1 #include <iostream>
 2 #include <cstring>
 3 #include <cmath>
 4 using namespace std;
 5 
 6 typedef __int64 lld;
 7 const int maxn = 3;
 8 const int mod = 200907;
 9 struct Matrix
10 {
11     int matrix[maxn][maxn];
12 }mat, unit;
13 
14 int n;
15 
16 
17 
18 Matrix Mul(Matrix a, Matrix b)
19 {
20     Matrix c;
21     int i, j, g;
22     for(i=0; i<3; ++i){
23         for(j=0; j<3; ++j){
24             c.matrix[i][j] = 0;
25             for(g=0; g<3; ++g){
26                 c.matrix[i][j] = ((lld)c.matrix[i][j] + (lld)a.matrix[i][g] * b.matrix[g][j] ) % mod;
27             }
28         }
29     }
30     return c;
31 }
32 
33 
34 
35 int solve()
36 {
37     if(n <= 2) return n;
38     Matrix p, q;
39     p = unit;
40     q = mat;
41     n -= 2;
42     while(n){
43         if(n & 1) p = Mul(p, q);
44         n >>= 1;
45         q = Mul(q, q);
46     }
47     return (p.matrix[2][0] + p.matrix[2][1] + p.matrix[2][2] * 2) % mod;
48 }
49 
50 
51 void init()
52 {
53     int i, j;
54     memset(mat.matrix, 0, sizeof(mat.matrix));
55     memset(unit.matrix, 0, sizeof(unit.matrix));
56     mat.matrix[0][0] = mat.matrix[1][2] = mat.matrix[2][0] = mat.matrix[2][2] = 1;
57     mat.matrix[2][1] = 2;
58     unit.matrix[0][0] = unit.matrix[1][1] = unit.matrix[2][2] = 1;
59 }
60 
61 
62 
63 int main()
64 {
65     init();
66     while(scanf("%d", &n) == 1 && n){
67         printf("%d\n", solve());
68     }
69     return 0;
70 }
71 

 8.hdojLucky Coins Sequence

f[1]  = 0; f[2] = 0;

f[3] = 2;

f[4] = 6;

f[5] = 16;

f[6] = 38;

f[7] = 86;

f[8] = 188;

f[4] - f[3] * 2 = 2;

f[5] - f[4] * 2 = 4;

f[6] - f[5] * 2 = 6;

f[7] - f[6] * 2 = 10;

等号右边为一等差数列，从而可以构造矩阵来解决问题

代码
#include <iostream>
#include <cstring>
using namespace std;

const int maxn = 3;
const int mod = 10007;
struct Matrix
{
    int matrix[maxn][maxn];
}mat, unit;
int n;

Matrix Mul(Matrix a, Matrix b)
{
    int i, j, k;
    Matrix c;
    for(i=0; i<maxn; ++i){
        for(j=0; j<maxn; ++j){
            c.matrix[i][j] = 0;
            for(k=0; k<maxn; ++k){
                c.matrix[i][j] += (a.matrix[i][k] * b.matrix[k][j]) % mod;
                c.matrix[i][j] %= mod;
            }
        }
    }
    return c;
}


int solve()
{
    if(n == 1 || n == 2) return 0;
    if(n == 3) return 2;
    if(n == 4) return 6;
    int k = n - 4;
    Matrix p, q;
    p = unit;
    q = mat;
    while(k){
        if(k & 1) p = Mul(p, q);
        q = Mul(q, q);
        k >>= 1;
    }
    return (p.matrix[2][0] * 2 + p.matrix[2][1] * 4 + p.matrix[2][2] * 6) % mod;
}


void init()
{
    int i;
    memset(unit.matrix, 0, sizeof(unit.matrix));
    memset(mat.matrix, 0, sizeof(mat.matrix));
    for(i=0; i<3; ++i){
        unit.matrix[i][i] = 1;
    }
    mat.matrix[0][1] = mat.matrix[1][0] = mat.matrix[1][1] = mat.matrix[2][1] = 1;
    mat.matrix[2][2] = 2;
}


int main()
{
    init();
    while(scanf("%d", &n)!=EOF){
        printf("%d\n", solve());
    }
    return 0;
}

 

 

 

posted on 2010-11-12 12:02  CrazyAC  阅读(293)  评论(0)  编辑  收藏  举报