Overview
| SHA1 Hash: | 7d2fe88cb2c970ed27cca0f03708a6da748725e6 |
|---|---|
| Date: | 2012-09-08 04:27:24 |
| User: | kinaba |
| Comment: | Cleaned up matrix multication library. |
| Timelines: | family | ancestors | descendants | both | trunk |
| Downloads: | Tarball | ZIP archive |
| Other Links: | files | file ages | manifest |
Tags And Properties
- branch=trunk inherited from [9165bd3629]
- sym-trunk inherited from [9165bd3629]
Changes
Modified lib/numeric/matrix.cpp from [2c92483914a46a56] to [ba1b2d78835ef59e].
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 +*/
Modified lib/numeric/modArith.cpp from [b6522f93bcade295] to [065fcd1fef379fbc].
1 1 2 2 //------------------------------------------------------------- 3 3 // Modulo Arithmetics 4 4 // 5 5 // Verified by 6 6 // - TCO10 R3 LV3 7 -// - SRM 545 LV2 7 +// - SRM 545 Div1 LV2 8 +// - SRM 554 Div1 LV3 8 9 //------------------------------------------------------------- 9 10 10 11 static const unsigned MODVAL = 1000000007; 11 12 struct mint 12 13 { 13 14 unsigned val; 14 15 mint():val(0){} ................................................................................ 23 24 mint operator-(mint x, mint y) { return x-=y; } 24 25 mint operator*(mint x, mint y) { return x*=y; } 25 26 26 27 mint POW(mint x, LL e) { mint v=1; for(;e;x*=x,e>>=1) if(e&1) v*=x; return v; } 27 28 mint& operator/=(mint& x, mint y) { return x *= POW(y, MODVAL-2); } 28 29 mint operator/(mint x, mint y) { return x/=y; } 29 30 30 -// nCk :: O(log MODVAL) time, O(n) space. 31 31 vector<mint> FAC_(1,1); 32 32 mint FAC(LL n) { while( FAC_.size()<=n ) FAC_.push_back( FAC_.back()*FAC_.size() ); return FAC_[n]; } 33 + 34 +// nCk :: O(log MODVAL) time, O(n) space. 33 35 mint C(LL n, LL k) { return k<0 || n<k ? 0 : FAC(n) / (FAC(k) * FAC(n-k)); } 34 36 35 37 // nCk :: O(1) time, O(n^2) space. 36 38 vector< vector<mint> > CP_; 37 39 mint C(LL n, LL k) { 38 40 while( CP_.size() <= n ) { 39 41 int nn = CP_.size(); ................................................................................ 40 42 CP_.push_back(vector<mint>(nn+1,1)); 41 43 for(int kk=1; kk<nn; ++kk) 42 44 CP_[nn][kk] = CP_[nn-1][kk-1] + CP_[nn-1][kk]; 43 45 } 44 46 return k<0 || n<k ? 0 : CP_[n][k]; 45 47 } 46 48 47 -template<typename T> 48 -vector<T> MATMUL(const vector< vector<T> >& a, const vector<T>& v) 49 +// O(log MODVAL), MODVAL must be prime: k^b + k^b+1 + ... + k^e 50 +mint GSS(mint k, LL b, LL e) 49 51 { 50 - int N = a.size(); 51 - vector<T> u(N); 52 - for(int i=0; i<N; ++i) 53 - for(int j=0; j<N; ++j) 54 - u[i] += a[i][j]*v[j]; 55 - return u; 52 + if( b > e ) return 0; 53 + if( k.val <= 1 ) return k*(e-b+1); 54 + return (POW(k, e+1) - POW(k,b)) / (k-1); 56 55 } 57 - 58 -template<typename T> 59 -vector< vector<T> > MATMUL(const vector< vector<T> >& a, const vector< vector<T> >& b) 60 -{ 61 - int N = a.size(); 62 - vector< vector<T> > c(N, vector<T>(N)); 63 - for(int i=0; i<N; ++i) 64 - for(int j=0; j<N; ++j) 65 - for(int k=0; k<N; ++k) 66 - c[i][j] += a[i][k]*b[k][j]; 67 - return c; 68 -} 69 - 70 -template<typename T> 71 -vector< vector<T> > MATPOW(vector< vector<T> > a, LL e) 72 -{ 73 - int N = a.size(); 74 - vector< vector<T> > c(N, vector<T>(N)); 75 - for(int i=0; i<N; ++i) c[i][i] = 1; 76 - for(; e; e>>=1) { 77 - if(e&1) 78 - c = MATMUL(c, a); 79 - a = MATMUL(a, a); 80 - } 81 - return c; 82 -} 83 - 84 -/* 85 -// MODVAL must be a prime!! 86 -LL GSS(LL k, LL b, LL e) // k^b + k^b+1 + ... + k^e 87 -{ 88 - if( b > e ) return 0; 89 - if( k <= 1 ) return k*(e-b+1); 90 - return DIV(SUB(POW(k, e+1), POW(k,b)), k-1); 91 -} 92 - 93 -// works for non-prime MODVAL 94 -LL GEO(LL x_, LL e) // x^0 + x^1 + ... + x^e-1 95 -{ 96 - vector< vector<LL> > v(2, vector<LL>(2)); 97 - vector< vector<LL> > x(2, vector<LL>(2)); 98 - v[0][0] = v[1][1] = 1; 99 - x[0][0] = x_; x[0][1] = 0; 100 - x[1][0] = 1 ; x[1][1] = 1; 101 - for(;e;x=MATMUL(x,x),e>>=1) 102 - if(e&1) 103 - v = MATMUL(v, x); 104 - return v[1][0]; 105 -} 106 - 107 -// works for non-prime MODVAL 108 -LL HYP(LL x_, LL e) // e x^0 + (e-1) x^1 + ... + 1 x^e-1 = GEO(x,1)+GEO(x,2)+...+GEO(x,e) 109 -{ 110 - vector< vector<LL> > v(3, vector<LL>(3)); 111 - vector< vector<LL> > x(3, vector<LL>(3)); 112 - v[0][0] = v[1][1] = v[2][2] = 1; 113 - x[0][0] = x_; x[0][1] = 0; x[0][2] = 0; 114 - x[1][0] = 1 ; x[1][1] = 1; x[1][2] = 0; 115 - x[2][0] = 0 ; x[2][1] = 1; x[2][2] = 1; 116 - e++; 117 - for(;e;x=MATMUL(x,x),e>>=1) 118 - if(e&1) 119 - v = MATMUL(v, x); 120 - return v[2][0]; 121 -} 122 -*/