mayoko’s diary

プロコンとかいろいろ。

AOJ 2638 Hyperrectangle

解法

まず, このページを見ると d 次元の立体における体積の求め方が分かります。各辺の長さが 1 の三角錐(?) の体積は 1/d! なんですね。
単体 (数学) - Wikipedia

で, 後は包除原理です。例えば 各辺の長さが {2, 3, 5}, S = 9 の場合を考えます。このとき,

  • 普通に全体を考えると 9^3 / 6 の体積を持つ
  • これだと余計に体積を考えているので体積を間引く。はみ出ているものとして,
  • これだと引きすぎているので体積を足す。はみ出ているものとして,

3 次元の図を上から見るとこんな感じです。もうちょっと足し匹するグループがありますがイメージは分かると思います。
f:id:mayokoex:20160811223744j:plain

普通に包除原理をすると O(2^d) かかりますが, 今回の包除原理は足し引き以外重要でないので, 「S から辺の長さをいくつ引いたか」の偶奇のみを覚えておけば良いです。つまり,

dp[d][rest][flag] = (d 番目の時点で, 残りの単体の辺の長さは rest で, S から辺を引いた回数の偶奇は flag であるときの体積) というような dp をやります。この考え方は過去の TC med にあったような気がします。
mayokoex.hatenablog.com

const int MAXD = 303;
const int MOD = 1e9+7;
int D, S;
int l[MAXD];
int dp[MAXD][MAXD*MAXD][2];

ll mod_pow(ll x, ll p, int MOD) {
	if (x == 0) return 0;
    ll a = 1;
    while (p) {
        if (p%2) a = a*x%MOD;
        x = x*x%MOD;
        p/=2;
    }
    return a;
}

int dfs(int d, int rest, int flag) {
	int& ret = dp[d][rest][flag];
	if (ret >= 0) return ret;
	if (d == D) {
		ret = mod_pow(rest, D, MOD);
		if (flag) ret = -ret;
		ret = (ret%MOD+MOD)%MOD;
	} else {
		ret = dfs(d+1, rest, flag);
		if (rest >= l[d]) {
			ret += dfs(d+1, rest-l[d], flag^1);
		}
		ret %= MOD;
	}

	return ret;
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
	cin >> D;
	for (int i = 0; i < D; i++)
		cin >> l[i];
	cin >> S;
	memset(dp, -1, sizeof(dp));
	cout << dfs(0, S, 0) << endl;
    return 0;
}

AOJ 2538 Stack Maze

問題

Stack Maze | Aizu Online Judge

H*W の迷路が与えられる。この迷路には小文字のアルファベットで表される宝石と, 大文字のアルファベットで表されるそれに対応する宝石を埋め込む穴がある(例えば 'a' とそれに対応する 'A' がある, という感じ)。

あなたは (0, 0) から (H-1, W-1) に (y, x) -> (y+1, x), (y, x) -> (y, x+1) という移動のみを許して進んでいきたい。また, この際になるべく多くの宝石を対応する穴に埋め込んでいきたい。最大でいくつの宝石を埋め込めるか。

解法

dp[y1][x1][y2][x2] = ((y1, x1) から (y2, x2) に移動するまでに最大でいくつの宝石を埋め込めるか) という区間 dp っぽいものを考えます。

遷移は

  • (y1, x1) を素通りする場合
    • (y1, x1) から下 or 右に移動してその最大値を答えとする
  • (y1, x1) が宝石の場合
    • それに対応する穴を y1 <= y <= y2, x1 <= x <= x2 の間で探す。で, dp[y1+dy][x1+dx][y-dy][x-dx] + dp[y][x][y2][x2] + 1 というようなものを考える。
    • ただ abs(y-y1) + abs(x-x1) == 1 の場合, 上のような遷移だと「(y1+dy, x1+dx) から (y-dy, x-dx) には行けないのでスルー」となってしまうので, この場合だけは例外的に扱います。

