34 int nn = CP_.size(); 34 int nn = CP_.size(); 35 CP_.push_back(vector<mint>(nn+1,1)); 35 CP_.push_back(vector<mint>(nn+1,1)); 36 for(int kk=1; kk<nn; ++kk) 36 for(int kk=1; kk<nn; ++kk) 37 CP_[nn][kk] = CP_[nn-1][kk-1] + CP_[nn-1][kk]; 37 CP_[nn][kk] = CP_[nn-1][kk-1] + CP_[nn-1][kk]; 38 } 38 } 39 return k<0 || n<k ? 0 : CP_[n][k]; 39 return k<0 || n<k ? 0 : CP_[n][k]; 40 } 40 } > 41 > 42 template<typename T> > 43 vector<T> MATMUL(const vector< vector<T> >& a, const vector<T>& v) > 44 { > 45 int N = a.size(); > 46 vector<T> u(N); > 47 for(int i=0; i<N; ++i) > 48 for(int j=0; j<N; ++j) > 49 u[i] += a[i][j]*v[j]; > 50 return u; > 51 } > 52 > 53 template<typename T> > 54 vector< vector<T> > MATMUL(const vector< vector<T> >& a, const vector< vector<T> > 55 { > 56 int N = a.size(); > 57 vector< vector<T> > c(N, vector<T>(N)); > 58 for(int i=0; i<N; ++i) > 59 for(int j=0; j<N; ++j) > 60 for(int k=0; k<N; ++k) > 61 c[i][j] += a[i][k]*b[k][j]; > 62 return c; > 63 } > 64 > 65 template<typename T> > 66 vector< vector<T> > MATPOW(vector< vector<T> > a, LL e) > 67 { > 68 int N = a.size(); > 69 vector< vector<T> > c(N, vector<T>(N)); > 70 for(int i=0; i<N; ++i) c[i][i] = 1; > 71 for(; e; e>>=1) { > 72 if(e&1) > 73 c = MATMUL(c, a); > 74 a = MATMUL(a, a); > 75 } > 76 return c; > 77 } 41 78 42 /* 79 /* 43 // MODVAL must be a prime!! 80 // MODVAL must be a prime!! 44 LL GSS(LL k, LL b, LL e) // k^b + k^b+1 + ... + k^e 81 LL GSS(LL k, LL b, LL e) // k^b + k^b+1 + ... + k^e 45 { 82 { 46 if( b > e ) return 0; 83 if( b > e ) return 0; 47 if( k <= 1 ) return k*(e-b+1); 84 if( k <= 1 ) return k*(e-b+1); 48 return DIV(SUB(POW(k, e+1), POW(k,b)), k-1); 85 return DIV(SUB(POW(k, e+1), POW(k,b)), k-1); 49 } 86 } 50 87 51 LL Cpascal(LL n, LL k) < 52 { < 53 vector< vector<LL> > c(n+1, vector<LL>(k+1)); < 54 for(LL nn=1; nn<=n; ++nn) < 55 for(LL kk=0; kk<=min(nn,k); ++kk) < 56 c[nn][kk] = kk==0 || kk==nn ? 1 : ADD(c[nn-1][kk-1], c[n < 57 return c[n][k]; < 58 } < 59 < 60 vector< vector<LL> > MATMUL(vector< vector<LL> >& a, vector< vector<LL> >& b) < 61 { < 62 int N = a.size(); < 63 vector< vector<LL> > c(N, vector<LL>(N)); < 64 for(int i=0; i<N; ++i) < 65 for(int j=0; j<N; ++j) < 66 for(int k=0; k<N; ++k) < 67 c[i][j] = ADD(c[i][j], MUL(a[i][k],b[k][j])); < 68 return c; < 69 } < 70 < 71 // works for non-prime MODVAL 88 // works for non-prime MODVAL 72 LL GEO(LL x_, LL e) // x^0 + x^1 + ... + x^e-1 89 LL GEO(LL x_, LL e) // x^0 + x^1 + ... + x^e-1 73 { 90 { 74 vector< vector<LL> > v(2, vector<LL>(2)); 91 vector< vector<LL> > v(2, vector<LL>(2)); 75 vector< vector<LL> > x(2, vector<LL>(2)); 92 vector< vector<LL> > x(2, vector<LL>(2)); 76 v[0][0] = v[1][1] = 1; 93 v[0][0] = v[1][1] = 1; 77 x[0][0] = x_; x[0][1] = 0; 94 x[0][0] = x_; x[0][1] = 0;