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 // Matrix Operations 2 // Matrix Operations
4 // 3 //
5 // Verified by 4 // Verified by
6 // - SRM342 Div1 LV3 5 // - SRM342 Div1 LV3
7 // - SRM341 Div1 LV3 6 // - SRM341 Div1 LV3
8 // - SRM338 Div1 LV2 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>
> 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(); | 36 for(; e; e>>=1, a=MATMUL(a,a))
14 <
> 37 if(e&1)
15 vector<LL> C(n); | 38 v = MATMUL(a, v);
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; | 39 return v;
20 } 40 }
21 41
22 vector< vector<LL> > mMul( const vector< vector<LL> >& A, const vector< vector<L | 42 template<typename T>
> 43 vector< vector<T> > MATPOW(vector< vector<T> > a, LL e)
23 { 44 {
24 const int n = A.size(); | 45 int N = a.size();
25 <
26 vector< vector<LL> > C(n, vector<LL>(n)); | 46 vector< vector<T> > c(N, vector<T>(N));
27 for(int i=0; i<n; ++i) | 47 for(int i=0; i<N; ++i) c[i][i] = 1;
> 48 for(; e; e>>=1) {
28 for(int j=0; j<n; ++j) { | 49 if(e&1)
29 LL Cij = 0; | 50 c = MATMUL(c, a);
30 for(int k=0; k<n; ++k) | 51 a = MATMUL(a, a);
31 Cij += A[i][k] * B[k][j]; <
32 C[i][j] = Cij; <
33 } | 52 }
34 return C; | 53 return c;
35 } 54 }
36 55
37 vector< vector<LL> > mPow( vector< vector<LL> > M, LL t ) // t>0 | 56 /****************************************************************************
> 57 NOTES ON THE USAGE
> 58 ****************************************************************************
38 { | 59
39 vector< vector<LL> > R; <
40 for(; t; t>>=1, M=mMul(M,M)) <
41 if( t&1 ) | 60 A[to][from] = #transition
42 R = (R.empty() ? M : mMul(R,M)); <
43 return R; <
> 61 --------------------------------------------
44 } | 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<L | 112 Application: e x^0 + (e-1) x^1 + ... + 1 x^e-1
47 { <
48 const int n = A.size(); <
> 113 --------------------------------------------
49 114
50 vector< vector<LL> > C(n, vector<LL>(n)); | 115 (x 0 0)^e (x) = (x^e+1 )
51 for(int i=0; i<n; ++i) | 116 (1 1 0) (1) (x^0 + x^1 + ... + x^e )
52 for(int j=0; j<n; ++j) | 117 (0 1 1) (0) ( of = e x^0 + (e-1) x^1 + ... + 1 x^e-1)
53 C[i][j] = A[i][j] + B[i][j]; <
54 return C; <
55 } <
> 118 */