HackerRank Boxes for Toys
問題概要
N個の箱が一列に並んでいて, それぞれ大きさは a_i * b_i * c_i . これらは自由に回転させられる.
これらのうちランダムな区間を選び, 含まれる箱のどれでも包含することができるような箱の体積の最小値の期待値 mod 1e9+7 を求めよ.
解法
期待値なので任意の区間の値の総和を求めて N*(N+1)/2 で割る.
まず, 区間に含まれる箱をどれでも包含できる箱の体積の最小値を求める方法を考える.
各箱は自由に回転させられるが, この３つの辺を１番長いもの、２番目に長い辺、３番目に長い辺でそれぞれ独立に長さの最大値を求めてその積が最小になる.
これはすぐ証明できる.
しかし, これでもまだ O(N^2) かかる.
そこで, 分割統治で考えてみる.
区間が長さ１のとき, その体積を返す.
区間が長さ２以上のとき, 真ん中の２つを含むような任意の区間を数える.
左を一つずつ伸ばしながら固定すると, 右端は４つに分類できる.
- 各辺の左側の最大値が右側の最大値より３つとも大きいもの.
- 各辺の左側の最大値が右側の最大値より２つ大きいもの.
- 各辺の左側の最大値が右側の最大値より１つ大きいもの.
- 各辺の左側の最大値が右側の最大値より全て小さいもの.
これらはそれぞれ連続する区間になっていて, 各区間の右端での最小の体積の和はそれぞれの辺を使うか使わないかの bit 全列挙で全て求めておくと, いい感じに求まる.
ソースコード
#include <bits/stdc++.h> using namespace std; // ライブラリ部分は省略 (mint : mod 取り構造体) const int MAX = 3e5; int a[3][MAX]; int ms[3][MAX]; mint calc(int lb, int ub) { if (ub - lb == 1) return mint(a[0][lb]) * a[1][lb] * a[2][lb]; int m = (lb + ub) >> 1; mint res = calc(lb, m) + calc(m, ub); int rid[3] = { m, m, m }; int lm[3] = {}; int rm[3] = {}; mint sum[4][8]; mint all = 0; { int tm[3] = {}; for (int i = m; i < ub; i++) { for (int j = 0; j < 3; j++) { tm[j] = max(tm[j], a[j][i]); ms[j][i] = tm[j]; } } } for (int i = m; i < ub; i++) { all += mint(ms[0][i]) * ms[1][i] * ms[2][i]; } for (int l = m - 1; l >= lb; l--) { for (int i = 0; i < 3; i++) { lm[i] = max(lm[i], a[i][l]); } for (int i = 0; i < 3; i++) { while (rid[i] < ub && rm[i] < lm[i] && a[i][rid[i]] < lm[i]) { for (int j = 0; j < 8; j++) { mint v = 1; for (int k = 0; k < 3; k++) if (j & (1 << k)) { v *= ms[k][rid[i]]; } sum[i][j] += v; } rm[i] = max(rm[i], a[i][rid[i]++]); } } int vs[3] = { 0, 1, 2 }; sort(vs, vs + 3, [&](int i1, int i2) { return rid[i1] < rid[i2]; }); int s = 0; for (int i = 0; i < 4; i++) { mint kei = 1; for (int j = 0; j < 3; j++) if ((~s) & (1 << j)) { kei *= lm[j]; } res += ((i < 3 ? sum[vs[i]][s] : all) - (i ? sum[vs[i - 1]][s] : 0)) * kei; if (i < 3) s |= 1 << vs[i]; } } return res; } int main() { ios::sync_with_stdio(false), cin.tie(0); int n; cin >> n; for (int i = 0; i < n; i++) { cin >> a[0][i] >> a[1][i] >> a[2][i]; if (a[0][i] < a[1][i]) swap(a[0][i], a[1][i]); if (a[0][i] < a[2][i]) swap(a[0][i], a[2][i]); if (a[1][i] < a[2][i]) swap(a[1][i], a[2][i]); } mint res = calc(0, n) * 2 / n / (n + 1); cout << res << endl; return 0; }
感想
かなりよくできてて面白い問題だと思った.
分割統治と決めてかかってもかなり難しくて, これをまともな時間で通せる気がしない.
2017-2018 Petrozavodsk Winter Training Camp, Saratov SU Contest F. GCD
問題概要
長さ n の数列のうち、k個を削除して残ったものの GCD の最大値を求めよ.
- n <= 10^5
- k <= n / 2
- a_i <= 10^18
解法
まず, 残すもののうちの１つを乱択で決め打ちする.
k が n / 2 以下であるという制約からこれは確率 1/2 なので10回くらい試せばある程度安心.
こうして１つを決めうちする (その値を X とする) と, 他の要素は X との GCD であると見なせるので, Xの約数のみを考えればよくなる.
とりあえず X を素因数分解する必要があるが, 愚直にやると間に合わないので ミラーラビン法 + ポラード・ロー法 を使うと O(X^(1/4) logX) くらいで求まる.
あとは, 倍数のゼータ変換 (?) を行えばよくて, これは各素因数ごとに dfs して大きい方から足し合わせていけば O( (Xの約数の個数) * (X の素因数の種類数) ) くらいになる.
大きそうに見えるが, 10^18 程度の値なら約数の個数は案外小さいので, なんとか間に合う.
ソースコード
#include <bits/stdc++.h> using namespace std; using ll = long long; ll mod_prod(ll a, ll b, ll md) { ll res = (a * b - (ll)((long double)a / md * b) * md) % md; return res < 0 ? res + md : res; } ll mod_pow(ll x, ll k, ll m) { ll res = 1; while (k) { if (k & 1) res = mod_prod(res, x, m); k /= 2; x = mod_prod(x, x, m); } return res; } bool miller_rabin_test(ll n) { const static vector<ll> sprp = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31 }; if (n < 2) return false; ll d = n - 1, s = 0; while ((d & 1) == 0) d >>= 1, ++s; for (ll a : sprp) if (a <= n) { if (a == n) return true; if (mod_pow(a, d, n) != 1) { bool ok = true; for (ll r = 0; r < s; r++) { if (mod_pow(a, d << r, n) == n - 1) { ok = false; break; } } if (ok) return false; } } return true; } ll pollard_rho(ll n) { static mt19937 mt(random_device{}()); if (n == 4) return 2; ll c = mt() % n; ll x = mt() % n, y = x, d = 1; while (d == 1) { x = (mod_prod(x, x, n) + c) % n; y = (mod_prod(y, y, n) + c) % n; y = (mod_prod(y, y, n) + c) % n; d = __gcd(abs(x - y), n); } return d; } template <typename T> vector<pair<T, int>> factorize(T n) { if (n <= 1) return{}; vector<T> st = { n }, ps; while (!st.empty()) { T t = st.back(); st.pop_back(); if (t <= 1) continue; if (miller_rabin_test(t)) { ps.push_back(t); continue; } T x; while (x = pollard_rho(t), x == t); st.push_back(x); st.push_back(t / x); } sort(ps.begin(), ps.end()); vector<pair<T, int>> res; for (auto p : ps) { if (res.empty() || res.back().first != p) res.emplace_back(p, 1); else res.back().second++; } return res; } mt19937 mt(random_device{}()); int main() { int n, k; cin >> n >> k; vector<ll> a(n); for (int i = 0; i < n; i++) { cin >> a[i]; } ll res = 0; vector<int> used(n); for (int cc = 0; cc < min(n, 10); cc++) { int id; while (used[id = mt() % n]); used[id] = 1; unordered_map<ll, int> cnt; for (int i = 0; i < n; i++) { cnt[__gcd(a[i], a[id])] += 1; } auto fs = factorize(a[id]); int sz = fs.size(); for (auto pp : fs) { ll p = pp.first; function<void(int, ll)> calc = [&](int it, ll val) { if (it == sz) { for (int i = 0; i < pp.second; i++) { val *= p; } for (int i = pp.second - 1; i >= 0; i--) { cnt[val / p] += cnt[val]; if (i > 0) val /= p; } return; } if (fs[it].first == p) { calc(it + 1, val); return; } for (int i = 0; i <= fs[it].second; i++) { calc(it + 1, val); if (i + 1 <= fs[it].second) val *= fs[it].first; } }; calc(0, 1); } for (auto p : cnt) if (p.second >= n - k) { res = max(res, p.first); } } cout << res << endl; return 0; }
感想
久々にためになる問題だと思った.
10^18 くらいの値は定数個くらいなら素因数分解できるのとか, 約数の個数がある程度小さいのとかは典型っぽい気がする.
CodeChef January Lunchtime 2019 Wanderer
問題概要
N 頂点 M 辺の無向単純グラフ上を K 回移動する. 始点は頂点1. Q個の条件が次のように与えられる.
(a_i, b_i) : b_i 回目の移動後に頂点 a_i にいる.
このとき, 全ての条件を満たすような動きの場合の数 mod 1e9+7 を求めよ.
- 1 <= N, M, K, Q <= 9000
解法
dp[i][j] : i 回の移動後に頂点 j にいる場合の数
を基本的に求めていく.
各 b_i 回目の移動後に 頂点 a_i 以外の値を 0 で埋めると条件を満たすもののみ求まる.
このとき, 遷移を隣接行列の掛け算でやると O(KN^2) になってしまうが, 冷静に考えると遷移は O(M) で十分なので毎回各辺を回すだけで O(KN+KM) になる.
CodeChef SPC19F1 Lockied Away
問題概要
N が与えられるので、積がN以上となるような非負整数のマルチセットのうち, 和の最小値を求めよ.
Nは桁数が 2*10^6 以下の非負整数.
解法
4 = 2 * 2 , 5 < 2 * 3 , 6 < 3 * 3 という風に考えていくと, 4以上の値は使う必要がないことがわかる.
当然１もいらなくて (ただしN=1を除く) , 2*2*2 < 3*3 なので, 2も2個以下しか使う必要がない.
そうすると, {3,3,3, ... ,3,3} or {2,3,3, ... ,3,3} or {2,2,3, ... ,3,3} の場合のみ考えれば十分であることがわかる.
あとは long double なんかで適当にやると誤差で無理なので, 実際に整数で計算しないといけない.
値が非常に大きいので繰り返し二乗法でかつ, 10進数で持って積は FFT を使う.
桁数から大体近い3のべき乗を見積もって, 生成したあとは N を超えるまで3をかけていく.
あと, logN 回程度 FFT をするけど小さいサイズから始めて毎回 resize をしていくと計算量が O( Σ (2^i)log(2^i) ) = O(XlogX) (Xは10進数での桁数) になる.
あと NTT (高速剰余変換) を使うと誤差もなく高速でできるので, TLが異常に厳しいけどなんとか間に合った.
ソースコード
https://www.codechef.com/viewsolution/22165745
#include <bits/stdc++.h> using namespace std; using ll = long long; int mod_inv(ll a, ll m) { ll b = m, u = 1, v = 0; while (b > 0) { ll t = a / b; a -= t * b; swap(a, b); u -= t * v; swap(u, v); } return (u % m + m) % m; } ll mod_pow(ll x, ll y, ll md) { ll res = 1; while (y) { if (y & 1) res = res * x % md; x = x * x % md; y >>= 1; } return res; } template <int Mod, int PrimitiveRoot> class fast_modulo_transform { public: static vector<int> fft(vector<int> v, bool inv) { const int N = v.size(); assert((N ^ (N & -N)) == 0); int ww = mod_pow(PrimitiveRoot, (Mod - 1) / N, Mod); if (inv) ww = mod_inv(ww, Mod); for (int m = N; m >= 2; m >>= 1) { const int mh = m >> 1; int w = 1; for (int i = 0; i < mh; ++i) { for (int j = i; j < N; j += m) { const int k = j + mh; int x = v[j] - v[k]; if (x < 0) x += Mod; v[j] += -Mod + v[k]; if (v[j] < 0) v[j] += Mod; v[k] = (1LL * w * x) % Mod; } w = (1LL * w * ww) % Mod; } ww = (1LL * ww * ww) % Mod; } int i = 0; for (int j = 1; j < N - 1; ++j) { for (int k = N >> 1; k > (i ^= k); k >>= 1); if (j < i) swap(v[i], v[j]); } if (inv) { const int inv_n = mod_inv(N, Mod); for (auto& x : v) { x = (1LL * x * inv_n) % Mod; assert(0 <= x && x < Mod); } } return move(v); } static vector<int> convolution(vector<int> f, vector<int> g) { int sz = 1; const int m = f.size() + g.size() - 1; while (sz < m) sz *= 2; f.resize(sz), g.resize(sz); f = fast_modulo_transform::fft(move(f), false); g = fast_modulo_transform::fft(move(g), false); for (int i = 0; i < sz; ++i) { f[i] = (1LL * f[i] * g[i]) % Mod; } return fast_modulo_transform::fft(move(f), true); } }; const ll fmt_mod = 998244353; using fmt = fast_modulo_transform<998244353, 3>; vector<int> prod(const vector<int>& a, const vector<int>& b) { auto v = fmt::convolution(a, b); int x = v.size(); int up = 0; for (int i = 0; i < x; i++) { up += v[i]; v[i] = up % 10; up /= 10; } while (up > 0) { v.push_back(up % 10); up /= 10; } while (v.back() == 0) v.pop_back(); return move(v); } vector<int> calc_pow(int n) { vector<int> res = { 1 }, tmp = { 3 }; while (n) { if (n & 1) res = prod(res, tmp); n >>= 1; if (n) tmp = prod(tmp, tmp); } return move(res); } vector<int> get_v(string s) { vector<int> v; reverse(s.begin(), s.end()); for (auto c : s) { v.push_back(c - '0'); } return v; } void kakeru(vector<int>& v, int d) { int n = v.size(); v.push_back(0); int up = 0; for (int i = 0; i <= n; i++) { up += v[i] * d; v[i] = up % 10; up /= 10; } if (v.back() == 0) v.pop_back(); } bool comp(const vector<int>& a, const vector<int>& b) { if (a.size() != b.size()) return a.size() < b.size(); int n = a.size(); for (int i = n - 1; i >= 0; i--) if (a[i] != b[i]) { return a[i] < b[i]; } return false; } int main() { ios::sync_with_stdio(false), cin.tie(0); string S; cin >> S; int keta = S.size(); if (keta <= 6) { int N = stoi(S); if (N == 1) { cout << 1 << endl; return 0; } int res = N; for (int i = 2; i <= 4; i++) { int val = i, sum = i; while (val < N) val *= 3, sum += 3; res = min(res, sum); } cout << res << endl; return 0; } auto vec = get_v(S); int x = int(2.0959 * (keta - 1)) - 2; auto pw = calc_pow(x); int res = 1e9; for (int i = 2; i <= 4; i++) { auto v = pw; kakeru(v, i); int sum = x * 3 + i; while (comp(v, vec)) { kakeru(v, 3); sum += 3; } res = min(res, sum); } cout << res << endl; return 0; }
感想
これコンテスト中に誰も解いてなかったけど面白かった.
2018-2019 ICPC, NEERC, Northern Eurasia Finals K King Kog's Reception
問題概要
以下のクエリを Q(<=3e5) 個処理する
- join t d : 時刻 t に開始して時間 d かかるタスクが予約される (joinでのt は全て異なる)
- cancel i : i 番目のクエリで追加されたタスクを削除
- query t : 現在予約されているタスクのうち時刻 t 以前(tも含む)に開始するタスクが全て終了する時刻 - t を求めよ
ただし, タスクは同時に1つしか行えない.
解法
追加と削除のクエリがあるので, 時系列をセグ木っぽく管理して追加クエリだけにする手法が脳裏をよぎる.
が, かなり大変な上に log が２個付くので TL 2s に収めるのはかなり厳しい.
実は, このタスクは時系列上のモノイドに落とすことができる.
具体的には, 区間に暇な時間がいくらあるかと, 区間終了時に終わっていないタスクがいくらあるか のペア.
なので, ただのセグメントツリーに載せるだけで解ける.
ソースコード
https://codeforces.com/contest/1089/submission/46514491
#include <bits/stdc++.h> using namespace std; using ll = long long; namespace segment_tree { using T = pair<ll, ll>; const T id(-1, -1); T op(T a, T b) { if (a == id) return b; if (b == id) return a; ll pa, ea; tie(pa, ea) = a; ll pb, eb; tie(pb, eb) = b; if (ea == 0) ++pa; else --ea; return pb <= ea ? make_pair(pa, ea - pb + eb) : make_pair(pa + pb - ea, eb); } const int MAX = 1 << 20; T seg[MAX * 2]; void init_seg() { for (int i = MAX; i < MAX * 2; ++i) seg[i] = make_pair(0, 0); for (int i = MAX - 1; i > 0; --i) seg[i] = op(seg[i << 1], seg[(i << 1) | 1]); } void update(int p, T val) { seg[p += MAX] = val; while (p >>= 1) seg[p] = op(seg[p << 1], seg[(p << 1) | 1]); } T find(int l, int r) { T res1 = id, res2 = id; for (l += MAX, r += MAX; l < r; l >>= 1, r >>= 1) { if (l & 1) res1 = op(res1, seg[l++]); if (r & 1) res2 = op(seg[--r], res2); } return op(res1, res2); } } using namespace segment_tree; int main() { ios::sync_with_stdio(false), cin.tie(0); int q; cin >> q; vector<int> t(q); init_seg(); for (int i = 0; i < q; i++) { char c; cin >> c >> t[i]; --t[i]; if (c == '+') { int d; cin >> d; update(t[i], make_pair(0, d)); } else if (c == '-') { update(t[t[i]], make_pair(0, 0)); } else { printf("%lld\n", find(0, t[i] + 1).second); } } return 0; }
動的な Segment Tree のテクニック
区間に(初期値を)代入
普通は区間代入といえば遅延伝搬を使いますが, 遅延なしでも初期値を代入する場合は実は簡単にできます.
方法は, 代入したい区間の各部分木を null に置き換えるだけです.
一応もう少し工夫すると, 区間代入にも一般化できます.
2つの列の同じ区間 [l, r] を swap
これは一見, 平衡二分木の split 操作が必要そうに見えます.
しかしよく考えると, 代入したい区間の各部分木を swap するとできます.
永続化すればコピー操作にもできますね.
bool値の区間 flip
上記の区間 swap を使うと, 遅延伝搬なしでbool値の区間 flip 操作が可能になります.
やり方は, 元の列と各値が反転した列の2つの列を持っておいて, 反転操作ではその区間でswapするとできます.
使える問題：AOJ RUPC2017 Day3 E
さらに, 整数に対してビット独立にこれを使うと区間XORができることもあります.
使える問題：HackerRank HourRank28 C
空間計算量を O(NlogX) から O(N) に減らす
基数木(Radix Tree/Patricia Trie)でも同等のことができますが, こっちの方が若干実装が楽そう(?)
点更新の場合
通常は null のノードにたどり着いても, 葉までノードを作りながら降りていきますが, これをすぐにやめるようにします.
各ノードには持たせるものは, 位置とその位置の値と部分木全体の値です.
更新時には更新される位置を持ったノードが見つかるまで降りていって更新しますが, 見つからなければ新しくノードを作成します.
取得クエリでは普通と同じように降りていきますが, 途中のノードも目的の区間に含まれていれば答えと合わせます.
あと, ノードの重複を避けるために適当な条件を満たすようにしておきましょう.
この方法を使うと平均的に深さが減ったり, メモリ確保の回数も減るので定数倍が速くなることが多い気がします.
実装例
template <typename M> // Mは 型type, 単位元id(), 演算op() を持つモノイド class dynamic_segment_tree { using T = typename M::type; struct node { int pos; T val, all; node *l, *r; node(int p, T v) : pos(p), val(v), all(v), l(nullptr), r(nullptr) {} }; private: const int n; node *root; public: dynamic_segment_tree(int n_) : n(1 << (int)ceil(log2(n_))), root(nullptr) {} void update(int p, T val) { root = change(root, p, val, 0, n); } T find(int l, int r) { return get(l, r, root, 0, n); } private: T value(node *t) { return t ? t->all : M::id(); } T get(int l, int r, node* t, int lb, int ub) { if (!t || ub <= l || r <= lb) return M::id(); if (l <= lb && ub <= r) return t->all; int c = (lb + ub) / 2; T res = get(l, r, t->l, lb, c); if (l <= t->pos && t->pos < r) res = M::op(res, t->val); return M::op(res, get(l, r, t->r, c, ub)); } node *change(node* t, int p, T val, int lb, int ub) { if (!t) return new node(p, val); if (t->pos == p) { t->val = val; t->all = M::op(value(t->l), M::op(t->val, value(t->r))); return t; } int c = (lb + ub) / 2; if (p < c) { if (p > t->pos) swap(p, t->pos), swap(val, t->val); t->l = change(t->l, p, val, lb, c); } else { if (p < t->pos) swap(p, t->pos), swap(val, t->val); t->r = change(t->r, p, val, c, ub); } t->all = M::op(value(t->l), M::op(t->val, value(t->r))); return t; } };
使える問題：AtCoder KUPC2018 M
(↑この問題の想定解はこのテクなんですが, 使わなくても通っちゃうみたいです. 悲しいね.)
ICPC Seoul Regional 2018 参加記
はじめに
ICPC ソウル地区大会に Zerokan_Sunshine というチームで参加してきました.
結果は 12 問中 7 完で 17 位でした.
10 完早解きじゃないと WF にいけないらしくてとてもつらい.
本番の流れ
うろ覚えで書いてるのでたまに時系列おかしいかも.
(0:00-1:00)
- 開始直後に問題(なぜか1部しか配られなかった)をバラバラにして眺めるとDがタイピング問題っぽいのでsuibakaに投げる
- suibakaと一緒に問題を眺めてるとサンプルがおかしいことに気付く
- 問題文を読むと, サンプル output が間違っているという結論に至ったのでそのままsuibakaが書くとAC
- えーなんだこれ
- 次にCが見覚えのある問題っぽかったのでよく読まずに書き始める
- 書き終えてサンプルが合わないのでちゃんと読み直すと, 順序を並び替えて良いことに気付く
- 小大小大…みたいな感じでやると良いんでは？とか適当に嘘解法を生やしたけどサンプルすら合わないので一旦保留
- nakanoから雑にH, E, Kの概要を聞く
- Hは明らかに不可能かつ幾何なので後回し
- Eはmaxの部分をなぜかsumだと勘違いする
- 軽く考えたけど解けそうない
- Kは条件が3つずつあるのと混乱して, 色も3つあると勘違いする
- 2色なら2-SATだけどなあ・・・3色ならわからん, となる
(1:00-2:00)
- suibaka に「書ける問題あったら投げて」と言われるけど無いと答える
- PCが空いてることに耐えかねたsuibakaがFの構文解析を書くけどWA
- 聞いた問題の中で解けそうな問題がなくて, ここから長時間虚無になる
- 考察が全然進まないから昼ご飯食べてた気がする
- nakanoがsuibakaにKを伝えてるのを聞いていると, 色が2色しかないことに気付いたので, 書き始める
- 2-SATなのでSCCのライブラリをsuibakaに要求すると, 2-SATのライブラリ(復元付き)が渡されたので写経しまくる
- 写経失敗したけどsuibakaに修正してもらう
- やっと2個目のAC
(2:00-3:00)
(3:00-4:00)
(4:00-5:00)
- Cをnakanoに伝えると, 右端を全通り試して残りは左から小大小大…をするといけるかも？と言われたので書く
- 嘘っぽいけど流石にサンプル１くらいは合うはずだが, バグりまくって全く合わないので絶望
- この程度の処理でこんなにバグるのヤバくないか？？？
- nakanoからEはにぶたんで解けると言われて, そもそもmaxをsumと勘違いしてることに気付く
- それなら解けるじゃん・・・
- suibakaがめちゃくちゃJを書いていたので, 交代交代で書くけどどっちも間に合わないまま終了
感想
いくらなんでも勘違いが多すぎるのでなんとかしたい.
でもうまくいってたとしても9完どまりくらいな気がするし, やっぱり韓国勢はヤバい.
