Differences From Artifact [b6522f93bcade295]:
- File
lib/numeric/modArith.cpp
- 2012-09-02 11:59:33 - part of checkin [0ea7f42c5c] on branch trunk - 554 (user: kinaba) [annotate]
To Artifact [065fcd1fef379fbc]:
- File
lib/numeric/modArith.cpp
- 2012-09-08 04:27:24 - part of checkin [7d2fe88cb2] on branch trunk - Cleaned up matrix multication library. (user: kinaba) [annotate]
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 */ <