とすれば良いです。事前に can[y1][x1][y2][x2] = ( (y1, x1) から (y2, x2) に移動可能か) というのを事前計算しておけば上記の dp は O(A) (A は迷路の中にある同じ種類の穴の最大の数) でできます。問題文によると A <= 10 なのでこれで十分高速に答えを得られます。

const int INF = 1e9;
const int MAX = 55;

string board[MAX];
int dp[MAX][MAX][MAX][MAX];
bool can[MAX][MAX][MAX][MAX];
vector<pii> pnts[33];
int H, W;

int dfs(int y1, int x1, int y2, int x2) {
	int& ret = dp[y1][x1][y2][x2];
	if (ret != -1) return ret;
	if (y1 == y2 && x1 == x2) return ret = 0;
	ret = -INF;
	char c = board[y1][x1];
	if (c == '#') return ret;
	// 今いる点はスルーする
	if (y1 + 1 <= y2 && board[y1 + 1][x1] != '#') {
		ret = max(ret, dfs(y1 + 1, x1, y2, x2));
	}
	if (x1 + 1 <= x2 && board[y1][x1 + 1] != '#') {
		ret = max(ret, dfs(y1, x1 + 1, y2, x2));
	}
	if ('a' <= c && c <= 'z') {
		// 今いる点を選んでみる
		int index = c - 'a';
		for (auto p : pnts[index]) {
			int ny = p.first, nx = p.second;
			if (ny > y2 || nx > x2) continue;
			if (abs(ny - y1) + abs(nx - x1) == 1) {
				ret = max(ret, 1 + dfs(ny, nx, y2, x2));
			}
			else {
				for (int i = 0; i < 2; i++) for (int j = 0; j < 2; j++) {
					int ny1 = y1 + dy[i], nx1 = x1 + dx[i];
					if (ny1 > y2 || nx1 > x2) continue;
					int bny = ny - dy[j], bnx = nx - dx[j];
					if (bny < y1 || bnx < x1) continue;
					if (!can[ny1][nx1][bny][bnx]) continue;
					ret = max(ret, 1 + dfs(ny1, nx1, bny, bnx) + dfs(ny, nx, y2, x2));
				}
			}
		}
	}
	return ret;
}

int main() {
	cin.tie(0);
	ios::sync_with_stdio(false);
	while (cin >> H >> W) {
		if (H == 0 && W == 0) break;
		for (int i = 0; i < H; i++)
			cin >> board[i];
		for (int i = 0; i < 30; i++)
			pnts[i].clear();
		for (int i = 0; i < H; i++) for (int j = 0; j < W; j++) {
			char c = board[i][j];
			if ('A' <= c && c <= 'Z') {
				pnts[c - 'A'].emplace_back(i, j);
			}
		}
		memset(can, false, sizeof(can));
		for (int y = H - 1; y >= 0; y--) for (int x = W - 1; x >= 0; x--) {
			if (board[y][x] == '#') continue;
			can[y][x][y][x] = true;
			if (y + 1 < H && board[y + 1][x] != '#') {
				for (int y1 = y + 1; y1 < H; y1++) for (int x1 = x; x1 < W; x1++) {
					can[y][x][y1][x1] |= can[y + 1][x][y1][x1];
				}
			}
			if (x + 1 < W && board[y][x + 1] != '#') {
				for (int y1 = y; y1 < H; y1++) for (int x1 = x + 1; x1 < W; x1++) {
					can[y][x][y1][x1] |= can[y][x + 1][y1][x1];
				}
			}
		}
		memset(dp, -1, sizeof(dp));
		int ans = dfs(0, 0, H - 1, W - 1);
		if (ans < 0) ans = -1;
		cout << ans << endl;
	}
	return 0;
}

stack みたいな問題は順列っぽい -> 区間 dp?みたいな反射がある。

AOJ 2296 Quest of Merchant

Merchant, かわいい(小並感)

解法

すごく複雑そうな問題に見えますが, 落ち着いて問題を分解すると一個一個は大したことありません。

