SRM 636

https://community.topcoder.com/stat?c=round_overview&rd=16079

Easy
累積和で前処理しておけばいいだけです。

struct ChocolateDividingEasy {
  vector<string> chocolate;
  int S[2510][2510];
  int findBest(vector<string> _chocolate) {
    chocolate = _chocolate;
	int n = sz(chocolate), m = sz(chocolate[0]);
	rep(i, 0, n) {
		rep(j, 0, m) {
			S[i + 1][j + 1] = S[i + 1][j] + S[i][j + 1] - S[i][j] + chocolate[i][j] - '0';
		}
	}
	int res = 0;
	for(int i = 1; i < n; i++) {
		for(int j = i + 1; j < n; j++) {
			for(int k = 1; k < m; k++) {
				for(int l = k + 1; l < m; l++) {
					vector<pi> x(3), y(3);
					x[0] = pi(0, i); y[0] = pi(0, k);
					x[1] = pi(i, j); y[1] = pi(k, l);
					x[2] = pi(j, n); y[2] = pi(l, m);
					vector<int> vec;
					rep(p, 0, 3) {
						rep(q, 0, 3) {
							int x1 = x[p].fst, x2 = x[p].sec;
							int y1 = y[q].fst, y2 = y[q].sec;
							vec.pb(S[x2][y2] - S[x2][y1] - S[x1][y2] + S[x1][y1]);
						}
					}
					sort(all(vec));
					MAX(res, vec[0]);
				}
			}
		}
	}
    return res;
  }
};

Medium
期待値について理解が深まりました。
よく観察すると、二重辺の数=連結成分の数になることがわかります。なぜなら構築したグラフはさきにループを含むDAGになり、またループは3点以上で構成できないからです。
よって期待値の加法定理より、求める値は、
Σ(e∈E)f(e)となります。ここでEは二点のすべてのとり方の集合、fはその二点が二重辺になる確率です。
各eにつきO(N^2)でfが求められるので、全体でO(N^6)となり、N<=20なので間に合います。でもTopCoderは最適化オプションつけないでコンパイルしているみたいで実行時間結構遅かった。

struct ClosestRabbit {
  vector<string> board;
  int r;
  double getExpected(vector<string> _board, int _r) {
    board = _board, r = _r;
	int n = sz(board), m = sz(board[0]);
	int cnt = 0;
	rep(i, 0, n) 
		rep(j, 0, m) 
			if(board[i][j] == '.') cnt++;
	double ans = 0.0;
	rep(i, 0, n) {
		rep(j, 0, m) {
			rep(k, 0, n) {
				rep(l, 0, m) {
					if(board[i][j] == '#' || board[k][l] == '#') continue;
					if(i == k && j == l) continue;
					double tmp = 1.0 * r * (r - 1) / (cnt * (cnt - 1));
					double dist = sqrt(1.0 * (i - k) * (i - k) + (j - l) * (j - l));
					int cnta = 0;
					rep(x, 0, n) {
						rep(y, 0, m) {
							if(board[x][y] == '#') continue;
							double d1 = sqrt(1.0 * (i - x) * (i - x) + (j - y) * (j - y));
							double d2 = sqrt(1.0 * (k - x) * (k - x) + (l - y) * (l - y));
							if(dist > d1 + eps) continue;
							if(dist > d2 + eps) continue;
							if(d1 - eps < dist && dist < d1 + eps) {
								if(x < k) continue;
								else if(x == k) {
									if(y < l) continue;
								}
							}
							if(d2 - eps < dist && dist < d2 + eps) {
								if(x < i) continue;
								else if(x == i) {
									if(y < j) continue;
								}
							}
							cnta++;
						}
					}
					if(cnta < r - 2) continue;
					rep(t, 0, r - 2) {
						tmp = tmp * 1.0 * (cnta - t) / (cnt - 2 - t);
					}
					ans += tmp;
				}
			}
		}
	}
    return ans / 2;
  }
};

数え上げDPについて

まずこの問題を考えてみましょう。
文字列S,Tが与えられる。最長共通部分列の長さを求めよ。

これは蟻本にあるようにdp[i][j]:=Sのi字目とTのj字目までの最長共通部分列の長さとすれば良くて
dp[i][j] = max(dp[i - 1][j], dp[i][j - 1]) (S[i]!=T[j])
dp[i][j] = dp[i - 1][j - 1] + 1 (S[i]==T[j])
でdp tableを更新すれば解くことができます。

