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 -*/