まず, 1 回の出張で訪れる町の集合がわかったら, その町の集合へ訪れる最短時間は巡回セールスマン問題を解けば簡単にわかります。また, 訪れる町の集合がわかると, 買うことのできるアイテム, およびそのアイテムを買うために必要な値段の最小値がわかります。ゆえに, そのアイテムの集合の中で得られる最大利益を計算できます。

これらを計算すると, dp[t] = (t の時間を消費して得ることのできる利益の最大値)という dp が普通にできます。普通に遷移を考えるとこの dp には O(T^2) かかるような気がしますが, 町への訪れ方は 2^N 通りしかないので, 各状態の遷移の仕方も 2^N 通りしかありません。よって十分時間内に答えが求められます。

int N, M, W, T;
const int MAX = 8;
const int MAXW = 11111;
string S[MAX];
int V[MAX], P[MAX];
int X[MAX], Y[MAX], L[MAX];
int R[MAX][MAX], Q[MAX][MAX];
ll dp[MAXW];

const int INF = 1e9;
int minTime[1 << MAX];
namespace TSP {
	int dp[1 << MAX][MAX];
	int dfs(int state, int now, int all) {
		int& ret = dp[state][now];
		if (ret >= 0) return ret;
		ret = INF;
		if (state == all) return ret = abs(X[now]) + abs(Y[now]);
		for (int i = 0; i < N; i++) if ((all >> i) & 1) {
			if ((state >> i) & 1) continue;
			ret = min(ret, abs(X[now] - X[i]) + abs(Y[now] - Y[i]) + dfs(state | (1 << i), i, all));
		}
		return ret;
	}
}

int best[1 << MAX];
namespace Item {
	ll dp[MAXW];
	ll calc(int state, int M, int W, vi costs) {
		memset(dp, 0, sizeof(dp));
		for (int i = 0; i < M; i++) {
			if ((state>>i)&1) {
				for (int j = 0; j+V[i] <= W; j++) {
					dp[j+V[i]] = max(dp[j+V[i]], dp[j] + P[i]-costs[i]);
				}
			}
		}
		ll ans = 0;
		for (int i = 0; i <= W; i++)
			ans = max(ans, dp[i]);
		return ans;
	}
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
	cin >> N >> M >> W >> T;
	for (int i = 0; i < M; i++) {
		cin >> S[i] >> V[i] >> P[i];
	}
	for (int i = 0; i < N; i++) {
		cin >> L[i] >> X[i] >> Y[i];
		for (int j = 0; j < L[i]; j++) {
			string s;
			cin >> s >> Q[i][j];
			for (int k = 0; k < M; k++) {
				if (S[k] == s) R[i][j] = k;
			}
		}
	}
	// ある集合の町を訪れる際に必要な時間の最小値をすべて求める
	for (int s = 1; s < 1 << N; s++) {
		memset(TSP::dp, -1, sizeof(TSP::dp));
		int mini = INF;
		for (int i = 0; i < N; i++) if ((s >> i) & 1) {
			mini = min(mini, abs(X[i]) + abs(Y[i]) + TSP::dfs(1 << i, i, s));
		}
		minTime[s] = mini;
	}
	map<int, ll> mp;
	for (int s = 1; s < 1<<N; s++) {
		if (minTime[s] > T) continue;
		// この町の訪れ方で得られるアイテムの集合
		int item = 0;
		vector<int> unko(M, INF);
		for (int j = 0; j < N; j++) if ((s>>j)&1) {
			for (int k = 0; k < L[j]; k++) {
				int l = R[j][k];
				item |= 1<<l;
				unko[l] = min(unko[l], Q[j][k]);
			}
		}
		ll tmp = Item::calc(item, M, W, unko);
		//cout << s << " " << tmp << endl;
		mp[minTime[s]] = max(mp[minTime[s]], tmp);
	}
	// dp
	for (int i = 0; i <= T; i++) {
		for (auto p : mp) {
			if (i+p.first <= T) {
				dp[i+p.first] = max(dp[i+p.first], dp[i] + p.second);
			}
		}
	}
	ll ans = 0;
	for (int i = 0; i <= T; i++)
		ans = max(ans, dp[i]);
	cout << ans << endl;
    return 0;
}