では次の問題
文字列S,Tが与えられる。共通部分列の数を求めよ。
これもdp[i][j]:=Sのi字目とTのj字目までの共通部分列の数として
dp[i][j] = dp[i - 1][j] + dp[i][j - 1] (S[i]!=T[j])
dp[i][j] = dp[i - 1][j - 1] (S[i]==T[j])
で更新すればいいように思えます。

でもこれではだめで、例えば
S="abc"
T="adc"を考えると
acという共通部分列がdp[3][2]とdp[2][3]から遷移する2通りがあるのでダブルカウントしてしまっていることがわかります。

このように数え上げDPはダブルカウントに気を付ける必要があり、普通のDPよりも難易度が増していることがわかると思います。

それでいろいろ考えていたら、DPってそもそもなんだ?状態とか遷移とかってなんだ?という話になったので、
自分なりに整理してみました。


DPの本質は探索にあって、探索はある操作を行うことで進行します。

操作の定義は簡単で、ある対象に対して何かしらの作用を加えることです。また操作を何回か行って得られたものを結果と定義します。

この際、DPを行える探索は3つの性質を満たしていなければなりません。

1.DAGである。すなわち、有限回の操作で完了する。
期待値DPやゲームDPだと自己ループなどがあったりしますが、結局DAGに帰着できます。

2.求めたいものがある結果に対応している。
つまり求めたいものが探索で網羅できているということです。

3.状態というものがあって、操作の結果はただ一つの状態に属していなければならない。
状態は好きなもので良くて、例えば上の例ではSのi字目とTのj字目でしたが、別にSのi字目だけでも状態になりえます。
ただ一つの状態に属していなければならない、は状態同士が互いに素ということですが、自然な状態を持ってくれば勝手に互いに素になってくれるのでそんな気にしなくていいです。
状態から状態への移り変わりを遷移と定義すると、DPはこの遷移を理解するのがポイントになることがわかると思います。
さっきSのi字目だけを状態にすると、遷移がややこしいことになってDPしにくいです。また1つの状態に求めたいものと求めたくないものが混在していると意味が無いです。このようになんでも状態にはできますが、計算しやすいものは限られています。

ここまでが普通のDPの満たすべき性質ですが、数え上げDPは2を

2'.求めたいものはただひとつの結果と対応する

に変更すればいいです。

例えば最長共通部分列を求めるときは操作を
Sのi字目とTのj字目を見て、等しかったらi+1,j+1に移動して、そうでなければi+1,jまたはi,j+1に移動する
と定義できて、それで全共通文字列を網羅できます。しかし、i,jが異なっていて、i+1,j+1に移動する方法を考えた時、i+1,jに先に移動するか、i,j+1に先に移動するかで求めているものは変わりませんが、結果は2つになっています。よってこれでは数え上げDPができないことがわかります。

では共通部分列の数を求めるときは
i,j以降でS[i']==T[j']となるi',j'に移動するにすれば良いです。計算量はO(|S|^2|T|^2)となりますが、累積和を使えばO(|S||T|)になります。

上の例からわかるように2'は
任意の2つの操作を入れ替えると求めているものが変わる。
に言い換えられます。

このように数え上げDPは遷移が少し複雑になります。計算量を落とす際は累積和を使えばよさそうですね。

一般的な数え上げでも、この操作から考える方法でもれなくだぶりなく数え上げることができます。


二週間ずっと悩み続けた…

SRM 371

Topcoderの環境をUbuntuで作ったので試してみました。

昔のSRMだとあまり考察がなくて解きやすく、実装力をつけるのにもってこいですね。

Easy
一周するとwidth-2, height-2の場合に帰着できて、あとはwidth==2or1の場合分けをやる。

struct SpiralRoute {
  int width;
  int length;
  pi loop(int a, int b) {
	  if(a == 1) return pi(0, b - 1);
	  if(b == 1) return pi(a - 1, 0);
	  if(a == 2 || b == 2) return pi(0, 1);
	  pi res = loop(a - 2, b - 2);
	  res.fst++; res.sec++;
	  return res;
  }
  vector<int> thronePosition(int _width, int _length) {
    width = _width, length = _length;
	pi p = loop(width, length);
	vector<int> res;
	res.pb(p.fst); res.pb(p.sec);
    return res;
  }
};

Middle
最小費用流流すだけ。

struct ChessMatchup {
  vector<int> us;
  vector<int> them;
  int maximumScore(vector<int> _us, vector<int> _them) {
    us = _us, them = _them;
	int n = sz(us);
	MCF::init(2 * n + 2);
	int s = 2 * n, t = 2 * n + 1;
	rep(i, 0, n) {
		MCF::add_edge(s, i, 1, 0);
		MCF::add_edge(i + n, t, 1, 0);
	}
	rep(i, 0, n) {
		rep(j, 0, n) {
			if(us[i] > them[j]) MCF::add_edge(i, j + n, 1, -2);
			else if(us[i] == them[j]) MCF::add_edge(i, j + n, 1, -1);
			else MCF::add_edge(i, j + n, 1, 0);
		}
	}
    return -MCF::get(s, t, n);
  }
};

