ADDED   SRM/554/1A.cpp
Index: SRM/554/1A.cpp
==================================================================
--- SRM/554/1A.cpp
+++ SRM/554/1A.cpp
@@ -0,0 +1,98 @@
+#include <iostream>
+#include <sstream>
+#include <iomanip>
+#include <vector>
+#include <string>
+#include <map>
+#include <set>
+#include <algorithm>
+#include <numeric>
+#include <iterator>
+#include <functional>
+#include <complex>
+#include <queue>
+#include <stack>
+#include <cmath>
+#include <cassert>
+using namespace std;
+typedef long long LL;
+typedef long double LD;
+typedef complex<LD> CMP;
+
+class TheBrickTowerEasyDivOne { public:
+	int find(int redCount, int redHeight, int blueCount, int blueHeight)
+	{
+		if(redHeight != blueHeight)
+		{
+			int ans = 0;
+			ans += min(redCount, blueCount);
+			ans += min(redCount-1, blueCount) + 1;
+			ans += min(redCount, blueCount-1) + 1;
+			return ans;
+		}
+		return min(redCount,blueCount)*2 + (redCount != blueCount);
+	}
+};
+
+// BEGIN CUT HERE
+#include <ctime>
+double start_time; string timer()
+ { ostringstream os; os << " (" << int((clock()-start_time)/CLOCKS_PER_SEC*1000) << " msec)"; return os.str(); }
+template<typename T> ostream& operator<<(ostream& os, const vector<T>& v)
+ { os << "{ ";
+   for(typename vector<T>::const_iterator it=v.begin(); it!=v.end(); ++it)
+   os << '\"' << *it << '\"' << (it+1==v.end() ? "" : ", "); os << " }"; return os; }
+void verify_case(const int& Expected, const int& Received) {
+ bool ok = (Expected == Received);
+ if(ok) cerr << "PASSED" << timer() << endl;  else { cerr << "FAILED" << timer() << endl;
+ cerr << "\to: \"" << Expected << '\"' << endl << "\tx: \"" << Received << '\"' << endl; } }
+#define CASE(N) {cerr << "Test Case #" << N << "..." << flush; start_time=clock();
+#define END	 verify_case(_, TheBrickTowerEasyDivOne().find(redCount, redHeight, blueCount, blueHeight));}
+int main(){
+
+CASE(0)
+	int redCount = 1; 
+	int redHeight = 2; 
+	int blueCount = 3; 
+	int blueHeight = 4; 
+	int _ = 4; 
+END
+CASE(1)
+	int redCount = 4; 
+	int redHeight = 4; 
+	int blueCount = 4; 
+	int blueHeight = 7; 
+	int _ = 12; 
+END
+CASE(2)
+	int redCount = 7; 
+	int redHeight = 7; 
+	int blueCount = 4; 
+	int blueHeight = 4; 
+	int _ = 13; 
+END
+CASE(3)
+	int redCount = 47; 
+	int redHeight = 47; 
+	int blueCount = 47; 
+	int blueHeight = 47; 
+	int _ = 94; 
+END
+/*
+CASE(4)
+	int redCount = ; 
+	int redHeight = ; 
+	int blueCount = ; 
+	int blueHeight = ; 
+	int _ = ; 
+END
+CASE(5)
+	int redCount = ; 
+	int redHeight = ; 
+	int blueCount = ; 
+	int blueHeight = ; 
+	int _ = ; 
+END
+*/
+}
+// END CUT HERE