AOJ 2168 Luigi's Tavern

問題

Luigi's Tavern | Aizu Online Judge

H 人の勇者, W 人の戦士, C 人の僧侶, M 人の魔法使いがいる。以下の条件を満たすパーティは最大いくつ作れるか。

  • パーティには必ず勇者がいなければならない。
  • 同じパーティの勇者と戦士は仲良くなければならない。
  • 同じパーティの戦士と僧侶は仲良くなければならない。
  • 同じパーティの僧侶と魔法使いは仲良くなければならない。
  • パーティには必ずしも 4 種のジョブがそろってなくても良い。ただし, 僧侶がいない場合は戦士, 魔法使いはいなければならない。
  • パーティには戦士, 僧侶, 魔法使いがいなくても良いが, それぞれ Nw, Nc, Nm 個以下のパーティでそうでなければならない。
解法

見た目からしてマッチング -> 最大フローという流れが目に浮かびます。具体的には以下のようなグラフを作れば良いです。

f:id:mayokoex:20160807211106j:plain

汚くてすみません。
個人的に引っかかったのですが, Nw とかの辺を意識するのはすぐわかったんですが普通の辺に対してもバッファみたいのを付けないとその頂点からフローが 2 以上流れる可能性があるのを忘れてました。

あと頂点のラベル付けがめっちゃめんどくさいのである程度まとまりをつけて変数化すべきですね。

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
	int H, W, C, M, Nw, Nc, Nm;
	while (cin >> H >> W >> C >> M >> Nw >> Nc >> Nm) {
		if (H == -1) break;
		int s = H + 2 * (W + 1) + 2 * (C + 1) + 2 * (M + 1);
		int t = s + 1, V = t + 1;
		Dinic dinic(V);
		int w1 = H + W, w2 = w1 + W + 1;
		int c1 = w2 + C + 1, c2 = c1 + C + 1;
		int m1 = c2 + M + 1, m2 = m1 + M + 1;

		for (int i = 0; i < H; i++) {
			dinic.add_edge(s, i, 1);
			dinic.add_edge(i, w1, 1);
		}

		for (int i = 0; i < W; i++) {
			// 上 -> 下
			dinic.add_edge(H + i, w1 + 1 + i, 1);
			// 次のプールにつなげる
			dinic.add_edge(w1 + 1 + i, c1, 1);
		}
		dinic.add_edge(w1, w2, Nw);
		
		for (int i = 0; i < C; i++) {
			// 上 -> 下
			dinic.add_edge(w2 + 1 + i, c1 + 1 + i, 1);
			// 次のプールにつなげる
			dinic.add_edge(c1 + 1 + i, m1, 1);
			// 前のプールからつなげる
			dinic.add_edge(w2, w2 + 1 + i, 1);
		}
		dinic.add_edge(c1, c2, Nc);

		for (int i = 0; i < M; i++) {
			// 上 -> 下
			dinic.add_edge(c2 + 1 + i, m1 + 1 + i, 1);
			// 前のプールからつなげる
			dinic.add_edge(c2, c2 + 1 + i, 1);
			// t につなげる
			dinic.add_edge(m1 + 1 + i, t, 1);
		}
		dinic.add_edge(m1, m2, Nm);
		dinic.add_edge(m2, t, Nm);

		for (int i = 0; i < W; i++) {
			int N;
			cin >> N;
			for (int j = 0; j < N; j++) {
				int h;
				cin >> h; h--;
				dinic.add_edge(h, H + i, 1);
			}
		}
		for (int i = 0; i < C; i++) {
			int N;
			cin >> N;
			for (int j = 0; j < N; j++) {
				int w;
				cin >> w; w--;
				dinic.add_edge(w + H, w2 + 1 + i, 1);
			}
		}
		for (int i = 0; i < M; i++) {
			int N;
			cin >> N;
			for (int j = 0; j < N; j++) {
				int c;
				cin >> c; c--;
				dinic.add_edge(c1 + 1 + c, m1 + 1 + i, 1);
			}
		}
		cout << dinic.max_flow(s, t) << endl;
	}
    return 0;
}