Hard
連結成分ごとにトポロジカルソートの数をbitdpで求めてあげればいいだけなのだが、かなりてこずった。

struct RaceOrdering {
  int n;
  vector<int> first;
  vector<int> second;
  vector<int> G[50];
  vector<int> C[50];
  int cnt[50];
  int pre[50];
  bool used[50];
  ll dp[(1 << 20)];

  void dfs(int v, int k) {
	  used[v] = true;
	  C[k].pb(v);
	  rep(i, 0, sz(G[v])) {
		  int n = G[v][i];
		  if(!used[n]) dfs(n, k);
	  }
  }

  int countOrders(int _n, vector<int> _first, vector<int> _second) {
    n = _n, first = _first, second = _second;
	C_init(n + 10);
	int m = sz(first);
	rep(i, 0, m) {
		G[second[i]].pb(first[i]);
		G[first[i]].pb(second[i]);
		cnt[second[i]]++;
	}
	rep(i, 0, n) debug(i, cnt[i]);
	int k = 0;
	rep(i, 0, n) {
		if(!used[i]) dfs(i, k++);
	}
	rep(i, 0, n) G[i].clear();
	rep(i, 0, m) {
		G[first[i]].pb(second[i]);
	}
	memset(used, 0, sizeof(used));
	rep(i, 0, k) {
		while(true) {
			bool found = false;
			rep(j, 0, sz(C[i])) {
				if(cnt[C[i][j]]) found = true;
			}
			if(!found) break;
			found = false;
			rep(j, 0, sz(C[i])) {
				int n = C[i][j];
				if(cnt[n] == 0 && !used[n]) {
					found = true;
					used[n] = true;
					int at = find(all(C[i]), n) - C[i].begin();
					rep(l, 0, sz(G[n])) {
						int nn = G[n][l];
						cnt[nn]--;
						pre[nn] |= (1 << at);
					}
				}
			}
			if(!found) return 0;
		}
	}
	ll ans = 1;
	rep(i, 0, n) debug(bitset<4>(pre[i]));
	rep(l, 0, k) {
		memset(dp, 0, sizeof(dp));
		dp[0] = 1;
		int nn = sz(C[l]);
		rep(bit, 0, (1 << nn)) {
			debug(bitset<4>(bit), dp[bit]);
			rep(i, 0, nn) {
				if(bit & (1 << i)) continue;
				if((pre[C[l][i]] & bit) == pre[C[l][i]]) {
					ADD(dp[bit | (1 << i)], dp[bit]);
				}
			}
		}
		MUL(ans, dp[(1 << nn) - 1]);
	}
	rep(l, 0, k) {
		MUL(ans, fiv[sz(C[l])]);
	}
	MUL(ans, fac[n]);
    return ans;
  }
};

ブログの更新が空いちゃって良くなかったね。CodeFestival予選CのEFをいろいろ考えてたんですが…。

AtCoder Regular Contest 070E: NarrowRectangles

http://arc070.contest.atcoder.jp/tasks/arc070_c

関数dpとかいうやつ。僕のいうところの「愚直に値を持つ」dpと同じです。
愚直にx座標すべて持ってみることを考えて、実は本当に持つ必要があるものは少なくて、更新も簡単にできる、みたいな感じです。
|x-c|を足すと、x=cで傾きが2変わることを失念してWA生やしました。
あとsetで最後の要素を消すのはS.erase(--S.end())です。

dpというとこのパターンと、考察して状態量を減らすパターンがあるのかな?前者はなんとかなりそうだけど、後者は一般的な解き方がなくて難しそう。

int N;
multiset<ll> L, R;
int A[MAX_N], B[MAX_N];