ADDED   SRM/554/1B.cpp
Index: SRM/554/1B.cpp
==================================================================
--- SRM/554/1B.cpp
+++ SRM/554/1B.cpp
@@ -0,0 +1,126 @@
+#include <iostream>
+#include <sstream>
+#include <iomanip>
+#include <vector>
+#include <string>
+#include <map>
+#include <set>
+#include <algorithm>
+#include <numeric>
+#include <iterator>
+#include <functional>
+#include <complex>
+#include <queue>
+#include <stack>
+#include <cmath>
+#include <cassert>
+using namespace std;
+typedef long long LL;
+typedef long double LD;
+typedef complex<LD> CMP;
+
+class TheBrickTowerMediumDivOne { public:
+	vector <int> find(vector<int> heights)
+	{
+		const int N = heights.size();
+		const int best = best_for(heights);
+		vector<bool> used(N, false);
+		vector<int> answer;
+		vector<int> current;
+		int current_score = 0;
+		for(int i=0; i<N; ++i)
+		{
+			for(int k=0; k<N; ++k) if(!used[k])
+			{
+				vector<int> hh;
+				for(int a=0; a<N; ++a) if(!used[a] && a!=k)
+					hh.push_back(heights[a]);
+				int ss = current_score;
+				if(!current.empty()) ss += max(current.back(), heights[k]);
+				ss += best_for(heights[k], hh);
+				if( ss == best ) {
+					answer.push_back(k);
+					if(!current.empty()) current_score += max(current.back(), heights[k]);
+					current.push_back(heights[k]);
+					used[k] = true;
+					break;
+				}
+			}
+		}
+		return answer;
+	}
+
+	int best_for(vector<int> h)
+	{
+		sort(h.begin(), h.end());
+		int best = 0;
+		for(int i=0; i+1<h.size(); ++i)
+			best += max(h[i], h[i+1]);
+		return best;
+	}
+
+	int best_for(int h0, vector<int> h)
+	{
+		if(h.empty())
+			return 0;
+
+		sort(h.begin(), h.end());
+		int best = max(h0, h[0]);
+		for(int i=0; i+1<h.size(); ++i)
+			best += max(h[i], h[i+1]);
+		return best;
+	}
+};
+
+// BEGIN CUT HERE
+#include <ctime>
+double start_time; string timer()
+ { ostringstream os; os << " (" << int((clock()-start_time)/CLOCKS_PER_SEC*1000) << " msec)"; return os.str(); }
+template<typename T> ostream& operator<<(ostream& os, const vector<T>& v)
+ { os << "{ ";
+   for(typename vector<T>::const_iterator it=v.begin(); it!=v.end(); ++it)
+   os << '\"' << *it << '\"' << (it+1==v.end() ? "" : ", "); os << " }"; return os; }
+void verify_case(const vector <int>& Expected, const vector <int>& Received) {
+ bool ok = (Expected == Received);
+ if(ok) cerr << "PASSED" << timer() << endl;  else { cerr << "FAILED" << timer() << endl;
+ cerr << "\to: " << Expected << endl << "\tx: " << Received << endl; } }
+#define CASE(N) {cerr << "Test Case #" << N << "..." << flush; start_time=clock();
+#define END	 verify_case(_, TheBrickTowerMediumDivOne().find(heights));}
+int main(){
+
+CASE(0)
+	int heights_[] = {4, 7, 5};
+	  vector <int> heights(heights_, heights_+sizeof(heights_)/sizeof(*heights_)); 
+	int __[] = {0, 2, 1 };
+	  vector <int> _(__, __+sizeof(__)/sizeof(*__)); 
+END
+CASE(1)
+	int heights_[] = {4, 4, 4, 4, 4, 4, 4};
+	  vector <int> heights(heights_, heights_+sizeof(heights_)/sizeof(*heights_)); 
+	int __[] = {0, 1, 2, 3, 4, 5, 6 };
+	  vector <int> _(__, __+sizeof(__)/sizeof(*__)); 
+END
+CASE(2)
+	int heights_[] = {2, 3, 3, 2};
+	  vector <int> heights(heights_, heights_+sizeof(heights_)/sizeof(*heights_)); 
+	int __[] = {0, 3, 1, 2 };
+	  vector <int> _(__, __+sizeof(__)/sizeof(*__)); 
+END
+CASE(3)
+	int heights_[] = {13, 32, 38, 25, 43, 47, 6};
+	  vector <int> heights(heights_, heights_+sizeof(heights_)/sizeof(*heights_)); 
+	int __[] = {0, 6, 3, 1, 2, 4, 5 };
+	  vector <int> _(__, __+sizeof(__)/sizeof(*__)); 
+END
+CASE(4)
+	int heights_[] = {28,36,27,33,30,18,16,23,39,31,4,40,41,41,30,43,45,16,5,26,35,19,2,40,27,3,37,46,5,34,22,25,19,6,45,6,45,36,40,13,29,25,10,23,6,3,21};
+	  vector <int> heights(heights_, heights_+sizeof(heights_)/sizeof(*heights_)); 
+	  vector <int> _;
+END
+CASE(5)
+	int heights_[] = {7,34,5,15,43,8,36,5,42,22,37,45,46,38,40,2,3,22,35,23,14,31,26,1,43,11,24,41,21,8,11,28,46,21,42,5,9,18,4,3,47,15,32,16,26,45,34};
+	  vector <int> heights(heights_, heights_+sizeof(heights_)/sizeof(*heights_)); 
+	  vector <int> _;
+END
+}
+// END CUT HERE

