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 // 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 */
Modified lib/numeric/modArith.cpp from [b6522f93bcade295] to [065fcd1fef379fbc].
1 1 2 //------------------------------------------------------------- 2 //------------------------------------------------------------- 3 // Modulo Arithmetics 3 // Modulo Arithmetics 4 // 4 // 5 // Verified by 5 // Verified by 6 // - TCO10 R3 LV3 6 // - TCO10 R3 LV3 7 // - SRM 545 LV2 | 7 // - SRM 545 Div1 LV2 > 8 // - SRM 554 Div1 LV3 8 //------------------------------------------------------------- 9 //------------------------------------------------------------- 9 10 10 static const unsigned MODVAL = 1000000007; 11 static const unsigned MODVAL = 1000000007; 11 struct mint 12 struct mint 12 { 13 { 13 unsigned val; 14 unsigned val; 14 mint():val(0){} 15 mint():val(0){} ................................................................................................................................................................................ 23 mint operator-(mint x, mint y) { return x-=y; } 24 mint operator-(mint x, mint y) { return x-=y; } 24 mint operator*(mint x, mint y) { return x*=y; } 25 mint operator*(mint x, mint y) { return x*=y; } 25 26 26 mint POW(mint x, LL e) { mint v=1; for(;e;x*=x,e>>=1) if(e&1) v*=x; return v; } 27 mint POW(mint x, LL e) { mint v=1; for(;e;x*=x,e>>=1) if(e&1) v*=x; return v; } 27 mint& operator/=(mint& x, mint y) { return x *= POW(y, MODVAL-2); } 28 mint& operator/=(mint& x, mint y) { return x *= POW(y, MODVAL-2); } 28 mint operator/(mint x, mint y) { return x/=y; } 29 mint operator/(mint x, mint y) { return x/=y; } 29 30 30 // nCk :: O(log MODVAL) time, O(n) space. < 31 vector<mint> FAC_(1,1); 31 vector<mint> FAC_(1,1); 32 mint FAC(LL n) { while( FAC_.size()<=n ) FAC_.push_back( FAC_.back()*FAC_.size() 32 mint FAC(LL n) { while( FAC_.size()<=n ) FAC_.push_back( FAC_.back()*FAC_.size() > 33 > 34 // nCk :: O(log MODVAL) time, O(n) space. 33 mint C(LL n, LL k) { return k<0 || n<k ? 0 : FAC(n) / (FAC(k) * FAC(n-k)); } 35 mint C(LL n, LL k) { return k<0 || n<k ? 0 : FAC(n) / (FAC(k) * FAC(n-k)); } 34 36 35 // nCk :: O(1) time, O(n^2) space. 37 // nCk :: O(1) time, O(n^2) space. 36 vector< vector<mint> > CP_; 38 vector< vector<mint> > CP_; 37 mint C(LL n, LL k) { 39 mint C(LL n, LL k) { 38 while( CP_.size() <= n ) { 40 while( CP_.size() <= n ) { 39 int nn = CP_.size(); 41 int nn = CP_.size(); ................................................................................................................................................................................ 40 CP_.push_back(vector<mint>(nn+1,1)); 42 CP_.push_back(vector<mint>(nn+1,1)); 41 for(int kk=1; kk<nn; ++kk) 43 for(int kk=1; kk<nn; ++kk) 42 CP_[nn][kk] = CP_[nn-1][kk-1] + CP_[nn-1][kk]; 44 CP_[nn][kk] = CP_[nn-1][kk-1] + CP_[nn-1][kk]; 43 } 45 } 44 return k<0 || n<k ? 0 : CP_[n][k]; 46 return k<0 || n<k ? 0 : CP_[n][k]; 45 } 47 } 46 48 47 template<typename T> | 49 // O(log MODVAL), MODVAL must be prime: k^b + k^b+1 + ... + k^e 48 vector<T> MATMUL(const vector< vector<T> >& a, const vector<T>& v) | 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> < 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)+... < 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 */ <