Differences From Artifact [2c92483914a46a56]:
- File
_lib/numeric/matrix.cpp
- 2011-02-23 09:21:16 - part of checkin [4fd800b3a8] on branch trunk - Copied from private svn repository. (user: kinaba) [annotate]
- File
lib/numeric/matrix.cpp
- 2011-02-23 11:18:09 - part of checkin [23dfcca431] on branch trunk - renamed _lib to lib (user: kinaba) [annotate]
To Artifact [ba1b2d78835ef59e]:
- File
lib/numeric/matrix.cpp
- 2012-09-08 04:27:24 - part of checkin [7d2fe88cb2] on branch trunk - Cleaned up matrix multication library. (user: kinaba) [annotate]
1 -
2 1 //-------------------------------------------------------------
3 2 // Matrix Operations
4 3 //
5 4 // Verified by
6 5 // - SRM342 Div1 LV3
7 6 // - SRM341 Div1 LV3
8 7 // - SRM338 Div1 LV2
9 8 //-------------------------------------------------------------
10 9
11 -vector<LL> vMul( const vector< vector<LL> >& A, const vector<LL>& B )
10 +template<typename T>
11 +vector<T> MATMUL(const vector< vector<T> >& a, const vector<T>& v)
12 +{
13 + int N = a.size();
14 + vector<T> u(N);
15 + for(int i=0; i<N; ++i)
16 + for(int j=0; j<N; ++j)
17 + u[i] += a[i][j]*v[j];
18 + return u;
19 +}
20 +
21 +template<typename T>
22 +vector< vector<T> > MATMUL(const vector< vector<T> >& a, const vector< vector<T> >& b)
23 +{
24 + int N = a.size();
25 + vector< vector<T> > c(N, vector<T>(N));
26 + for(int i=0; i<N; ++i)
27 + for(int j=0; j<N; ++j)
28 + for(int k=0; k<N; ++k)
29 + c[i][j] += a[i][k]*b[k][j];
30 + return c;
31 +}
32 +
33 +template<typename T>
34 +vector<T> MATPOWMUL(vector< vector<T> > a, LL e, vector<T> v)
12 35 {
13 - const int n = A.size();
14 -
15 - vector<LL> C(n);
16 - for(int i=0; i<n; ++i)
17 - for(int j=0; j<n; ++j)
18 - C[i] += A[i][j] * B[j];
19 - return C;
36 + for(; e; e>>=1, a=MATMUL(a,a))
37 + if(e&1)
38 + v = MATMUL(a, v);
39 + return v;
20 40 }
21 41
22 -vector< vector<LL> > mMul( const vector< vector<LL> >& A, const vector< vector<LL> >& B )
42 +template<typename T>
43 +vector< vector<T> > MATPOW(vector< vector<T> > a, LL e)
23 44 {
24 - const int n = A.size();
25 -
26 - vector< vector<LL> > C(n, vector<LL>(n));
27 - for(int i=0; i<n; ++i)
28 - for(int j=0; j<n; ++j) {
29 - LL Cij = 0;
30 - for(int k=0; k<n; ++k)
31 - Cij += A[i][k] * B[k][j];
32 - C[i][j] = Cij;
33 - }
34 - return C;
45 + int N = a.size();
46 + vector< vector<T> > c(N, vector<T>(N));
47 + for(int i=0; i<N; ++i) c[i][i] = 1;
48 + for(; e; e>>=1) {
49 + if(e&1)
50 + c = MATMUL(c, a);
51 + a = MATMUL(a, a);
52 + }
53 + return c;
35 54 }
36 55
37 -vector< vector<LL> > mPow( vector< vector<LL> > M, LL t ) // t>0
38 -{
39 - vector< vector<LL> > R;
40 - for(; t; t>>=1, M=mMul(M,M))
41 - if( t&1 )
42 - R = (R.empty() ? M : mMul(R,M));
43 - return R;
44 -}
56 +/****************************************************************************
57 + NOTES ON THE USAGE
58 + ****************************************************************************
59 +
60 +A[to][from] = #transition
61 +--------------------------------------------
62 +
63 + frm
64 + (. . .) (q0)
65 +to (. . .) (q1)
66 + (. . .) (q2)
67 +
68 +
69 +
70 + qi
71 +--------------------------------------------
72 +
73 + frm
74 + (. . . 0) (.)
75 +to (. . . 0) (.)
76 + (. . . 0) (.)
77 + (0 1 0 1) (0)
78 + ^
79 + i
80 +Then,
81 + q[i]@0 = (A^1 v)[$-1]
82 + q[i]@0 + q[i]@1 = (A^2 v)[$-1]
83 + q[i]@0 + ... + q[i]@k = (A^(k+1) v)[$-1]
84 +
85 +
86 +
87 + q_all
88 +--------------------------------------------
89 +
90 + frm
91 + (. . . 0) (.)
92 +to (. . . 0) (.)
93 + (. . . 0) (.)
94 + (1 1 1 1) (0)
95 +
96 +Then,
97 + q @0 = (A^1 v)[$-1]
98 + q @0 + ... + q @k = (A^(k+1) v)[$-1]
99 +
100 +
101 +
102 +
103 +
104 +Application: x^0 + x^1 + ... + x^e-1
105 +--------------------------------------------
106 +
107 + (x 0)^e (1) = (x^e )
108 + (1 1) (0) (x^0 + x^1 + ... + x^e-1)
109 +
110 +
45 111
46 -vector< vector<LL> > mAdd( const vector< vector<LL> >& A, const vector< vector<LL> >& B )
47 -{
48 - const int n = A.size();
112 +Application: e x^0 + (e-1) x^1 + ... + 1 x^e-1
113 +--------------------------------------------
49 114
50 - vector< vector<LL> > C(n, vector<LL>(n));
51 - for(int i=0; i<n; ++i)
52 - for(int j=0; j<n; ++j)
53 - C[i][j] = A[i][j] + B[i][j];
54 - return C;
55 -}
115 + (x 0 0)^e (x) = (x^e+1 )
116 + (1 1 0) (1) (x^0 + x^1 + ... + x^e )
117 + (0 1 1) (0) ( of = e x^0 + (e-1) x^1 + ... + 1 x^e-1)
118 +*/