ADDED   SRM/554/1C.cpp
Index: SRM/554/1C.cpp
==================================================================
--- SRM/554/1C.cpp
+++ SRM/554/1C.cpp
@@ -0,0 +1,241 @@
+#include <iostream>
+#include <sstream>
+#include <iomanip>
+#include <vector>
+#include <string>
+#include <map>
+#include <set>
+#include <algorithm>
+#include <numeric>
+#include <iterator>
+#include <functional>
+#include <complex>
+#include <queue>
+#include <stack>
+#include <cmath>
+#include <cassert>
+using namespace std;
+typedef long long LL;
+typedef long double LD;
+typedef complex<LD> CMP;
+
+static const int MODVAL = 1234567891;
+struct mint
+{
+	unsigned val;
+	mint():val(0){}
+	mint(int      x):val(x%MODVAL) {}
+	mint(unsigned x):val(x%MODVAL) {}
+	mint(LL       x):val(x%MODVAL) {}
+};
+mint& operator+=(mint& x, mint y) { return x = x.val+y.val; }
+mint& operator-=(mint& x, mint y) { return x = x.val-y.val+MODVAL; }
+mint& operator*=(mint& x, mint y) { return x = LL(x.val)*y.val; }
+mint operator+(mint x, mint y) { return x+=y; }
+mint operator-(mint x, mint y) { return x-=y; }
+mint operator*(mint x, mint y) { return x*=y; }
+
+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< 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;
+}
+
+class TheBrickTowerHardDivOne { public:
+	int find(int C, int K, long long H)
+	{
+		vector< vector<int> > patterns;
+		{
+			vector<int> p;
+			vector<bool> used(4);
+			generate_patterns(true, 4, used, p, 0, patterns);
+		}
+		const vector<int> compress = calc_compress_table(patterns).first;
+		const vector<bool> is_repr = calc_compress_table(patterns).second;
+
+		const int PN = patterns.size();
+		const int CPN = *max_element(compress.begin(), compress.end())+1;
+		const int D = CPN * (K+1);
+
+		vector< vector<mint> > m(D+1, vector<mint>(D+1));
+		for(int pa=0; pa<PN; ++pa) if(is_repr[pa])
+		for(int pb=0; pb<PN; ++pb) {
+			vector<mint> tt = transition(patterns[pa], patterns[pb], K, C);
+			for(int ka=0; ka<=K; ++ka)
+			for(int kt=0; kt<=K; ++kt) {
+				int kb = ka + kt + internal(patterns[pb]);
+				if(kb <= K)
+					m[compress[pb]*(K+1)+kb][compress[pa]*(K+1)+ka] += tt[kt];
+			}
+		}
+		for(int x=0; x<D+1; ++x)
+			m[D][x] = 1;
+
+		vector<mint> v(D+1);
+		for(int pa=0; pa<PN; ++pa) {
+			int ka = internal(patterns[pa]);
+			if(ka <= K)
+				v[compress[pa]*(K+1)+ka] += counting(patterns[pa], C);
+		}
+
+		v = MATMUL(MATPOW(m, H-1), v);
+		return accumulate(v.begin(), v.end(), mint(0)).val;
+	}
+
+	// TODO: auto generate.
+	pair< vector<int>, vector<bool> >
+		calc_compress_table(const vector< vector<int> >& patterns)
+	{
+		// {0}         00,00
+		// {1,2,5,9}   00,01
+		// {3,6}       00,11
+		// {4,7,12,13} 00,12
+		// {8}         01,10
+		// {10,11}     01,12
+		// {14}        01,23
+		vector<int> repr(patterns.size());
+		vector<bool> is_repr(patterns.size(), false);
+		is_repr[0] = true; repr[0] = 0;
+		is_repr[1] = true; repr[1] = repr[2] = repr[5] = repr[9] = 1;
+		is_repr[3] = true; repr[3] = repr[6] = 2;
+		is_repr[4] = true; repr[4] = repr[7] = repr[12] = repr[13] = 3;
+		is_repr[8] = true; repr[8] = 4;
+		is_repr[10] = true; repr[10] = repr[11] = 5;
+		is_repr[14] = true; repr[14] = 6;
+		return make_pair(repr, is_repr);
+	}
+
+	mint counting(const vector<int>& p, int C)
+	{
+		return counting(set<int>(p.begin(), p.end()).size(), C);
+	}
+
+	mint counting(int N, int C)
+	{
+		if(C<0)
+			return 0;
+		mint v = 1;
+		for(int i=0; i<N; ++i)
+			v *= (C-i);
+		return v;
+	}
+
+	int internal(const vector<int>& p)
+	{
+		return (p[0]==p[1])+(p[1]==p[3])+(p[3]==p[2])+(p[2]==p[0]);
+	}
+
+	vector<mint> transition(const vector<int>& pa, const vector<int>& pb, int K, int C)
+	{
+		const int AN = set<int>(pa.begin(), pa.end()).size();
+		const int BN = set<int>(pb.begin(), pb.end()).size();
+
+		vector< vector<int> > patterns;
+		{
+			vector<int> p;
+			vector<bool> used(AN+BN, false);
+			generate_patterns(false, BN, used, p, AN, patterns);
+		}
+
+		vector<mint> answer(K+1);
+		for(int i=0; i<patterns.size(); ++i) {
+			const vector<int>& p = patterns[i];
+			int k = 0;
+			for(int x=0; x<4; ++x) k += (pa[x]==p[pb[x]]);
+			
+			if(k <= K) {
+				int N = max(*max_element(p.begin(), p.end()) - (AN-1), 0);
+				answer[k] += counting(N, C-AN);
+			}
+		}
+		return answer;
+	}
+
+	void generate_patterns(bool ALLOW_DUP, int LEN, vector<bool>& used,
+	                       vector<int>& p, int next, vector< vector<int> >& patterns)
+	{
+		if(p.size() == LEN)
+			patterns.push_back(p);
+		else
+			for(int v=0; v<=next; ++v) if(ALLOW_DUP || !used[v]) {
+				p.push_back(v);
+				used[v] = true;
+				generate_patterns(ALLOW_DUP, LEN, used, p, max(v+1, next), patterns);
+				used[v] = false;
+				p.pop_back();
+			}
+	}
+};
+
+// BEGIN CUT HERE
+#include <ctime>
+double start_time; string timer()
+ { ostringstream os; os << " (" << int((clock()-start_time)/CLOCKS_PER_SEC*1000) << " msec)"; return os.str(); }
+template<typename T> ostream& operator<<(ostream& os, const vector<T>& v)
+ { os << "{ ";
+   for(typename vector<T>::const_iterator it=v.begin(); it!=v.end(); ++it)
+   os << '\"' << *it << '\"' << (it+1==v.end() ? "" : ", "); os << " }"; return os; }
+void verify_case(const int& Expected, const int& Received) {
+ bool ok = (Expected == Received);
+ if(ok) cerr << "PASSED" << timer() << endl;  else { cerr << "FAILED" << timer() << endl;
+ cerr << "\to: \"" << Expected << '\"' << endl << "\tx: \"" << Received << '\"' << endl; } }
+#define CASE(N) {cerr << "Test Case #" << N << "..." << flush; start_time=clock();
+#define END	 verify_case(_, TheBrickTowerHardDivOne().find(C, K, H));}
+int main(){
+CASE(0)
+	int C = 2; 
+	int K = 0; 
+	long long H = 2LL; 
+	int _ = 4; 
+END
+CASE(2)
+	int C = 2; 
+	int K = 3; 
+	long long H = 1LL; 
+	int _ = 14; 
+END
+CASE(3)
+	int C = 4; 
+	int K = 7; 
+	long long H = 47LL; 
+	int _ = 1008981254; 
+END
+CASE(4)
+	int C = 4747; 
+	int K = 7; 
+	long long H = 474747474747474747LL; 
+	int _ = 473182063; 
+END
+}
+// END CUT HERE

