4fd800b3a8 2011-02-23 kinaba: 4fd800b3a8 2011-02-23 kinaba: //------------------------------------------------------------- 4fd800b3a8 2011-02-23 kinaba: // Modulo Arithmetics 4fd800b3a8 2011-02-23 kinaba: // 4fd800b3a8 2011-02-23 kinaba: // Verified by 4fd800b3a8 2011-02-23 kinaba: // - SRM 397 Div1 LV2 4fd800b3a8 2011-02-23 kinaba: // - SRM 428 Div1 LV2 4fd800b3a8 2011-02-23 kinaba: // - CodeCraft 2010 CNTINT DRAW 4fd800b3a8 2011-02-23 kinaba: //------------------------------------------------------------- 4fd800b3a8 2011-02-23 kinaba: 4fd800b3a8 2011-02-23 kinaba: static const LL MODVAL = 1000000007; // must fit in 32-bits 4fd800b3a8 2011-02-23 kinaba: 4fd800b3a8 2011-02-23 kinaba: LL ADD(LL x, LL y) { return (x+y)%MODVAL; } 4fd800b3a8 2011-02-23 kinaba: LL SUB(LL x, LL y) { return (x-y+MODVAL)%MODVAL; } 4fd800b3a8 2011-02-23 kinaba: LL MUL(LL x, LL y) { return (x*y)%MODVAL; } 4fd800b3a8 2011-02-23 kinaba: LL POW(LL x, LL e) { 4fd800b3a8 2011-02-23 kinaba: LL v = 1; 4fd800b3a8 2011-02-23 kinaba: for(;e;x=MUL(x,x),e>>=1) 4fd800b3a8 2011-02-23 kinaba: if(e&1) 4fd800b3a8 2011-02-23 kinaba: v = MUL(v, x); 4fd800b3a8 2011-02-23 kinaba: return v; 4fd800b3a8 2011-02-23 kinaba: } 4fd800b3a8 2011-02-23 kinaba: 4fd800b3a8 2011-02-23 kinaba: // MODVAL must be a prime!! 4fd800b3a8 2011-02-23 kinaba: LL DIV(LL x, LL y) { return MUL(x, POW(y, MODVAL-2)); } 4fd800b3a8 2011-02-23 kinaba: 4fd800b3a8 2011-02-23 kinaba: // MODVAL must be a prime!! 4fd800b3a8 2011-02-23 kinaba: LL C(LL n, LL k) { 4fd800b3a8 2011-02-23 kinaba: LL v = 1; 4fd800b3a8 2011-02-23 kinaba: for(LL i=1; i<=k; ++i) 4fd800b3a8 2011-02-23 kinaba: v = DIV(MUL(v, n-i+1), i); 4fd800b3a8 2011-02-23 kinaba: return v; 4fd800b3a8 2011-02-23 kinaba: } 4fd800b3a8 2011-02-23 kinaba: 4fd800b3a8 2011-02-23 kinaba: // MODVAL must be a prime!! 4fd800b3a8 2011-02-23 kinaba: LL GSS(LL k, LL b, LL e) // k^b + k^b+1 + ... + k^e 4fd800b3a8 2011-02-23 kinaba: { 4fd800b3a8 2011-02-23 kinaba: if( b > e ) return 0; 4fd800b3a8 2011-02-23 kinaba: if( k <= 1 ) return k*(e-b+1); 4fd800b3a8 2011-02-23 kinaba: return DIV(SUB(POW(k, e+1), POW(k,b)), k-1); 4fd800b3a8 2011-02-23 kinaba: } 4fd800b3a8 2011-02-23 kinaba: 4fd800b3a8 2011-02-23 kinaba: 4fd800b3a8 2011-02-23 kinaba: 4fd800b3a8 2011-02-23 kinaba: 4fd800b3a8 2011-02-23 kinaba: LL Cpascal(LL n, LL k) 4fd800b3a8 2011-02-23 kinaba: { 4fd800b3a8 2011-02-23 kinaba: vector< vector<LL> > c(n+1, vector<LL>(k+1)); 4fd800b3a8 2011-02-23 kinaba: for(LL nn=1; nn<=n; ++nn) 4fd800b3a8 2011-02-23 kinaba: for(LL kk=0; kk<=min(nn,k); ++kk) 4fd800b3a8 2011-02-23 kinaba: c[nn][kk] = kk==0 || kk==nn ? 1 : ADD(c[nn-1][kk-1], c[nn-1][kk]); 4fd800b3a8 2011-02-23 kinaba: return c[n][k]; 4fd800b3a8 2011-02-23 kinaba: } 4fd800b3a8 2011-02-23 kinaba: 4fd800b3a8 2011-02-23 kinaba: vector< vector<LL> > MATMUL(vector< vector<LL> >& a, vector< vector<LL> >& b) 4fd800b3a8 2011-02-23 kinaba: { 4fd800b3a8 2011-02-23 kinaba: int N = a.size(); 4fd800b3a8 2011-02-23 kinaba: vector< vector<LL> > c(N, vector<LL>(N)); 4fd800b3a8 2011-02-23 kinaba: for(int i=0; i<N; ++i) 4fd800b3a8 2011-02-23 kinaba: for(int j=0; j<N; ++j) 4fd800b3a8 2011-02-23 kinaba: for(int k=0; k<N; ++k) 4fd800b3a8 2011-02-23 kinaba: c[i][j] = ADD(c[i][j], MUL(a[i][k],b[k][j])); 4fd800b3a8 2011-02-23 kinaba: return c; 4fd800b3a8 2011-02-23 kinaba: } 4fd800b3a8 2011-02-23 kinaba: 4fd800b3a8 2011-02-23 kinaba: // works for non-prime MODVAL 4fd800b3a8 2011-02-23 kinaba: LL GEO(LL x_, LL e) // x^0 + x^1 + ... + x^e-1 4fd800b3a8 2011-02-23 kinaba: { 4fd800b3a8 2011-02-23 kinaba: vector< vector<LL> > v(2, vector<LL>(2)); 4fd800b3a8 2011-02-23 kinaba: vector< vector<LL> > x(2, vector<LL>(2)); 4fd800b3a8 2011-02-23 kinaba: v[0][0] = v[1][1] = 1; 4fd800b3a8 2011-02-23 kinaba: x[0][0] = x_; x[0][1] = 0; 4fd800b3a8 2011-02-23 kinaba: x[1][0] = 1 ; x[1][1] = 1; 4fd800b3a8 2011-02-23 kinaba: for(;e;x=MATMUL(x,x),e>>=1) 4fd800b3a8 2011-02-23 kinaba: if(e&1) 4fd800b3a8 2011-02-23 kinaba: v = MATMUL(v, x); 4fd800b3a8 2011-02-23 kinaba: return v[1][0]; 4fd800b3a8 2011-02-23 kinaba: } 4fd800b3a8 2011-02-23 kinaba: 4fd800b3a8 2011-02-23 kinaba: // works for non-prime MODVAL 4fd800b3a8 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) 4fd800b3a8 2011-02-23 kinaba: { 4fd800b3a8 2011-02-23 kinaba: vector< vector<LL> > v(3, vector<LL>(3)); 4fd800b3a8 2011-02-23 kinaba: vector< vector<LL> > x(3, vector<LL>(3)); 4fd800b3a8 2011-02-23 kinaba: v[0][0] = v[1][1] = v[2][2] = 1; 4fd800b3a8 2011-02-23 kinaba: x[0][0] = x_; x[0][1] = 0; x[0][2] = 0; 4fd800b3a8 2011-02-23 kinaba: x[1][0] = 1 ; x[1][1] = 1; x[1][2] = 0; 4fd800b3a8 2011-02-23 kinaba: x[2][0] = 0 ; x[2][1] = 1; x[2][2] = 1; 4fd800b3a8 2011-02-23 kinaba: e++; 4fd800b3a8 2011-02-23 kinaba: for(;e;x=MATMUL(x,x),e>>=1) 4fd800b3a8 2011-02-23 kinaba: if(e&1) 4fd800b3a8 2011-02-23 kinaba: v = MATMUL(v, x); 4fd800b3a8 2011-02-23 kinaba: return v[2][0]; 4fd800b3a8 2011-02-23 kinaba: }