ARC058 E - 和風いろはちゃん
解法
X, Y, Z が小さいので bitDP でやります．
dp[i][S] := i 番目まで見たときに，直前の和が X + Y + Z 以下になるところまでの状態が S であるような場合の数
とします．
たとえば，X = 5, Y = 7, Z = 5 のときに，"575" という状態を "10000 1000000 10000" という (X + Y + Z) ビットからなる列として表します．
XYZ であるか，という状態は，よく考えると下から Z ビット目，(Z + Y) ビット目，(X + Y + Z) ビット目が 1 になっているか，という問と同じなので，そういうマスクを作っておけば判定は簡単です．
過去に XYZ を含んだ数列は次からは任意の数を取って良いので，その数列を特別な場所に保存しておくと単に 10 倍するだけで済みます．
ソースコード
#include <bits/stdc++.h> using namespace std; using ll = long long; constexpr int mod = 1e9 + 7; int main() { int N, X, Y, Z; cin >> N >> X >> Y >> Z; int const M = 1 << (X + Y + Z); vector<vector<ll>> dp(N + 1, vector<ll>(M + 1)); int const mask = (1 << (Z - 1)) | (1 << (Z + Y - 1)) | (1 << (Z + Y + X - 1)); dp[0][0] = 1; for(int i = 0; i < N; ++i) { for(int S = 0; S < M; ++S) { if(dp[i][S] == 0) { continue; } for(int j = 1; j <= 10; ++j) { int nS = ((S << j) | (1 << (j - 1))) & (M - 1); if((nS & mask) == mask) { (dp[i + 1][M] += dp[i][S]) %= mod; } else { (dp[i + 1][nS] += dp[i][S]) %= mod; } } } (dp[i + 1][M] += dp[i][M] * 10) %= mod; } cout << dp[N].back() << endl; }
感想
bitDP だとわかれば簡単（ほんまか？）
AOJ 1301 Malfatti Circles
解法
まず，一つの円を決め打ちします．
すると，残りの2円の位置は，二分探索などをして求めることが出来ます．
ここで，残りの2円が離れすぎている場合，最初に決め打ちした円が小さすぎるということになります．逆に，近すぎる場合（交わっている場合），最初に決め打ちした円が大きすぎるということになります．
なので，最初の円を二分探索でちょうどいいサイズにすると，答えが求まります．
二分探索の上限の円の中心位置は，三角形の重心にしておかないとめんどくさいことになりそうな気がします．
ソースコード
#include <bits/stdc++.h> using namespace std; using ld = long double; constexpr ld eps = 1e-8; // ライブラリ部分は長いので略． point rot(point p, ld ang) { ld x = real(p) * cos(ang) - imag(p) * sin(ang); ld y = real(p) * sin(ang) + imag(p) * cos(ang); return point(x, y); } circle calc(point cp, ld r, point const& limit, line l1, line l2) { point lb = l1.a, ub = limit; ld rad = 0; for(int i = 0; i < 60; ++i) { point m = 0.5l * (lb + ub); rad = dist_lp(l1, m); if(abs(cp - m) + eps < r + rad) { ub = m; } else { lb = m; } } return circle(lb, rad); } int main() { int x1, y1, x2, y2, x3, y3; while(cin >> x1 >> y1 >> x2 >> y2 >> x3 >> y3) { if(x1 == 0 && y1 == 0 && x2 == 0 && y2 == 0 && x3 == 0 && y3 == 0) { break; } point p1(x1, y1), p2(x2, y2), p3(x3, y3); point vec1 = rot(p2 - p1, (arg(p3 - p1) + arg(p2 - p1)) / 2 - arg(p2 - p1)); point vec2 = rot(p1 - p3, (arg(p1 - p3) + arg(p2 - p3)) / 2 - arg(p1 - p3)); point g = is_ll(line(p1, p1 + vec1), line(p3, p3 + vec2)); point lb = p1, ub = g; ld r1, r2, r3; for(int i = 0; i < 60; ++i) { point m = (ub + lb) * 0.5l; r1 = dist_lp(line(p1, p2), m); auto c2 = calc(m, r1, g, line(p2, p1), line(p2, p3)); auto c3 = calc(m, r1, g, line(p3, p1), line(p3, p2)); r2 = c2.r; r3 = c3.r; if(abs(c2.p - c3.p) + eps < c2.r + c3.r) { lb = m; } else { ub = m; } } cout << setprecision(10) << fixed; cout << r1 << ' ' << r2 << ' ' << r3 << endl; } }
感想
多分かなり遅い実装をしていて，ループ回数をそれぞれ100ずつにするとTLEギリギリ通るぐらいでした（座標のサイズ的に100回も回す必要はないんですけどね…）．
まあ通ればなんでもいいよね．速度的には一方の二分探索をやめて，2次方程式で直接求めるほうが早いでしょう．
ARC067 E - Grouping
解法
dp[i + 1][j] := i 人以下のグループのみを用いて，あと j 人グループ分けされてない人がいるときの場合の数
という dp で解けます．
遷移は，use == 0 または C <= use <= D のとき
$$\displaystyle dp[i + 1][j - use \times i] = dp[i][j] \times \frac{_jC_i \times _{j-i}C_i \times \ldots \times _{j-(use-1)\times i}C_i}{(use)!}$$
になります．これは，j 人から i 人を選ぶ，という操作を use 回繰り返すが，グループ自体の区別はつけないので，グループの順列で割るということです．
ソースコード
#include <bits/stdc++.h> using namespace std; // mod は略 int main() { int N, A, B, C, D; cin >> N >> A >> B >> C >> D; vector<vector<mod>> dp(B + 2, vector<mod>(N + 1)); dp[A][N] = 1; for(int i = A; i <= B; ++i) { for(int j = 0; j <= N; ++j) { if(dp[i][j] == 0) { continue; } mod t = dp[i][j]; for(int use = 0; use * i <= j; ++use) { if(use == 0 || C <= use && use <= D) { dp[i + 1][j - use * i] += t * fact<MOD>(use, false); } t *= comb<MOD>(j - use * i, i); } } } cout << (int)dp[B + 1][0] << endl; }
感想
数え上げる問題が苦手すぎる．中学高校で訓練しておけばよかった…．
AOJ 2309 Vector Compression
問題文
まず問題読解が難しかった．
http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2309
問題概要
N次元のベクトルが M 個与えられる．これを v[1], ..., v[M] とする．
M 個のベクトルを好きな順番に並べ替え，その順に記録していく．記録する際にコストがかかる．コストは，|v[i]|^2 か，それ以前に記録していたベクトルの中から好きなもの v[j] を選んで，|v[i] - r * v[j]|^2 のどちらかを選択できる．ただし，r は記録のたびに好きな値を適当にとって良いとする．
この時，すべてのベクトルを記録するのにかかるコストを最小化せよ．
・制約
1 <= N, M <= 100
解法
v[i] を記録するときに v[j] を使う，というのを有向辺 (j, i) （向きに注意）であらわし，そのコストを |v[i] - r * v[j]|^2 の最小値とする．
また，v[i] を直接使う場合は，有向辺 (M, i) で表し，そのコストを |v[i]|^2 とする．
すると，この問題はこうして出来たグラフの最小有向全域木を求めよ，と言い換えることができる．
r の求め方だが，幸い |v[i] - r * v[j]|^2 は明らかに下に凸であり，かつ r で微分するのも簡単なので，極値が簡単に求められる（平方完成してもいい）．
計算量は O(nm^2)．
ソースコード
#include <bits/stdc++.h> using namespace std; using weight = double; constexpr weight INF = 1e9; constexpr weight eps = 1e-8; struct edge { int from, to; weight cost; bool operator<(edge const& e) const { return cost < e.cost; } }; using edges = std::vector<edge>; using graph = std::vector<edges>; void add_edge(graph& g, int from, int to, weight cost) { g[from].push_back(edge{from, to, cost}); } int scc(std::vector<std::vector<int>>& G, std::vector<int>& cmp) { int V = G.size(); std::vector<std::vector<int>> g(V), rg(V); std::vector<bool> used(V, false); std::vector<int> vs; cmp.resize(V); for(int i = 0; i < V; ++i) { for(auto to : G[i]) { g[i].push_back(to); rg[to].push_back(i); } } std::function<void(int)> dfs = [&g, &vs, &used, &dfs](int v) { used[v] = true; for(auto i : g[v]) { if(!used[i]) { dfs(i); } } vs.push_back(v); }; std::function<void(int, int)> rdfs = [&rg, &cmp, &used, &rdfs](int v, int k) { used[v] = true; cmp[v] = k; for(int i : rg[v]) { if(!used[i]) { rdfs(i, k); } } }; for(int v=0; v<V; ++v) { if(!used[v]) { dfs(v); } } std::fill(used.begin(), used.end(), false); int k = 0; for(int i=vs.size()-1; i>=0; --i) { if(!used[vs[i]]) { rdfs(vs[i], k++); } } return k; } std::vector<std::vector<int>> build_graph(std::vector<std::vector<int>> const& g, std::vector<int> const& cmp, int K) { int V = g.size(); std::vector<std::set<int>> s(K); std::vector<std::vector<int>> res(K); for(int i = 0; i < V; ++i) { for(auto to : g[i]) { s[cmp[i]].insert(cmp[to]); } } for(int i = 0; i < K; ++i) { for(auto j : s[i]) { if(i != j) { res[i].push_back(j); } } } return res; } weight minimum_spanning_rooted_arborescence(graph const& g, int root, weight sum = 0) { int const n = g.size(); std::vector<int> rev(n, -1); std::vector<weight> cost(n, INF); for(int i = 0; i < n; ++i) { for(auto& e : g[i]) { if(e.cost < cost[e.to]) { cost[e.to] = e.cost; rev[e.to] = i; } } } for(int i = 0; i < n; ++i) { if(i != root && rev[i] == -1) { // not exists return INF; } } std::vector<std::vector<int>> g2(n); for(int i = 0; i < n; ++i) { if(root == i) { continue; } g2[rev[i]].push_back(i); sum += cost[i]; } std::vector<int> cmp(n); int const K = scc(g2, cmp); auto nxt = build_graph(g2, cmp, K); if(nxt.size() == n) { return sum; } graph ng(nxt.size()); for(int i = 0; i < n; ++i) { for(auto& e : g[i]) { if(cmp[i] == cmp[e.to]) { continue; } ng[cmp[i]].push_back(edge{cmp[i], cmp[e.to], e.cost - cost[e.to]}); } } return minimum_spanning_rooted_arborescence(ng, cmp[root], sum); } int main() { int N, M; cin >> N >> M; vector<vector<double>> v; for(int i = 0; i < M; ++i) { double t = 0; vector<weight> w(N); for(int j = 0; j < N; ++j) { cin >> w[j]; t += w[j] * w[j]; } if(t > eps) { v.push_back(w); } } M = v.size(); graph g(M + 1); int const root = M; for(int i = 0; i < M; ++i) { for(int j = 0; j < M; ++j) { if(i == j) { continue; } double prod = 0, dom = 0; for(int k = 0; k < N; ++k) { prod += v[i][k] * v[j][k]; dom += v[j][k] * v[j][k]; } double r = prod / dom; double cost = 0; for(int k = 0; k < N; ++k) { cost += (v[i][k] - r * v[j][k]) * (v[i][k] - r * v[j][k]); } add_edge(g, j, i, cost); } double d = 0; for(int j = 0; j < N; ++j) { d += v[i][j] * v[i][j]; } add_edge(g, root, i, d); } cout << fixed << setprecision(10) << minimum_spanning_rooted_arborescence(g, root) << endl; }
感想
最初フローかなあと思っていて全然解けず困っていた．
最小有向全域木を書いたのははじめてなのでいい勉強になりました．
AOJ 2446 Enumeration
解法
簡単のために，試しに2つの数字 a, b を用意して考えてみます．a, b が選ばれる確率をそれぞれ p, q とします．
まず，a のみを選んだ場合の期待値は，
- (m / a) * p * (1 - q)
となります．
b のみを選んだ場合は，
- (m / b) * (1 - p) * q
となります．
最後に a と b の両方を選んだ場合は，
- (m / a + m / b - m / lcm(a, b)) * p * q
となります．
これらを足し合わせてみると，以下のようなきれいな式になります．
- (m / a) * p + (m / b) * q - (m / lcm(a, b)) * p * q
これは一般の場合でも成り立つので，結局包除原理で簡単に求めることが出来ます． 2^n 通り試してやればよいです．
lcm のオーバーフローに注意．
ソースコード
#include <bits/stdc++.h> using namespace std; using ll = long long; constexpr ll INF = 2e18; constexpr ll invalid = -1; ll lcm(ll a, ll b) { ll g = __gcd(a, b); a /= g; if(a > INF / b) { return invalid; } return a * b; } int main() { int n; ll m; cin >> n >> m; vector<ll> a(n); vector<double> p(n); for(int i = 0; i < n; ++i) { cin >> a[i]; } for(int i = 0; i < n; ++i) { cin >> p[i]; p[i] /= 100; } double res = 0; for(int S = 1; S < 1 << n; ++S) { int cnt = 0; double q = 1; for(int i = 0; i < n; ++i) { if((S >> i) & 1) { cnt++; q *= p[i]; } } ll t = 1; for(int i = 0; i < n; ++i) { if((S >> i) & 1) { t = lcm(t, a[i]); } if(t == invalid) { break; } } if(t != invalid) { if(cnt & 1) { res += q * (m / t); } else { res -= q * (m / t); } } } cout << fixed << setprecision(10) << res << endl; }
Codeforces Round #230 (Div.1) B. Tower of Hanoi
問題概要
ハノイの塔のパズルの最小手数は 2^n - 1 であることは有名な事実であるが，今回は以下の問題を考えるとする．
ある段を位置 i から位置 j に移すのにコストが t[i][j] かかるとするとき( 1 <= i, j <= 3 )，最初位置 0 にある n 段の塔を，位置 2 に移すときのトータルのコストを最小化せよ．手数は最小化しなくともよい．
・制約
1 <= n <= 40
解法
最小手数が 2^n - 1 である，ということの漸化式の証明が頭にあれば，すぐに dp の式を導けると思います．
dp[i][from][to] := i 段の塔を位置 from から to に移すときの最小コスト
とします．また，もう一つの位置を other とします．
dp[i][from][to] は，以下の2つのどちらかです．
- i - 1 段の塔を from から other に移す．その後，i 段目を from から to に移す．最後に i - 1 段の塔を other から to に移す．
- i - 1 段の塔を from から to に移す．次に i 段目を other に移す．その後 i - 1 段の塔をまた from に戻す．そして i 段目を other から to に移す．最後に i - 1 段の塔を from から to に移す．
ソースコード
#include <bits/stdc++.h> using namespace std; using ll = long long; int main() { int t[3][3] = {}; int n; for(int i = 0; i < 3; ++i) { for(int j = 0; j < 3; ++j) { cin >> t[i][j]; } } cin >> n; ll dp[41][3][3] = {}; for(int i = 1; i <= n; ++i) { for(int from = 0; from < 3; ++from) { for(int to = 0; to < 3; ++to) { int other = 3 - from - to; dp[i][from][to] = t[from][to] + dp[i - 1][from][other] + dp[i - 1][other][to]; ll tmp = 2 * dp[i - 1][from][to] + dp[i - 1][to][from] + t[from][other] + t[other][to]; dp[i][from][to] = min(dp[i][from][to], tmp); } } } cout << dp[n][0][2] << endl; }
AOJ 2611 Ordering
解法
問題文を丁寧に読まなければならない問題です．
「それぞれの塔は自分自身よりも大きな塔とは高々1回しか比較されていない」とあるので，与えられた関係 (a, b) を有向辺とみるとこれは木になっています．
この時点で木DPかなあという見通しが立ちます．
ある頂点のそれぞれの部分木は独立に考えることができるので，あとはそれらをいい感じにマージしていくことになります．
引き分けがなければ単純に挿入していくDPなんですが，今回は引き分けがあるのでうまく処理しなければなりません．
（高さの異なる塔を高さ順に並べていくイメージだとすると，引き分けは同じところに複数の塔を置くということになります．）
ここで以下のDPを考えます．
dp2[i][j][k] := 高さが異なる i 個の集合 A と，高さが異なる j 個の集合 B をマージして，最終的に高さが異なる k 個の集合にする場合の数
遷移は以下の通り，3通りあります．
- A の i 番目の塔が B の j 番目よりも高い => dp2[i - 1][j][k - 1]
- A の i 番目と B の j 番目の高さが等しい => dp2[i - 1][j - 1][k - 1]
- A の i 番目より B の j 番目のほうが高い => dp2[i][j - 1][k - 1]
dp2[i][j][k] はこれらの和になります．
このDPテーブルは前もって計算できるのでそうします．
あとはマージしていくだけで，これは難しくありません．以下のDPでできます．
dp[v][i] := 根が v の部分木で，高さが異なる塔が i 個あるような関係の作り方の総数
この dp テーブルを先程の dp2 を使って計算していきます．
今構築中の高さの関係を
tmp[i][j] := i 番目の子まで見たときに，高さが異なる塔が j 個であるような総数
とすると，
tmp[i + 1][l] += tmp[i][j] * dp[ch[i]][k] * dp2[j][k][l]
という式が書けるので，各 i, j, k, l についてループを回します．
各頂点でこれをやるので，ぱっと見た感じでは（ch の i は無視できるので）j, k, l の3重ループと合わせて O(N^4) っぽいですが，計算量の漸化式を書いてみると実は O(N^3) になっています．
ソースコード
#include <bits/stdc++.h> using namespace std; using ll = long long; constexpr ll M = 1e9 + 7; ll dp2[256][256][256]; int sz[256]; int dfs(int v, vector<vector<int>>& g, vector<vector<ll>>& dp) { if(g[v].size() == 0) { dp[v][1] = 1; return sz[v] = 1; } sz[v] = 1; for(auto ch : g[v]) { sz[v] += dfs(ch, g, dp); } vector<vector<ll>> tmp(g[v].size() + 1, vector<ll>(sum)); tmp[0][0] = 1; for(int i = 0; i < (int)g[v].size(); ++i) { int ch = g[v][i]; for(int l = 0; l < sz[v]; ++l) { for(int j = 0; j <= l; ++j) { for(int k = 1; k <= sz[ch]; ++k) { (tmp[i + 1][l] += (tmp[i][j] * dp[ch][k] % M) * dp2[j][k][l]) %= M; } } } } for(int i = 0; i < sum; ++i) { dp[v][i + 1] = tmp[g[v].size()][i]; } return sz[v]; } int main() { for(int i = 0; i < 201; ++i) { dp2[i][0][i] = 1; dp2[0][i][i] = 1; } for(int k = 1; k < 201; ++k) { for(int i = 1; i < 201; ++i) { for(int j = 1; j < 201; ++j) { (dp2[i][j][k] += dp2[i - 1][j - 1][k - 1] + dp2[i - 1][j][k - 1] + dp2[i][j - 1][k - 1]) %= M; } } } int N; cin >> N; vector<vector<int>> g(N); vector<bool> in(N); for(int i = 0; i < N - 1; ++i) { int a, b; cin >> a >> b; g[a].push_back(b); in[b] = true; } int root = -1; for(int i = 0; i < N; ++i) { if(!in[i]) { root = i; break; } } vector<vector<ll>> dp(N + 1, vector<ll>(N + 1)); dfs(root, g, dp); ll res = 0; for(int i = 0; i <= N; ++i) { (res += dp[root][i]) %= M; } cout << res << endl; }
感想
O(N^4) から落ちへんやんけ！と言いながらヤケクソで投げたら通り，計算量の見積もりが甘かったと気付かされた．
実はオーダーが落ちて（？）いるというのは，結構面白かったし今後も使えそうなので覚えておきたい．