Index: lib/numeric/modArith.cpp
==================================================================
--- lib/numeric/modArith.cpp
+++ lib/numeric/modArith.cpp
@@ -5,32 +5,37 @@
 // Verified by
 //   - TCO10 R3 LV3
 //   - SRM 545 LV2
 //-------------------------------------------------------------
 
-static const int MODVAL = 1000000007;
+static const unsigned MODVAL = 1000000007;
 struct mint
 {
-	int val;
+	unsigned val;
 	mint():val(0){}
-	mint(int    x):val(x%MODVAL) {}
-	mint(size_t x):val(x%MODVAL) {}
-	mint(LL     x):val(x%MODVAL) {}
+	mint(int      x):val(x%MODVAL) {}
+	mint(unsigned x):val(x%MODVAL) {}
+	mint(LL       x):val(x%MODVAL) {}
 };
 mint& operator+=(mint& x, mint y) { return x = x.val+y.val; }
 mint& operator-=(mint& x, mint y) { return x = x.val-y.val+MODVAL; }
 mint& operator*=(mint& x, mint y) { return x = LL(x.val)*y.val; }
-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; }
 mint operator-(mint x, mint y) { return x-=y; }
 mint operator*(mint x, mint y) { return x*=y; }
+
+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]; }
 mint C(LL n, LL k) { return k<0 || n<k ? 0 : FAC(n) / (FAC(k) * FAC(n-k)); }
-vector< vector<mint> > CP_; // Pascal's triangle: if O(1) nCk is needed.
+
+// nCk :: O(1) time, O(n^2) space.
+vector< vector<mint> > CP_;
 mint C(LL n, LL k) {
 	while( CP_.size() <= n ) {
 		int nn = CP_.size();
 		CP_.push_back(vector<mint>(nn+1,1));
 		for(int kk=1; kk<nn; ++kk)