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