23dfcca431 2011-02-23 kinaba: 23dfcca431 2011-02-23 kinaba: //------------------------------------------------------------- 23dfcca431 2011-02-23 kinaba: // Modulo Arithmetics 23dfcca431 2011-02-23 kinaba: // 23dfcca431 2011-02-23 kinaba: // Verified by 23dfcca431 2011-02-23 kinaba: // - TCO10 R3 LV3 16685c1758 2012-06-07 kinaba: // - SRM 545 LV2 23dfcca431 2011-02-23 kinaba: //------------------------------------------------------------- 23dfcca431 2011-02-23 kinaba: 16685c1758 2012-06-07 kinaba: static const int MODVAL = 1000000007; 23dfcca431 2011-02-23 kinaba: struct mint 23dfcca431 2011-02-23 kinaba: { 23dfcca431 2011-02-23 kinaba: int val; 23dfcca431 2011-02-23 kinaba: mint():val(0){} 16685c1758 2012-06-07 kinaba: mint(int x):val(x%MODVAL) {} 16685c1758 2012-06-07 kinaba: mint(size_t x):val(x%MODVAL) {} 16685c1758 2012-06-07 kinaba: mint(LL x):val(x%MODVAL) {} 23dfcca431 2011-02-23 kinaba: }; 068d2337c1 2011-09-18 kinaba: mint& operator+=(mint& x, mint y) { return x = x.val+y.val; } 068d2337c1 2011-09-18 kinaba: mint& operator-=(mint& x, mint y) { return x = x.val-y.val+MODVAL; } 068d2337c1 2011-09-18 kinaba: mint& operator*=(mint& x, mint y) { return x = LL(x.val)*y.val; } 068d2337c1 2011-09-18 kinaba: mint POW(mint x, LL e) { mint v=1; for(;e;x*=x,e>>=1) if(e&1) v*=x; return v; } 068d2337c1 2011-09-18 kinaba: mint& operator/=(mint& x, mint y) { return x *= POW(y, MODVAL-2); } 068d2337c1 2011-09-18 kinaba: mint operator+(mint x, mint y) { return x+=y; } 068d2337c1 2011-09-18 kinaba: mint operator-(mint x, mint y) { return x-=y; } 068d2337c1 2011-09-18 kinaba: mint operator*(mint x, mint y) { return x*=y; } 068d2337c1 2011-09-18 kinaba: mint operator/(mint x, mint y) { return x/=y; } 23dfcca431 2011-02-23 kinaba: vector<mint> FAC_(1,1); 068d2337c1 2011-09-18 kinaba: mint FAC(LL n) { while( FAC_.size()<=n ) FAC_.push_back( FAC_.back()*FAC_.size() ); return FAC_[n]; } 068d2337c1 2011-09-18 kinaba: mint C(LL n, LL k) { return k<0 || n<k ? 0 : FAC(n) / (FAC(k) * FAC(n-k)); } 4318dd2827 2012-06-07 kinaba: vector< vector<mint> > CP_; // Pascal's triangle: if O(1) nCk is needed. 16685c1758 2012-06-07 kinaba: mint C(LL n, LL k) { 4318dd2827 2012-06-07 kinaba: while( CP_.size() <= n ) { 4318dd2827 2012-06-07 kinaba: int nn = CP_.size(); 4318dd2827 2012-06-07 kinaba: CP_.push_back(vector<mint>(nn+1,1)); 4318dd2827 2012-06-07 kinaba: for(int kk=1; kk<nn; ++kk) 4318dd2827 2012-06-07 kinaba: CP_[nn][kk] = CP_[nn-1][kk-1] + CP_[nn-1][kk]; 16685c1758 2012-06-07 kinaba: } 16685c1758 2012-06-07 kinaba: return k<0 || n<k ? 0 : CP_[n][k]; 16685c1758 2012-06-07 kinaba: } 16685c1758 2012-06-07 kinaba: c81888b809 2012-08-06 kinaba: template<typename T> c81888b809 2012-08-06 kinaba: vector<T> MATMUL(const vector< vector<T> >& a, const vector<T>& v) c81888b809 2012-08-06 kinaba: { c81888b809 2012-08-06 kinaba: int N = a.size(); c81888b809 2012-08-06 kinaba: vector<T> u(N); c81888b809 2012-08-06 kinaba: for(int i=0; i<N; ++i) c81888b809 2012-08-06 kinaba: for(int j=0; j<N; ++j) c81888b809 2012-08-06 kinaba: u[i] += a[i][j]*v[j]; c81888b809 2012-08-06 kinaba: return u; c81888b809 2012-08-06 kinaba: } c81888b809 2012-08-06 kinaba: c81888b809 2012-08-06 kinaba: template<typename T> c81888b809 2012-08-06 kinaba: vector< vector<T> > MATMUL(const vector< vector<T> >& a, const vector< vector<T> >& b) c81888b809 2012-08-06 kinaba: { c81888b809 2012-08-06 kinaba: int N = a.size(); c81888b809 2012-08-06 kinaba: vector< vector<T> > c(N, vector<T>(N)); c81888b809 2012-08-06 kinaba: for(int i=0; i<N; ++i) c81888b809 2012-08-06 kinaba: for(int j=0; j<N; ++j) c81888b809 2012-08-06 kinaba: for(int k=0; k<N; ++k) c81888b809 2012-08-06 kinaba: c[i][j] += a[i][k]*b[k][j]; c81888b809 2012-08-06 kinaba: return c; c81888b809 2012-08-06 kinaba: } c81888b809 2012-08-06 kinaba: c81888b809 2012-08-06 kinaba: template<typename T> c81888b809 2012-08-06 kinaba: vector< vector<T> > MATPOW(vector< vector<T> > a, LL e) c81888b809 2012-08-06 kinaba: { c81888b809 2012-08-06 kinaba: int N = a.size(); c81888b809 2012-08-06 kinaba: vector< vector<T> > c(N, vector<T>(N)); c81888b809 2012-08-06 kinaba: for(int i=0; i<N; ++i) c[i][i] = 1; c81888b809 2012-08-06 kinaba: for(; e; e>>=1) { c81888b809 2012-08-06 kinaba: if(e&1) c81888b809 2012-08-06 kinaba: c = MATMUL(c, a); c81888b809 2012-08-06 kinaba: a = MATMUL(a, a); c81888b809 2012-08-06 kinaba: } c81888b809 2012-08-06 kinaba: return c; c81888b809 2012-08-06 kinaba: } c81888b809 2012-08-06 kinaba: 23dfcca431 2011-02-23 kinaba: /* 23dfcca431 2011-02-23 kinaba: // MODVAL must be a prime!! 23dfcca431 2011-02-23 kinaba: LL GSS(LL k, LL b, LL e) // k^b + k^b+1 + ... + k^e 23dfcca431 2011-02-23 kinaba: { 23dfcca431 2011-02-23 kinaba: if( b > e ) return 0; 23dfcca431 2011-02-23 kinaba: if( k <= 1 ) return k*(e-b+1); 23dfcca431 2011-02-23 kinaba: return DIV(SUB(POW(k, e+1), POW(k,b)), k-1); 23dfcca431 2011-02-23 kinaba: } 23dfcca431 2011-02-23 kinaba: 23dfcca431 2011-02-23 kinaba: // works for non-prime MODVAL 23dfcca431 2011-02-23 kinaba: LL GEO(LL x_, LL e) // x^0 + x^1 + ... + x^e-1 23dfcca431 2011-02-23 kinaba: { 23dfcca431 2011-02-23 kinaba: vector< vector<LL> > v(2, vector<LL>(2)); 23dfcca431 2011-02-23 kinaba: vector< vector<LL> > x(2, vector<LL>(2)); 23dfcca431 2011-02-23 kinaba: v[0][0] = v[1][1] = 1; 23dfcca431 2011-02-23 kinaba: x[0][0] = x_; x[0][1] = 0; 23dfcca431 2011-02-23 kinaba: x[1][0] = 1 ; x[1][1] = 1; 23dfcca431 2011-02-23 kinaba: for(;e;x=MATMUL(x,x),e>>=1) 23dfcca431 2011-02-23 kinaba: if(e&1) 23dfcca431 2011-02-23 kinaba: v = MATMUL(v, x); 23dfcca431 2011-02-23 kinaba: return v[1][0]; 23dfcca431 2011-02-23 kinaba: } 23dfcca431 2011-02-23 kinaba: 23dfcca431 2011-02-23 kinaba: // works for non-prime MODVAL 23dfcca431 2011-02-23 kinaba: 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) 23dfcca431 2011-02-23 kinaba: { 23dfcca431 2011-02-23 kinaba: vector< vector<LL> > v(3, vector<LL>(3)); 23dfcca431 2011-02-23 kinaba: vector< vector<LL> > x(3, vector<LL>(3)); 23dfcca431 2011-02-23 kinaba: v[0][0] = v[1][1] = v[2][2] = 1; 23dfcca431 2011-02-23 kinaba: x[0][0] = x_; x[0][1] = 0; x[0][2] = 0; 23dfcca431 2011-02-23 kinaba: x[1][0] = 1 ; x[1][1] = 1; x[1][2] = 0; 23dfcca431 2011-02-23 kinaba: x[2][0] = 0 ; x[2][1] = 1; x[2][2] = 1; 23dfcca431 2011-02-23 kinaba: e++; 23dfcca431 2011-02-23 kinaba: for(;e;x=MATMUL(x,x),e>>=1) 23dfcca431 2011-02-23 kinaba: if(e&1) 23dfcca431 2011-02-23 kinaba: v = MATMUL(v, x); 23dfcca431 2011-02-23 kinaba: return v[2][0]; 23dfcca431 2011-02-23 kinaba: } 23dfcca431 2011-02-23 kinaba: */