yukicoder No.409 ダイエット

この問題実は btk さんから相談をもらっていたんですけど時間内に解けなかった…(主に蟻本を探していた)

解法

dp[i] = (i 日目に最後にドーナツを食べる際の最小体重)とします。この発想の時点で結構頭良い気がします(こういう風に書くと状態に「何日間ドーナツを食べていないか」というのを持たなくても良いところが good)。

これで普通に状態遷移を考えると,

dp[i] = min(dp[j] - (i-j-1)*A + (i-j)*(i-j-1)/2 * B + D[i])

となり, これをすべての j について調べることになります。これを普通にやると O(n^2) かかるんですが, 上式を変形して

dp[i] = D[i] - (i-1)*A + (i*i-i)/2*B + min(-B*j*i + j*A + (j*j+j)/2*B + dp[j])

というように変換すると, 結局各 j について考えるべきは, min の中身だけになります。さらに, min の中身は i の一次式の形をしているので, 蟻本に書いてある convex-hull trick を用いることで全体として O(n) の時間で各 i の min 値を求めることができます。

const int MAXN = 300300;
ll dp[MAXN];
int deq[MAXN];
int N, A, B, W;

bool check(int f1, int f2, int f3) {
    ll a1 = -(ll)B*f1, b1 = (ll)f1*A + ((ll)f1*f1+f1)/2 * B + dp[f1];
    ll a2 = -(ll)B*f2, b2 = (ll)f2*A + ((ll)f2*f2+f2)/2 * B + dp[f2];
    ll a3 = -(ll)B*f3, b3 = (ll)f3*A + ((ll)f3*f3+f3)/2 * B + dp[f3];
    return (a2-a1) * (b3-b2) >= (b2-b1) * (a3-a2);
}

ll f(int j, int x) {
    return -(ll)B*j*x + (ll)j*A + ((ll)j*j+j)/2*B + dp[j];
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    cin >> N >> A >> B >> W;
    vector<int> D(N);
    for (int i = 0; i < N; i++)
        cin >> D[i];
    dp[0] = 0;
    int s = 0, t = 1;
    deq[0] = 0;
    for (int i = 1; i <= N; i++) {
        while (s+1 < t && f(deq[s], i) >= f(deq[s+1], i)) s++;
        dp[i] = f(deq[s], i) + D[i-1] - (ll)(i-1)*A + ((ll)i*i-i)/2*B;
        // 末尾から最小値を取る可能性がなくなったものを取り除く
        while (s+1 < t && check(deq[t-2], deq[t-1], i)) t--;
        deq[t++] = i;
    }
    ll ans = 1ll<<60;
    for (int i = 0; i <= N; i++) {
        ll d = N-i;
        ans = min(ans, dp[i]+B*(d+1)*d/2 - A*d);
    }
    cout << ans+W << endl;
    return 0;
}

yukicoder No.408 五輪ピック

解法

「頂点 1 からの距離が 2 の頂点 a, b が辺でつながっていたら, 1 -> c -> a -> b -> d -> 1 みたいな感じでループにならない?」っていうのが基本的な発想です。

つまり, 各頂点について「1 からの距離が 2 になりうるか」というのをまず覚えておきます。

で, そのあと各辺を調べて, その辺を構成する頂点 a, b が 頂点 1 からの距離 2 になることができるかを判定し, 両方とも距離が 2 になりうるなら YES, とする方針です。

ただ, これだと例えば 1-c, c-a, 1-a, a-b と辺がある場合とかで勘違いをします(a, b に辺があって, 1-a と 1-b 間の距離は 2 になりうるけどこれは長さ 5 のループになっていない)。

