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
23dfcca431 2011-02-23        kinaba: //-------------------------------------------------------------
23dfcca431 2011-02-23        kinaba: 
23dfcca431 2011-02-23        kinaba: static const int MODVAL = 1000000007; // op/ depends on primarity of the MODVAL
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){}
23dfcca431 2011-02-23        kinaba: 	mint(int x):val(x%MODVAL) {}
23dfcca431 2011-02-23        kinaba: 	mint(LL  x):val(x%MODVAL) {}
23dfcca431 2011-02-23        kinaba: };
23dfcca431 2011-02-23        kinaba: 
23dfcca431 2011-02-23        kinaba: mint operator+(mint x, mint y) { return x.val+y.val; }
23dfcca431 2011-02-23        kinaba: mint operator-(mint x, mint y) { return x.val-y.val+MODVAL; }
23dfcca431 2011-02-23        kinaba: mint operator*(mint x, mint y) { return LL(x.val)*y.val; }
23dfcca431 2011-02-23        kinaba: mint POW(mint x, int e) {
23dfcca431 2011-02-23        kinaba: 	mint v = 1;
23dfcca431 2011-02-23        kinaba: 	for(;e;x=x*x,e>>=1)
23dfcca431 2011-02-23        kinaba: 		if(e&1)
23dfcca431 2011-02-23        kinaba: 			v=v*x;
23dfcca431 2011-02-23        kinaba: 	return v;
23dfcca431 2011-02-23        kinaba: }
23dfcca431 2011-02-23        kinaba: mint operator/(mint x, mint y) { return x * POW(y, MODVAL-2); }
23dfcca431 2011-02-23        kinaba: 
23dfcca431 2011-02-23        kinaba: vector<mint> FAC_(1,1);
23dfcca431 2011-02-23        kinaba: void FAC_INIT(int n) { for(int i=1; i<=n; ++i) FAC_.push_back( FAC_.back()*i ); }
23dfcca431 2011-02-23        kinaba: mint FAC(mint x)       { return FAC_[x.val]; }
23dfcca431 2011-02-23        kinaba: mint C(mint n, mint k) { return k.val<0 || n.val<k.val ? 0 : FAC(n) / (FAC(k) * FAC(n-k)); }
23dfcca431 2011-02-23        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: LL Cpascal(LL n, LL k)
23dfcca431 2011-02-23        kinaba: {
23dfcca431 2011-02-23        kinaba: 	vector< vector<LL> > c(n+1, vector<LL>(k+1));
23dfcca431 2011-02-23        kinaba: 	for(LL nn=1; nn<=n; ++nn)
23dfcca431 2011-02-23        kinaba: 		for(LL kk=0; kk<=min(nn,k); ++kk)
23dfcca431 2011-02-23        kinaba: 			c[nn][kk] = kk==0 || kk==nn ? 1 : ADD(c[nn-1][kk-1], c[nn-1][kk]);
23dfcca431 2011-02-23        kinaba: 	return c[n][k];
23dfcca431 2011-02-23        kinaba: }
23dfcca431 2011-02-23        kinaba: 
23dfcca431 2011-02-23        kinaba: vector< vector<LL> > MATMUL(vector< vector<LL> >& a, vector< vector<LL> >& b)
23dfcca431 2011-02-23        kinaba: {
23dfcca431 2011-02-23        kinaba:    int N = a.size();
23dfcca431 2011-02-23        kinaba:    vector< vector<LL> > c(N, vector<LL>(N));
23dfcca431 2011-02-23        kinaba:    for(int i=0; i<N; ++i)
23dfcca431 2011-02-23        kinaba:    for(int j=0; j<N; ++j)
23dfcca431 2011-02-23        kinaba:    for(int k=0; k<N; ++k)
23dfcca431 2011-02-23        kinaba:       c[i][j] = ADD(c[i][j], MUL(a[i][k],b[k][j]));
23dfcca431 2011-02-23        kinaba:    return c;
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: */