void solve() {
	cin >> N;
	rep(i, 0, N) cin >> A[i] >> B[i];
	ll a = A[0], b = A[0];
	L.insert(A[0]); R.insert(A[0]);
	ll ans = 0;
	rep(i, 1, N) {
		a -= B[i] - A[i];
		b += B[i - 1] - A[i - 1];
		ll l = *L.rbegin(), r = *R.begin();
		if(a <= A[i] && A[i] <= b) {
			L.insert(A[i] - (a - l));
			R.insert(A[i] - (b - r));
			a = A[i]; b = A[i];
		}
		else if(A[i] < a) {
			ans += abs(A[i] - a);
			L.insert(A[i] - (a - l));
			L.insert(A[i] - (a - l));
			ll le = *L.rbegin() + (a - l);
			L.erase(--L.end());
			R.insert(le - (b - r));
			a = *L.rbegin() + (a - l);
			b = le;
		}
		else {
			ans += abs(A[i] - b);
			R.insert(A[i] - (b - r));
			R.insert(A[i] - (b - r));
			ll rb = *R.begin() + (b - r);
			R.erase(R.begin());
			L.insert(rb - (a - l));
			a = rb;
			b = *R.begin() + (b - r);
		}
	}
	cout << ans << "\n";
}

JOI2012春day1-2: 魚

https://www.ioi-jp.org/camp/2012/2012-sp-tasks/

setを使う問題といえばコレ。直方体の体積を断面図を考えて求めていく問題です。

こういう問題でsetを使う際、
(1)番兵を使うと便利。
(2)lb,ubを覚えてまとめてerase(lb,ub)する。
(3)ある要素を見つけてそこからincrement or decrementする感じ。
かな?
ほとんどAlgoogleさんのコードと同じです。
JOI 春合宿 2012 Fish - Algoogle

struct cube {
	int p[3];
	cube(int x, int y, int z) { p[0] = x, p[1] = y, p[2] = z; }
	cube() { p[0] = p[1] = p[2] = 0;}
	friend ostream& operator<<(ostream& out, const cube& c){
		out << "(" << c.p[0] << ", " << c.p[1] << ", " << c.p[2] << ")"; return out; }
	bool operator<(const cube& c) { return p[2] > c.p[2]; }
};

int N;
char RGB[3] = {'R', 'G', 'B'};
pi P[MAX_N];
cube C[MAX_N];

void solve() {
	cin >> N;
	rep(i, 0, N) {
		int a; char c;
		cin >> a >> c;
		int b = 0;
		while(RGB[b] != c) b++;
		P[i] = pi(a, b);
	}
	sort(P, P + N);
	cube c(1, 1, 1);
	int r = 0;
	rep(i, 0, N) {
		while(r != N && P[r].fst - P[i].fst * 2 < 0) {
			c.p[P[r].sec]++;
			r++;
		}
		C[i] = c;
		c.p[P[i].sec]--;
	}
	C[N] = cube(0, 0, 0);
	sort(C, C + N + 1);

	ll ans = 0;
	ll area = 0, pz = inf;
	set<pi> s = {pi(0, inf), pi(inf, 0)};

	rep(i, 0, N + 1) {
		ll x = C[i].p[0], y = C[i].p[1], z = C[i].p[2];
		if(i != 0) ans += area * (pz - z);
		pz = z;
		set<pi>::iterator lb = s.lower_bound(pi(x, y)), ub = lb;
		int px = x, py = (*lb).sec;
		if(py >= y) continue;
		while(lb != s.begin()) {
			lb--;
			int nx = (*lb).fst, ny = (*lb).sec;
			area += 1LL * (px - nx) * (y - py);
			if(ny > y) break;
			px = nx; py = ny;
		}
		lb++;
		s.erase(lb, ub);
		s.insert(ub, pi(x, y));
	}
	cout << ans - 1 << "\n";
}

Codeforces Round #441E: Delivery Club

http://codeforces.com/contest/875/problem/E

ならし計算決めるアレです。
とりあえず距離kのとき行けるかどうかを判定する。
setでi番目まで処理し終わったとき、もう一人がどこにいるかを管理する。
i+1番目を考えて、もう一人の候補の座標としてありうるものを更新する。
一見O(N^2logN)だが、一度候補から外れると、もう候補となることはないのでならし計算量O(NlogN)となる。

結局全部値を持ってdpするだけですね。一見よくない戦略に見えて、すべての要素についてsetにぶち込まれる回数は定数倍回なのでならし計算量がうまくきまるという原理です。

よくよく考えてみるとstackやqueueを使う問題も同じです。単調性があってそのまま要素を末尾から足せばいいだけの場合はこれらを使ってO(N)とかにできます。

こういうパターンを僕は「愚直に値を持つdp」ということにします。

区間を扱う問題で単調性があったりするとこのパターンが利用できるかも?
dp高速化とかも愚直に全部値もってみて、初めてsegtreeなどデータ構造を使えることに気づくので発想は同じです。
ただしdp高速化は更新を高速にするだけで、データはすべて持ちます。「愚直に値をもつdp」はとりうる状態が限られていることを利用して計算量を減らします。