ここをうまくやるために, 距離が 2 になる各頂点 a について, 「1 - c - a と渡って 1-a 間の距離が 2 になった場合, その中間点 c」を覚えておくことにします(この c は 3 つ覚えておけば十分)。
この c がうまいこと条件を満たしていれば上と同じような考えで長さ 5 のループになっていることがわかります。

const int MAX = 20020;
vector<int> from[MAX];

void check(int v, const vector<vi>& G) {
    for (int ch : G[v]) if (ch != 0) {
        if (from[ch].size() < 3) from[ch].push_back(v);
    }
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    int N, M;
    cin >> N >> M;
    vector<vector<int> > G(N);
    for (int i = 0; i < M; i++) {
        int a, b;
        cin >> a >> b;
        a--; b--;
        G[a].emplace_back(b);
        G[b].emplace_back(a);
    }
    for (int v : G[0]) {
        check(v, G);
    }
    for (int i = 1; i < N; i++) {
        for (int ch : G[i]) {
            for (int j = 0; j < from[i].size(); j++) for (int k = 0; k < from[ch].size(); k++) {
                if (from[i][j] != ch && from[i][j] != from[ch][k] && from[ch][k] != i) {
                    cout << "YES" << endl;
                    return 0;
                }
            }
        }
    }
    cout << "NO" << endl;
    return 0;
}

AOJ 2305 Beautiful Currency

問題

Beautiful Currency | Aizu Online Judge

素数 n の数列 a が beautiful であるとは, 1 以上 n-1 以下の任意の i について, a[i+1] % a[i] == 0 となることである。

数列 a のいくつかの数を入れ替えて, a を beautiful にしたい。a[i] の値を b[i] に変更した際のコストは abs(a[i]-b[i])/a[i] で与えられる。数列 a を b に変更した際のコストは上記コストの最大値で与えられる。

数列 a を beautiful にするための最小コストを求めよ。

解法

まず考察ですが, 全ての数を 1 にすれば beautiful になるので, 答えは絶対 1 以下です。つまり, 各 a[i] について, その値は 1 ~ 2*a[i] の範囲のみを考えれば良いですね。

また, 2*10^5 以下の数では約数の数はそれほど多くならなそう(10^5 までだと 128 個が最高)です。とすると, [1, 200000] の範囲で各 i の約数の数の合計は まぁそれほど多くなさそうです。

そこで, 数列 a を反転させて, dp[i][j] = (i 番目を数 j にしたときの最小コスト) というものを考えます。各 j に対して, i+1 番目の要素は a[i] の約数になっていないといけないので, i+1 番目の数の候補は j の約数のみです。

あらかじめ [1, 200000] の約数を求めておけば, これは十分早く動作します。

const int MAXN = 22;
const int MAX = 202000;
double dp[MAXN][MAX];
vi dd[MAX];

vi divs(int n) {
	vi ret;
	for (int i = 1; i*i <= n; i++) {
		if (n%i == 0) {
			ret.push_back(i);
			if (i*i < n) ret.push_back(n / i);
		}
	}
	return ret;
}

int main() {
	for (int i = 1; i < MAX; i++)
		dd[i] = divs(i);
	int N;
	cin >> N;
	vector<int> A(N);
	for (int i = 0; i < N; i++)
		cin >> A[i];
	reverse(A.begin(), A.end());
	for (int i = 0; i < N; i++) for (int j = 0; j < MAX; j++)
		dp[i][j] = 1;
	for (int i = 1; i < 2 * A[0]; i++) {
		dp[0][i] = 1.*abs(i - A[0]) / A[0];
	}
	for (int i = 0; i < N-1; i++) {
		for (int j = 1; j < 2 * A[i]; j++) if (dp[i][j] < 1) {
			for (int d : dd[j]) {
				double tmp = max(dp[i][j], 1.*abs(A[i + 1] - d) / A[i + 1]);
				dp[i + 1][d] = min(dp[i + 1][d], tmp);
			}
		}
	}
	double ans = 1;
	for (int i = 0; i < MAX; i++)
		ans = min(ans, dp[N - 1][i]);
	printf("%.10lf\n", ans);
    return 0;
}