Index: lib/numeric/matrix.cpp
==================================================================
--- lib/numeric/matrix.cpp
+++ lib/numeric/matrix.cpp
@@ -1,55 +1,118 @@
-
 //-------------------------------------------------------------
 // Matrix Operations
 //
 // Verified by
 //   - SRM342 Div1 LV3
 //   - SRM341 Div1 LV3
 //   - SRM338 Div1 LV2
 //-------------------------------------------------------------
 
-vector<LL> vMul( const vector< vector<LL> >& A, const vector<LL>& B )
+template<typename T>
+vector<T> MATMUL(const vector< vector<T> >& a, const vector<T>& v)
+{
+	int N = a.size();
+	vector<T> u(N);
+	for(int i=0; i<N; ++i)
+	for(int j=0; j<N; ++j)
+		u[i] += a[i][j]*v[j];
+	return u;
+}
+
+template<typename T>
+vector< vector<T> > MATMUL(const vector< vector<T> >& a, const vector< vector<T> >& b)
+{
+	int N = a.size();
+	vector< vector<T> > c(N, vector<T>(N));
+	for(int i=0; i<N; ++i)
+	for(int j=0; j<N; ++j)
+	for(int k=0; k<N; ++k)
+		c[i][j] += a[i][k]*b[k][j];
+	return c;
+}
+
+template<typename T>
+vector<T> MATPOWMUL(vector< vector<T> > a, LL e, vector<T> v)
 {
-	const int n = A.size();
-
-	vector<LL> C(n);
-	for(int i=0; i<n; ++i)
-		for(int j=0; j<n; ++j)
-			C[i] += A[i][j] * B[j];
-	return C;
+	for(; e; e>>=1, a=MATMUL(a,a))
+		if(e&1)
+			v = MATMUL(a, v);
+	return v;
 }
 
-vector< vector<LL> > mMul( const vector< vector<LL> >& A, const vector< vector<LL> >& B )
+template<typename T>
+vector< vector<T> > MATPOW(vector< vector<T> > a, LL e)
 {
-	const int n = A.size();
-
-	vector< vector<LL> > C(n, vector<LL>(n));
-	for(int i=0; i<n; ++i)
-		for(int j=0; j<n; ++j) {
-			LL Cij = 0;
-			for(int k=0; k<n; ++k)
-				Cij += A[i][k] * B[k][j];
-			C[i][j] = Cij;
-		}
-	return C;
+	int N = a.size();
+	vector< vector<T> > c(N, vector<T>(N));
+	for(int i=0; i<N; ++i) c[i][i] = 1;
+	for(; e; e>>=1) {
+		if(e&1)
+			c = MATMUL(c, a);
+		a = MATMUL(a, a);
+	}
+	return c;
 }
 
-vector< vector<LL> > mPow( vector< vector<LL> > M, LL t ) // t>0
-{
-	vector< vector<LL> > R;
-	for(; t; t>>=1, M=mMul(M,M))
-		if( t&1 )
-			R = (R.empty() ? M : mMul(R,M));
-	return R;
-}
+/****************************************************************************
+   NOTES ON THE USAGE
+ ****************************************************************************
+
+A[to][from] = #transition
+--------------------------------------------
+
+     frm
+   (. . .) (q0)
+to (. . .) (q1)
+   (. . .) (q2)
+
+
+
+ƒ° qi
+--------------------------------------------
+
+     frm
+   (. . . 0) (.)
+to (. . . 0) (.)
+   (. . . 0) (.)
+   (0 1 0 1) (0)
+      ^
+      i
+Then,
+  q[i]@0                 = (A^1 v)[$-1]
+  q[i]@0 + q[i]@1        = (A^2 v)[$-1]
+  q[i]@0 + ... + q[i]@k  = (A^(k+1) v)[$-1]
+
+
+
+ƒ° q_all
+--------------------------------------------
+
+     frm
+   (. . . 0) (.)
+to (. . . 0) (.)
+   (. . . 0) (.)
+   (1 1 1 1) (0)
+
+Then,
+  ƒ°q @0                 = (A^1 v)[$-1]
+  ƒ°q @0 + ... + ƒ°q @k  = (A^(k+1) v)[$-1]
+
+
+
+
+
+Application: x^0 + x^1 + ... + x^e-1
+--------------------------------------------
+
+  (x 0)^e (1)  = (x^e                    )
+  (1 1)   (0)    (x^0 + x^1 + ... + x^e-1)
+
+
 
-vector< vector<LL> > mAdd( const vector< vector<LL> >& A, const vector< vector<LL> >& B )
-{
-	const int n = A.size();
+Application: e x^0 + (e-1) x^1 + ... + 1 x^e-1
+--------------------------------------------
 
-	vector< vector<LL> > C(n, vector<LL>(n));
-	for(int i=0; i<n; ++i)
-		for(int j=0; j<n; ++j)
-			C[i][j] = A[i][j] + B[i][j];
-	return C;
-}
+  (x 0 0)^e (x)  = (x^e+1                                       )
+  (1 1 0)   (1)    (x^0 + x^1 + ... + x^e                       )
+  (0 1 1)   (0)    (ƒ° of �ª = e x^0 + (e-1) x^1 + ... + 1 x^e-1)
+*/