最初考察したとき二次元plotしたんですが、それは悪手でした。今回は2人を区別する必要がないからです。2つ値を持つからってxy-plotはさすがに頭が悪かった。使ってみて見通しよくなかったらすぐ切り替えましょう。

N=10^5だと貪欲だと思いがちなので、dpにも頭を働かせるようにしましょう。

int N, S1, S2;
int A[MAX_N];

bool C(int m) {
	if(abs(S1 - S2) > m) return false;
	set<int> s;
	s.insert(S1); s.insert(S2);
	rep(i, 0, N) {
		while(!s.empty() && abs(*s.begin() - A[i]) > m) s.erase(s.begin());
		while(!s.empty() && abs(*s.rbegin() - A[i]) > m) s.erase(*s.rbegin());
		if(i != 0 && abs(A[i - 1] - A[i]) <= m) {
			s.insert(A[i - 1]);
		}
		if(s.empty()) return false;
	}
	return true;
}

void solve() {
	cin >> N >> S1 >> S2;
	rep(i, 0, N) cin >> A[i];
	int lv = -1, rv = 1000000000;
	while(rv - lv > 1) {
		int m = (lv + rv) / 2;
		if(C(m)) rv = m;
		else lv = m;
	}
	cout << rv << "\n";
}

Tenka1 Programmer Contest: ModularPowerEquation!! (AC解法(解法なんだからACなのは当たり前だろ))

http://tenka1-2017.contest.atcoder.jp/tasks/tenka1_2017_f

http://omochan.hatenablog.com/entry/2017/10/20/004232の続きです。

普通に蟻本通りにやればいいだけでした…。
x≡k(modm)
x≡k'(modm')を満たすxは
x=mt+kとして
mt≡k'-k(modm')をtについて解けばいいだけでした。k,m,k',m'がすべて10^9より小さければ、xはオーバーフローせずに求められます。
具体的にはmt'≡g(modm')(g=gcd(m,m'))となるt'(modm')を求めてからt=(k'-k)/g*t'とすればいいだけです。

modの基本定理のおさらい。

  • g|Mのときga≡gb(modM)<=>a≡b(modM/g)
  • gとMが互いに素の時,ga≡gb(modM)<=>a≡b(modM)

二つとも両辺に寄せてからg(a-b)=Mtとかすれば導出できます。これが一番modの中で基礎となる式かもしれない。modの積と商の変形ルールを述べています。
ちなみに和(と差)は当たり前ですが、
a≡b(modM)<=>a+c≡b+c(modM)です。

四則演算の変形+中国剰余定理でだいぶmodは使いやすくなった?

ll phi(ll M) {
	ll A = M;
	for(ll i = 2; i * i <= A; i++) {
		if(A % i == 0) {
			M = M / i * (i - 1);
			while(A % i == 0) A /= i;
		}
	}
	if(A != 1) M = M / A * (A - 1);
	return M;
}

ll gcd(ll a, ll b) {
	return b == 0 ? a : gcd(b, a % b);
}

ll extgcd(ll a, ll b, ll& x, ll& y) { //get(ax + by = gcd(a, b))
	if(b == 0) {
		x = 1; y = 0;
		return a;
	}
	ll res = extgcd(b, a % b, y, x);
	y -= (a / b) * x;
	return res;
}

ll mod_pow(ll a, ll n, ll mod) { //get a^n mod
	if(n == 0) return 1;
	ll res = mod_pow(a, n / 2, mod);
	if(n % 2 == 0) return res * res % mod;
	else return a * res % mod * res % mod;
}

ll mod_inv(ll a, ll m) { //get x such that ax=gcd(a, m)(modm)
	ll x, y;
	extgcd(a, m, x, y);
	return (x % m + m) % m;
}

int Q;
ll A, M;

ll loop(ll A, ll M) {
	if(M == 1) return 10000;
	ll K = loop(A, phi(M));
	ll tK = mod_pow(A, K, M);
	return tK + (tK < 100 ? M : 0);
}

void solve() {
	cin >> Q;
	while(Q--) {
		cin >> A >> M;
		ll pM = phi(M);
		ll K = loop(A, M) % M;
		ll pK = loop(A, pM) % pM;
		ll pt = mod_inv(M, pM);
		ll g = gcd(M, pM), l = M * pM / g;
		ll t = ((pK - K) / g * pt) % pM;
		//debug(M, pM, K, pK, pt, g, l, t);
		ll pans = ((t * M % l + K) % l + l) % l;
		ll ans = pans < 100 ? pans + l : pans;
		cout << ans << "\n";
		//debug(mod_pow(A, ans, M), ans % M);
	}
}