Index: lib/numeric/modArith.cpp
==================================================================
--- lib/numeric/modArith.cpp
+++ lib/numeric/modArith.cpp
@@ -2,11 +2,12 @@
 //-------------------------------------------------------------
 // Modulo Arithmetics
 //
 // Verified by
 //   - TCO10 R3 LV3
-//   - SRM 545 LV2
+//   - SRM 545 Div1 LV2
+//   - SRM 554 Div1 LV3
 //-------------------------------------------------------------
 
 static const unsigned MODVAL = 1000000007;
 struct mint
 {
@@ -25,13 +26,14 @@
 
 mint POW(mint x, LL e) { mint v=1; for(;e;x*=x,e>>=1) if(e&1) v*=x; return v; }
 mint& operator/=(mint& x, mint y) { return x *= POW(y, MODVAL-2); }
 mint operator/(mint x, mint y) { return x/=y; }
 
-// nCk :: O(log MODVAL) time, O(n) space.
 vector<mint> FAC_(1,1);
 mint FAC(LL n) { while( FAC_.size()<=n ) FAC_.push_back( FAC_.back()*FAC_.size() ); return FAC_[n]; }
+
+// nCk :: O(log MODVAL) time, O(n) space.
 mint C(LL n, LL k) { return k<0 || n<k ? 0 : FAC(n) / (FAC(k) * FAC(n-k)); }
 
 // nCk :: O(1) time, O(n^2) space.
 vector< vector<mint> > CP_;
 mint C(LL n, LL k) {
@@ -42,81 +44,12 @@
 			CP_[nn][kk] = CP_[nn-1][kk-1] + CP_[nn-1][kk];
 	}
 	return k<0 || n<k ? 0 : CP_[n][k];
 }
 
-template<typename T>
-vector<T> MATMUL(const vector< vector<T> >& a, const vector<T>& v)
+// O(log MODVAL), MODVAL must be prime: k^b + k^b+1 + ... + k^e
+mint GSS(mint k, LL b, LL e) 
 {
-	int N = a.size();
-	vector<T> u(N);
-	for(int i=0; i<N; ++i)
-	for(int j=0; j<N; ++j)
-		u[i] += a[i][j]*v[j];
-	return u;
+	if( b > e ) return 0;
+	if( k.val <= 1 ) return k*(e-b+1);
+	return (POW(k, e+1) - POW(k,b)) / (k-1);
 }
-
-template<typename T>
-vector< vector<T> > MATMUL(const vector< vector<T> >& a, const vector< vector<T> >& b)
-{
-	int N = a.size();
-	vector< vector<T> > c(N, vector<T>(N));
-	for(int i=0; i<N; ++i)
-	for(int j=0; j<N; ++j)
-	for(int k=0; k<N; ++k)
-		c[i][j] += a[i][k]*b[k][j];
-	return c;
-}
-
-template<typename T>
-vector< vector<T> > MATPOW(vector< vector<T> > a, LL e)
-{
-	int N = a.size();
-	vector< vector<T> > c(N, vector<T>(N));
-	for(int i=0; i<N; ++i) c[i][i] = 1;
-	for(; e; e>>=1) {
-		if(e&1)
-			c = MATMUL(c, a);
-		a = MATMUL(a, a);
-	}
-	return c;
-}
-
-/*
-// MODVAL must be a prime!!
-LL GSS(LL k, LL b, LL e) // k^b + k^b+1 + ... + k^e
-{
-	if( b >  e ) return 0;
-	if( k <= 1 ) return k*(e-b+1);
-	return DIV(SUB(POW(k, e+1), POW(k,b)), k-1);
-}
-
-// works for non-prime MODVAL
-LL GEO(LL x_, LL e) // x^0 + x^1 + ... + x^e-1
-{
-   vector< vector<LL> > v(2, vector<LL>(2));
-   vector< vector<LL> > x(2, vector<LL>(2));
-   v[0][0] = v[1][1] = 1;
-   x[0][0] = x_; x[0][1] = 0;
-   x[1][0] = 1 ; x[1][1] = 1;
-   for(;e;x=MATMUL(x,x),e>>=1)
-      if(e&1)
-         v = MATMUL(v, x);
-   return v[1][0];
-}
-
-// works for non-prime MODVAL
-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)
-{
-   vector< vector<LL> > v(3, vector<LL>(3));
-   vector< vector<LL> > x(3, vector<LL>(3));
-   v[0][0] = v[1][1] = v[2][2] = 1;
-   x[0][0] = x_; x[0][1] = 0; x[0][2] = 0;
-   x[1][0] = 1 ; x[1][1] = 1; x[1][2] = 0;
-   x[2][0] = 0 ; x[2][1] = 1; x[2][2] = 1;
-   e++;
-   for(;e;x=MATMUL(x,x),e>>=1)
-      if(e&1)
-         v = MATMUL(v, x);
-   return v[2][0];
-}
-*/