AtCoder 全国統一プログラミング王決定戦本戦/NIKKEI Programming Contest 2019 H - Homework Scheduling
問題概要
$ N $個の宿題があり、1日に1つずつ終わらせていく。各宿題には期限$ A_i $ が設定されており、宿題$ i $をDay $ A_i $までに終わらせると$ X_i $点、それ以後に終わらせると$ Y_i $点が得られる。$ 1\leq k \leq N $なる$ k $について、Day 1 ~ Day $ k $で得られる得点を最大化せよ。
制約
- $1\leq N\leq 2\times 10^5$
- $1 \leq A_i \leq N $
- $ 1\leq Y_i < X_i \leq N $
解法
はじめに$A_i$が昇順になるように宿題を並べ替える。
最小費用流への言い換え
グラフの構築
この問題のDay $k$に関する答えは、次のグラフにおける流量$ k $の$S\to T$最小費用流のコストの-1倍である。
ただし、$ Q_i \to T$の辺たちは次の方法によってはられる。
「$i = 1,2,\dots,N$のそれぞれに対して、$A_k< i$となる$ k $の個数を$c_i$とする。$Q_{c_i+1} \to T$に容量1の辺をはる。」
$ Q_i \to T$の辺に左から順に$r_1,\dots,r_n$と名前をつけておく。
グラフの例
問題の例2について、$\{A_i\}_{i=1}^6=\{1,2,2,3,4,4\}$なので$\{c_i\} _{i=1}^6= \{1,3,4,6,7,7\}$となる。 グラフは下のようになる。
グラフの最小費用流へ帰着できることの説明
グラフ上での$ S\to T$パスは、
( $j \leq i$)
- $S\to P_i \to Q_{N+1} \to Q_j \overset{r_k}{\rightarrow} T$ ( $j\leq N+1$)
のどちらかの形になる。
前者を宿題$i$を$k$日目に終わらせ、かつ$X_i$点を獲得した（すなわち$k\leq A_i$）ことに対応させる。また、後者を宿題$i$を$k$日目に終わらせ、かつ$Y_i$点を獲得した（すなわち$k> A_i$）ことに対応させる。 従って、下図青のパスは宿題2を1日目に終わらせ、$X_2$点を得たことに対応し、赤のパスは宿題3を4日目に終わらせ、$Y_3$点を得たこと に対応している。
ここで、
(A)任意のDay $k$までの宿題の終わらせ方がグラフ上の容量$k$のフローで表現でき、かつ( 宿題の得点の和 )＝( フローのコスト )$\times (-1) $であること
逆に(B)グラフ上の任意の容量$k$のフローがDay $k$までの宿題の終わらせ方を表し、かつ( 宿題の得点の和 )＝( フローのコスト )$\times (-1) $であること
がわかれば、最小費用流へ問題を言い換えてよいと言える。それぞれが正しいことは次のように考えるとわかる。
(A) Day $k$までのある宿題の終わらせ方が与えられたとする。日ごとに終わらせた宿題をグラフ上の$S\to T$パスに置き換えて、$k$本の$S\to T$パスを作る。$r_1,\dots,r_n$の作り方から、$k$本のパスに流量1を流したフローがグラフの容量制約をみたすことがわかる。
(B) グラフ上の容量$k$のフローが与えられたとする。辺$r_1,\dots,r_n$が左から優先して使われるようにパスを左に詰める。たとえば上図の容量2のフローの場合、赤いパスをにずらす。この変更はフローのコストを変えない。変更を行った後でそれぞれのパスを宿題の終わらせ方に置き換える。辺$r_1,\dots,r_n$の作り方から、制約をみたす宿題の終わらせ方であることがわかる。
最小費用流問題の解き方
ad-hocに工夫したPrimal-Dualによって解く。 （Primal-Dualの参考 Spaghetti Source - 最小費用流 (Primal-Dual) 最小費用流(Primal-Dual) | Luzhiled’s memo）
Primal-Dualは、残余グラフ上の最小コストの$S\to T$パスに流量1のフローを流す操作を繰り返していく方法である。この問題ではグラフの構造から最小コストの$S\to T$パスを求める部分が高速にできるのがポイントである。
なお、editorial中にある仮定
- 辺 Qi → T の容量が正であるすべての i に対して、辺 Qi−1 → Qi の容量が 0 である。
というのは、前項とは逆に辺$r_1,\dots,r_n$が右から優先して使われるようにパスを右に詰めることを意味している。この記事中でもそのように辺$r_1,\dots,r_n$を右から使うことにする。
残余グラフ上の$S\to T$パスは次の4種類に分けられる。
- A :
( $j \leq i$)
- B : $S\to P_i \to Q_{N+1} \to Q_j \overset{r_l}{\rightarrow} T$ ( $j\leq N+1$)
- C : $S\to P_i \to Q_i \to Q_j \to P_j \to Q_{N+1} \to Q_k \overset{r_l}{\rightarrow} T$ ( $i> j$)
- D : $S\to P_i \to Q_i \to Q_j \to P_j \to Q_{N+1} \to Q_k \overset{r_l}{\rightarrow} T$ ( $i< j$)
仮定の部分でも述べたように、$r_l$は残余グラフ上で空いている中でできるだけ右、つまり$l$の大きいものを使う。従って、A~Dのタイプの中での最小コストは$i,j$を決めることで一意に定まる。
C, Dタイプの最小コストを与える$i,j$は、解説にあるSegment Treeを実装すると求められる。Cタイプの最小コストは$C_i+D_j ~(i>j)$の最小値である。Dタイプの最小コストは$m\neq 0$ ならば単に$C_i+D_j ~(i<j)$の最小値であり、$m=0$ならば$C_i+D_j ~(i<j, E_i,E_j >m)$の最小値である。
Aタイプについては，点更新と区間minクエリに対応したSegmentTreeを用意して、 のパスがあるなら$seg[i]=$そのコスト、ないなら$seg[i]=\infty$となるように更新していけば求められる。Bタイプも同様。
実装の際には、辺$r_1,\dots,r_n$を右から順に使うという制約をみたすため、残余グラフ上の$ Q_i \to T$の辺のうち容量が正のものを管理しておく必要がある。
アルゴリズム
$1\leq k\leq N$について，
- A,B,C,Dそれぞれのタイプでの最小コストの$S\to T$パスを求める
- その中で最小コストの$S\to T$パスに流量1を流し、残余グラフを更新する
- 現時点でのコストの(-1)倍を答える
を繰り返す。
解答例
#include <bits/stdc++.h> #define DBG(...) ; using LL = long long; constexpr LL LINF = 334ll << 50; constexpr int INF = 15 << 26; namespace Problem { using namespace std; struct Min_Id { LL val; int id; Min_Id() : val(LINF), id(-1){}; Min_Id(LL v, const int id) : val(v), id(id){}; bool operator<(const Min_Id &b) const { return val < b.val; } bool operator==(const Min_Id &b) const { return val == b.val; } }; struct Min_IdLR { LL val; int idl, idr; Min_IdLR() : val(LINF * 2), idl(INF), idr(-1){}; Min_IdLR(LL v, const int l, const int r) : val(v), idl(l), idr(r){}; bool operator<(const Min_IdLR &b) const { return val < b.val; } bool operator==(const Min_IdLR &b) const { return val == b.val; } tuple<const LL &, const int &, const int &> ref() const { return tie(val, idl, idr); } }; struct Min_Range { LL val; int idl, idr; Min_Range() : val(LINF), idl(INF), idr(-1){}; Min_Range(LL v, const int l, const int r) : val(v), idl(l), idr(r){}; bool operator<(const Min_Range &b) const { return val < b.val; } bool operator==(const Min_Range &b) const { return val == b.val; } }; struct Node { Min_Range E; Min_Id CL, CR, DL, DR; Min_IdLR Dc, CdMinE, CdNoMinE; Node() : E(), CL(), CR(), DL(), DR(), Dc(), CdMinE(), CdNoMinE(){}; }; struct Operator { LL addC, addD; int addE; Operator() : addC(0), addD(0), addE(0){}; Operator(LL c, LL d, int e) : addC(c), addD(d), addE(e){}; }; struct RangeNodeMin { using T = Node; T operator()(const T &a, const T &b) { T ret = T(); if (a.E < b.E) { ret.E = a.E; ret.CL = a.CL; ret.CR = min({a.CR, b.CL, b.CR}); ret.DL = a.DL; ret.DR = min({a.DR, b.DL, b.DR}); auto CdMinETrans = Min_IdLR(a.CL.val + min(b.DL, b.DR).val, a.CL.id, min(b.DL, b.DR).id); ret.CdMinE = min(a.CdMinE, CdMinETrans); auto CdNoMinETrans = Min_IdLR(a.CR.val + min(b.DL, b.DR).val, a.CR.id, min(b.DL, b.DR).id); ret.CdNoMinE = min({a.CdNoMinE, b.CdNoMinE, b.CdMinE, CdNoMinETrans}); } else if (b.E < a.E) { ret.E = b.E; ret.CL = min({a.CL, a.CR, b.CL}); ret.CR = b.CR; ret.DL = min({a.DL, a.DR, b.DL}); ret.DR = b.DR; auto CdMinETrans = Min_IdLR(b.DR.val + min(a.CL, a.CR).val, min(a.CL, a.CR).id, b.DR.id); ret.CdMinE = min(b.CdMinE, CdMinETrans); auto CdNoMinETrans = Min_IdLR(b.DL.val + min(a.CL, a.CR).val, min(a.CL, a.CR).id, b.DL.id); ret.CdNoMinE = min({a.CdNoMinE, a.CdMinE, b.CdNoMinE, CdNoMinETrans}); } else { ret.E = a.E; ret.E.idr = b.E.idr; ret.CL = min({a.CL, a.CR, b.CL}); ret.CR = b.CR; ret.DL = a.DL; ret.DR = min({a.DR, b.DL, b.DR}); auto CdMinETrans = Min_IdLR(a.CL.val + min(b.DL, b.DR).val, a.CL.id, min(b.DL, b.DR).id); auto CdMinETrans2 = Min_IdLR(b.DR.val + min(a.CL, a.CR).val, min(a.CL, a.CR).id, b.DR.id); ret.CdMinE = min({a.CdMinE, b.CdMinE, CdMinETrans, CdMinETrans2}); auto CdNoMinETrans = Min_IdLR(a.CR.val + b.DL.val, a.CR.id, b.DL.id); ret.CdNoMinE = min({a.CdNoMinE, b.CdNoMinE, CdNoMinETrans}); } Min_IdLR minDcTrans = Min_IdLR(min(a.DL.val, a.DR.val) + min(b.CL.val, b.CR.val), min(a.DL, a.DR).id, min(b.CL, b.CR).id); ret.Dc = min({a.Dc, b.Dc, minDcTrans}); return ret; } static T id() { return T(); }; }; struct RangeAdd { using T = Operator; T operator()(const T &a, const T &b) { Operator ret = a; ret.addC += b.addC; ret.addD += b.addD; ret.addE += b.addE; return ret; }; static T id() { return T(); }; }; struct Min_Add { using V = RangeNodeMin; using A = RangeAdd; V::T operator()(const V::T &a, const A::T &b, const int &h) { V::T ret = a; ret.CL.val += b.addC; ret.CR.val += b.addC; ret.DL.val += b.addD; ret.DR.val += b.addD; ret.E.val += b.addE; return ret; }; }; struct RangeMin { using T = pair<long long, int>; T operator()(const T &a, const T &b) { return min(a, b); } static constexpr T id() { return T(LINF, -1); }; }; struct RangeSet { using T = long long; T operator()(const T &a, const T &b) { return b; }; static constexpr T id() { return T(LINF * 2); }; }; struct Min_Set { using V = RangeMin; using A = RangeSet; V::T operator()(const V::T &a, const A::T &b, const int &h) { if (b == LINF * 2) return a; else return V::T(b, a.second); }; }; template <typename M> struct LazySegmentTree { using Val = typename M::V::T; using Del = typename M::A::T; typename M::V calc; typename M::A composite; M act; vector<Val> val; vector<Del> delay; int n, height; explicit LazySegmentTree(int size) : n(1), height(0) { while (n < size) { n <<= 1; ++height; } val = vector<Val>(n << 1, calc.id()); delay = vector<Del>(n << 1, composite.id()); } explicit LazySegmentTree(const vector<Val> &v) : n(1), height(0) { while (n < v.size()) { n <<= 1; ++height; } val = vector<Val>(n << 1, calc.id()); delay = vector<Del>(n << 1, composite.id()); copy(v.begin(), v.end(), val.begin() + n); for (int i = n - 1; i > 0; --i) { val[i] = calc(val[i * 2], val[i * 2 + 1]); } } void propagate(int k, int h) { val[k] = act(val[k], delay[k], h); if (h > 0) { delay[2 * k] = composite(delay[2 * k], delay[k]); delay[2 * k + 1] = composite(delay[2 * k + 1], delay[k]); } delay[k] = composite.id(); } inline void update(int l, int r, Del v) { update(l, r, v, 1, 0, n, height); } void update(int l, int r, Del v, int k, int a, int b, int h) { if (l <= a && b <= r) { delay[k] = composite(delay[k], v); propagate(k, h); return; } propagate(k, h); if (b <= l || r <= a) { return; } // update(l, r, v, 2 * k, a, (a + b) / 2, h - 1); update(l, r, v, 2 * k + 1, (a + b) / 2, b, h - 1); val[k] = calc(val[2 * k], val[2 * k + 1]); } inline Val query(int l, int r) { return query(l, r, 1, 0, n, height); } Val query(int l, int r, int k, int a, int b, int h) { if (b <= l || r <= a) { return calc.id(); } // propagate(k, h); if (l <= a && b <= r) { return val[k]; } Val vall = query(l, r, 2 * k, a, (a + b) / 2, h - 1); Val valr = query(l, r, 2 * k + 1, (a + b) / 2, b, h - 1); return calc(vall, valr); } }; struct Edge { int to, cap; LL cost; int rev; }; class Solver { public: int n; vector<LL> a, x, y; LazySegmentTree<Min_Add> seg; //解説中にあるSegment Tree C,Dタイプの最小コストパスを求めるために使う LazySegmentTree<Min_Set> segA, segB; //A,Bタイプ multiset<int> Qpos; //残余グラフ上の$ Q\_i \to T$容量と同じ個数の i を含む多重集合 Solver(LL n) : n(n), a(n), x(n), y(n), seg(0), segA(0), segB(0){}; void solve() { input(); makeGraph(); //vertex //P:0~n-1 (P_i = i) //Q:n~2*n (Q_i = n + i) //S:2*n+1,T:2*n+2 //path //A:S->P->Q->Q->T //B:S->P->Q_{n}->Q->T //C:S->P->Q_i->Q_j->P->Q_{n}->Q->T (i>j) //D:S->P->Q_i->Q_j->P->Q_{n}->Q->T (i<j) LL mincost = 0; for (int i = 1; i <= n; ++i) { auto path = minCostAugment(); mincost += augmentPath(path); cout << -mincost << endl; } } vector<int> minCostAugment() { int idA = 0, idB = 0, idCl = 0, idCr = 0, idDl = 0, idDr = 0; LL costA = LINF, costB = LINF, costC = LINF, costD = LINF; tie(costA, idA) = segA.query(*Qpos.begin(), n); tie(costB, idB) = segB.query(0, n); Node node = seg.query(0, n + 1); tie(costC, idCl, idCr) = node.Dc.ref(); if (node.E.val == 0) { tie(costD, idDl, idDr) = node.CdNoMinE.ref(); } else { tie(costD, idDl, idDr) = min(node.CdMinE, node.CdNoMinE).ref(); } vector<int> ret; if (min({costA, costB, costC, costD}) == costA) { ret = {2 * n + 1, idA, n + idA}; //DBG('A', ret) } else if (min({costA, costB, costC, costD}) == costB) { ret = {2 * n + 1, idB, n + n}; //DBG('B', ret) } else if (min({costA, costB, costC, costD}) == costC) { ret = {2 * n + 1, idCr, idCr + n, idCl + n, idCl, n + n}; //DBG('C', ret) } else { ret = {2 * n + 1, idDl, idDl + n, idDr + n, idDr, n + n}; //DBG('D', ret) } connecttoT(ret); return ret; } LL augmentPath(vector<int> &p) { LL ret = 0; for (int i = 0; i < (int)p.size() - 1; ++i) { int u = p[i], v = p[i + 1]; if (u == 2 * n + 1) { //S->P segA.update(v, v + 1, LINF); segB.update(v, v + 1, LINF); seg.update(v, v + 1, Operator(x[v] + LINF, 0, 0)); } else if (v == 2 * n + 2) { //Q->T Qpos.erase(Qpos.lower_bound(u - n)); } else if (0 <= u && u < n && v == n + n) { //P -> Q_n ret -= y[u]; } else if (0 <= u && u < n && n <= v && v < n + n) { //P->Q seg.update(u, u + 1, Operator(0, -LINF / 2 + x[u] - y[u], 0)); ret -= x[u]; } else if (0 <= v && v < n && n <= u && u < n + n) { //Q->P seg.update(v, v + 1, Operator(0, +LINF / 2 - x[v] + y[v], 0)); ret += x[v]; } else if (n <= u && u <= n + n && n <= v && v <= n + n && u < v) { //Q_i->Q_j (i<j) seg.update(u - n, v - n, Operator(0, 0, -1)); } else if (n <= u && u <= n + n && n <= v && v <= n + n && u > v) { //Q_i->Q_j (i>j) seg.update(v - n, u - n, Operator(0, 0, 1)); } else if (n <= u && u <= n + n && n <= v && v <= n + n && u == v) { //Q_i->Q_i } else { assert(false); } } return ret; } void connecttoT(vector<int> &v) { auto i = *prev(Qpos.upper_bound(v.back() - n)); v.push_back(i + n); v.push_back(2 * n + 2); } void input() { vector<vector<LL>> tmp(n, vector<LL>(3)); for (int i = 0; i < n; ++i) { cin >> tmp[i][0] >> tmp[i][1] >> tmp[i][2]; } sort(tmp.begin(), tmp.end()); for (int i = 0; i < n; ++i) { a[i] = tmp[i][0] - 1; x[i] = tmp[i][1]; y[i] = tmp[i][2]; } } void makeGraph() { for (int i = 0; i < n; ++i) { auto lb = distance(a.begin(), lower_bound(a.begin(), a.end(), i)); Qpos.insert(lb); } vector<Node> nodes(n); vector<pair<LL, int>> costA(n), costB(n); for (int i = 0; i < n; ++i) { nodes[i].E = Min_Range(0, i, i); nodes[i].CL = Min_Id(-x[i], i); nodes[i].CR = Min_Id(LINF, -1); nodes[i].DL = Min_Id(LINF / 2, i); nodes[i].DR = Min_Id(LINF, -1); nodes[i].Dc = Min_IdLR(LINF * 2, -1, -1); nodes[i].CdMinE = Min_IdLR(LINF * 2, -1, -1); nodes[i].CdNoMinE = Min_IdLR(LINF * 2, -1, -1); costA[i] = {-x[i], i}; costB[i] = {-y[i], i}; } seg = LazySegmentTree<Min_Add>(nodes); segA = LazySegmentTree<Min_Set>(costA); segB = LazySegmentTree<Min_Set>(costB); } }; } } // namespace Problem int main() { std::cin.tie(0); std::ios_base::sync_with_stdio(false); long long n = 0; std::cin >> n; Problem::Solver sol(n); sol.solve(); return 0; }
note
新たなモノイドが手に入った。
AtCoder 全国統一プログラミング王決定戦予選/NIKKEI Programming Contest 2019 E - Weights on Vertices and Edges
問題概要
$ N $頂点$ M $辺のグラフ $ G $があり、頂点と辺には重みが与えられている。このグラフから辺を削除して、次の条件を満たすようにしたい。
条件：「削除されていない任意の辺について、その辺を含む連結成分の頂点の重みの総和が、その辺の重み以上である。」
削除する必要がある辺は最小で何本か。
制約
- $1\leq N\leq 10^5$
- $N-1\leq M\leq 10^5$
- グラフは連結
解法
辺が0本削除された状態で条件をみたしていない辺があるとする。グラフの他の辺が削除されるかどうかにかかわらず、その辺は絶対に削除しなければならない。辺を削除したことで新たに条件に違反する辺が現れたなら、新たに違反した辺も削除する必要がある。このように考えていくと、条件に違反する辺から次々と削除していき、違反する辺がなくなるまでに削除された本数を答えればよいことになる。
辺を1つ削除するたびにすべての辺について条件をみたすか確認しなければならないように思えるが、重みが大きいほど条件に抵触しやすいため重みの降順に一巡するだけでよい。
従って、重みの降順に辺を見ていき、
- 見ている辺Aが条件をみたしているのかを判定する
- みたしているなら辺Aを残し、みたしていないなら辺Aを削除する
をという操作を順次行えば解ける。これは
- 頂点が属する連結成分内の頂点重みの総和を答える
- 辺を削除したときに連結成分を分割する
の操作ができるデータ構造があれば可能。
この問題では、ある辺が削除されるときには連結成分内のそれより重い辺がすべて削除されている。そのため「$ N $頂点に$ G $の各頂点の重みを持たせた部分永続Union-Findを用意して、$ M $本の辺を軽い方から順にUniteしたもの」は要求をみたす。永続ではなく、操作をUndoできるUnion-Findでも十分。
解答例
#include <bits/stdc++.h> using LL = long long; constexpr LL LINF = 334ll << 53; namespace Problem { using namespace std; struct UnionFind { vector<pair<LL, LL>> memo; vector<LL> par; UnionFind(int n) { par = vector<LL>(n, -1); } LL root(int x) { if (par[x] < 0) return x; return par[x] = root(par[x]); } LL reversible_root(int x) { if (par[x] < 0) return x; return reversible_root(par[x]); } void unite(int x, int y) { x = root(x); y = root(y); if (x != y) { if (par[x] < par[y]) swap(x, y); par[y] += par[x]; par[x] = y; } return; } void reversible_unite(int x, int y) { x = reversible_root(x); y = reversible_root(y); if (x != y) { if (par[x] < par[y]) swap(x, y); memo.emplace_back(x, par[x]); memo.emplace_back(y, par[y]); par[y] += par[x]; par[x] = y; } return; } bool same(int x, int y) { return root(x) == root(y); } LL reversible_same(int x, int y) { return reversible_root(x) == reversible_root(y); } LL size(int x) { return -par[root(x)]; } LL reversible_size(int x) { return -par[reversible_root(x)]; } void revert_all() { for (int i = memo.size() - 1; i >= 0; --i) { par[memo[i].first] = memo[i].second; } memo.clear(); } void revert(int i) { par[memo[i].first] = memo[i].second; } }; class Solver { public: int n, m; UnionFind uf; vector<pair<pair<int, int>, pair<LL, int>>> edges; Solver(LL n, LL m) : n(n), m(m), uf(n){}; void solve() { //UFに頂点の重みを持たせる for (int i = 0; i < n; ++i) { int tmp; cin >> tmp; uf.par[i] = -tmp; } for (int i = 0; i < m; ++i) { int a, b, y; cin >> a >> b >> y; --a; --b; edges.push_back({{a, b}, {y, -1}}); } sort(edges.begin(), edges.end(), [&](auto xx, auto yy) { return xx.second.first < yy.second.first; }); //軽い辺から連結成分を繋いで、uniteしたらメモする for (int i = 0; i < m; ++i) { int a = edges[i].first.first; int b = edges[i].first.second; if (uf.reversible_root(a) != uf.reversible_root(b)) { edges[i].second.second = uf.memo.size(); uf.reversible_unite(a, b); } } //条件に違反する辺を削除して、uniteを巻き戻す int ans = 0; for (int i = m - 1; i >= 0; --i) { if (uf.reversible_size(edges[i].first.first) < edges[i].second.first) { ans++; if (edges[i].second.second != -1) { uf.revert(edges[i].second.second); uf.revert(edges[i].second.second + 1); } } } cout << ans << endl; } }; } // namespace Problem signed main() { std::cin.tie(0); std::ios_base::sync_with_stdio(false); long long n = 0, m; std::cin >> n >> m; Problem::Solver sol(n, m); sol.solve(); return 0; }
note
AtCoder Educational DP Contest V Subtree
問題概要
$ N $頂点の木(辺は$(x_1,y_1),\dots,(x_{N-1},y_{N-1})$)の頂点を白と黒に塗り分ける。ただし、黒い頂点どうしは相互に黒い頂点のみを通って到達可能でなければならない。そのような塗り方のうち、頂点$i$が黒く塗られるものは何通りあるか。$1\leq i\leq n$のそれぞれについて、$ M $で割った余りを答えよ。
制約
- $1\leq N\leq 10^5$
- $2\leq M\leq 10^9$
解法
特定の頂点$i$について答えるなら、頂点$i$を根とした木DPによって$\mathcal O(N)$で答えられる。この問題ではすべての$i$について答えなければならないが、全方位木DPをすればやはり$\mathcal O(N)$でできる。
辺で結ばれた頂点$i,j$について、
$dp[i][j]=$頂点$i$から見た$j$の部分木（色付きの部分）を塗る組み合わせのうち頂点$j$を黒く塗る方法の数
とする。実装上はこの値を辺の情報として持たせている。
1回目のDFS
適当な頂点を根とする。頂点$v$の親を$p$,子を$c_1,\dots,c_r$ ($ r $は$ v $の次数-1)とする。頂点$v$を黒く塗るとき、頂点$c_i$以下の部分木を塗る方法は、
- $c_i$を黒く塗る：$ dp[v][c_i] $通り
- $c_i$を白く塗る：1通り
となる。よってこのときの遷移は、
\begin{align} dp[p][v] = \prod_{k\in\{1,\dots,r\}}(dp[v][c_k]+1) \end{align}
2回目のDFS
頂点$v$に隣接する頂点を$u_1,\dots,u_q$ ($ q $は$ v $の次数)とする。$dp[v][u_1],\dots,dp[v][u_q]$が求まっていれば、
\begin{align} dp[u_i][v] = \prod_{k\in\{1,\dots,q\}-\{i\}}(dp[v][u_k]+1) \end{align}
で求められる。これを求めるときに \begin{align} &\prod_{k\in\{1,\dots,q\}-\{i\}}(dp[v][u_k]+1)\\ =& \frac{\prod_{k\in\{1,\dots,q\}}(dp[v][u_k]+1)}{dp[v][u_i]+1}\\ =&\prod_{k\in\{1,\dots,q\}}(dp[v][u_k]+1)\times (dp[v][u_i]+1)^{-1}~({\rm mod}~ M) \end{align} として計算しようとすると、$(dp[v][u_i]+1)=0 ~({\rm mod}~M)$のときに$0^{-1}~({\rm mod}~M)$が現れてしまい破綻する。この問題を回避するには、
\begin{align} &\prod_{k\in\{1,\dots,q\}-\{i\}}(dp[v][u_k]+1)\\ =&\prod_{k\in\{1,\dots,i-1\}}(dp[v][u_k]+1)\times \prod_{k\in\{i+1,\dots,q\}}(dp[v][u_k]+1) \end{align}
と考えればよい。前後両側からの累積積を前計算しておけば、$\prod_{k\in\{1,\dots,i-1\}}(dp[v][u_k]+1),\prod_{k\in\{i+1,\dots,q\}}(dp[v][u_k]+1)$を除算なしで求められる。
解答例
#include <bits/stdc++.h> using LL = long long; namespace Problem { using namespace std; class Solver { public: int n; LL mod; vector<LL> ans; struct Edge { int to, rev; LL val; }; vector<vector<Edge>> t; Solver(LL n, LL m) : n(n), mod(m), ans(n), t(n){}; void solve() { for (int i = 0; i < n - 1; ++i) { int a, b; cin >> a >> b; --a; --b; t[a].push_back({b, (int)t[b].size(), 1}); t[b].push_back({a, (int)t[a].size() - 1, 1}); } if (n == 1) { cout << 1 << endl; return; } dfs(0); dfs2(0); for (int i = 0; i < n; ++i) { cout << ans[i] << endl; } } void dfs(int v, int p = -1) { LL res = 1; int par; for (int i = 0; i < (int)t[v].size(); ++i) { if (t[v][i].to != p) { dfs(t[v][i].to, v); res *= t[v][i].val + 1; res %= mod; } else { par = i; } } if (p == -1) { ans[v] = res; } else { t[p][t[v][par].rev].val = res; } } void dfs2(int v, int p = -1) { LL res = 1; for (int i = 0; i < t[v].size(); ++i) { res *= t[v][i].val + 1; res %= mod; } ans[v] = res; //累積積の計算 vector<LL> prodl(t[v].size(), 1), prodr(t[v].size(), 1); prodl[0] = (t[v][0].val + 1) % mod; for (int i = 1; i < t[v].size(); ++i) { prodl[i] = prodl[i - 1] * (t[v][i].val + 1) % mod; } prodr[t[v].size() - 1] = (t[v][t[v].size() - 1].val + 1) % mod; for (int i = (int)t[v].size() - 2; i >= 0; --i) { prodr[i] = prodr[i + 1] * (t[v][i].val + 1) % mod; } for (int i = 0; i < t[v].size(); ++i) { if (t[v][i].to != p) { if (i > 0) { t[t[v][i].to][t[v][i].rev].val *= prodl[i - 1]; t[t[v][i].to][t[v][i].rev].val %= mod; } if (i < (int)t[v].size() - 1) { t[t[v][i].to][t[v][i].rev].val *= prodr[i + 1]; t[t[v][i].to][t[v][i].rev].val %= mod; } dfs2(t[v][i].to, v); } } } }; } // namespace Problem int main() { std::cin.tie(0); std::ios_base::sync_with_stdio(false); // std::cout << std::fixed << std::setprecision(12); long long n = 0, m; std::cin >> n >> m; Problem::Solver sol(n, m); sol.solve(); return 0; }
note
Educational DP Contest のうち、自分にとって最も教育的だった問題。
コンテスト中は「 $ M $が素数でなくて面倒だな。合成数でmodを計算させる教育的問題なのかな？」と思っていた。だが、この問題で注意すべきなのは2回目のDFSで$0^{-1}$をかける計算が発生しないようにすること。それは$ M $の値が素数でも合成数でも同じ。
codeforces Good Bye 2018 E. New Year and the Acquaintance Estimation
問題概要
頂点のグラフ$G$のうち、
頂点の次数が数列$a_1,\dots,a_N$として与えられている。$G$が単純無向グラフであるとき、残り1頂点の次数$a_{0}$としてありうるものをすべて答えよ。
制約
- $1\leq N\leq 5\times10^5$
- $0 \leq a_i \leq N$
解法
問題文に貼られたリンク Graph realization problem - Wikipedia から Erdős–Gallai theorem - Wikipedia を見つけることが重要だった。
Erdős–Gallai theorem
$d_0\geq\dots\geq d_N$がある$N+1$頂点単純無向グラフの次数の列となっていることは、次の2条件を満たすことと同値：
- $\displaystyle\sum_{i=0}^N d_i \equiv 0~({\rm mod}~2)$
- $\displaystyle\sum_{i=0}^k d_i\leq k(k+1) + \sum_{i=k+1}^N \min(d_i,k+1) ~~(k = 0,\dots, N)$.
Editorial解
Erdős–Gallai theoremを$a_{0}=0,\dots,N$についてすべて調べると$\mathcal O(N^2)$になってしまう。しかし実は$X<Y$の2数が答えに含まれるとき、$Z \equiv X~({\rm mod}~2)$かつ$X<Z<Y$なる$Z$は答えになる（Good Bye 2018 — Editorial - Codeforces）。従って二分探索の要領で$\mathcal O(\log N)$個の$a_{0}$について、Erdős–Gallai theoremの条件が成立するか確かめれば解ける。$a_{0}$を固定したとき$\mathcal O(N)$で確かめられるので、全体で$\mathcal O(N\log N)$。
Segment tree 解
答えが区間になることが分かっていなくても解ける。$a_0,\dots,a_N$を昇順に並べ替えた数列を$b_0\geq\dots\geq b_N$とする。また、$\displaystyle B_i :=i(i+1) +\sum_{j=i+1}^N \min(b_j,i+1)-\sum_{j=0}^i b_j$とする。すると$\min_{i=0}^N B_i\geq 0$のときErdős–Gallai theoremの条件を満たすこととなる。
例として$N=10, \{a_i\} = \{a_0, 2, 3, 4, 9, 8, 5, 2, 10, 5, 3\}$とし、$a_0$ の値を変えながら計算してみる。色付きは$a_0$を$j-1\to j$に変えたとき値が変化していることを示す。
- $a_0 = 0$
| i | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| $b_i$ | 10 | 9 | 8 | 5 | 5 | 4 | 3 | 3 | 2 | 2 | 0 |
| $\displaystyle\sum_{j=0}^i b_j$ | 10 | 19 | 27 | 32 | 37 | 41 | 44 | 47 | 49 | 51 | 51 |
| $\displaystyle i(i+1)$ | 0 | 2 | 6 | 12 | 20 | 30 | 42 | 56 | 72 | 90 | 110 |
| $\displaystyle\sum_{j=i+1}^N \min(b_j,i+1)$ | 9 | 16 | 19 | 18 | 14 | 10 | 7 | 4 | 2 | 0 | 0 |
| $\displaystyle B_i$ | -1 | -1 | -2 | -2 | -3 | -1 | 5 | 13 | 25 | 39 | 59 |
- $a_0 = 1$
| i | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| $b_i$ | 10 | 9 | 8 | 5 | 5 | 4 | 3 | 3 | 2 | 2 | $\color{red}{1}$ |
| $\displaystyle\sum_{j=0}^i b_j$ | 10 | 19 | 27 | 32 | 37 | 41 | 44 | 47 | 49 | 51 | $\color{red}{52}$ |
| $\displaystyle i(i+1)$ | 0 | 2 | 6 | 12 | 20 | 30 | 42 | 56 | 72 | 90 | 110 |
| $\displaystyle\sum_{j=i+1}^N \min(b_j,i+1)$ | $\color{red}{10}$ | $\color{red}{17}$ | $\color{red}{20}$ | $\color{red}{19}$ | $\color{red}{15}$ | $\color{red}{11}$ | $\color{red}{8}$ | $\color{red}{5}$ | $\color{red}{3}$ | $\color{red}{1}$ | 0 |
| $\displaystyle B_i$ | $\color{red}{0}$ | $\color{red}{0}$ | $\color{red}{-1}$ | $\color{red}{-1}$ | $\color{red}{-2}$ | $\color{red}{0}$ | $\color{red}{6}$ | $\color{red}{14}$ | $\color{red}{26}$ | $\color{red}{40}$ | $\color{blue}{58}$ |
- $a_0 = 2$
| i | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| $b_i$ | 10 | 9 | 8 | 5 | 5 | 4 | 3 | 3 | 2 | 2 | $\color{red}{2}$ |
| $\displaystyle\sum_{j=0}^i b_j$ | 10 | 19 | 27 | 32 | 37 | 41 | 44 | 47 | 49 | 51 | $\color{red}{53}$ |
| $\displaystyle i(i+1)$ | 0 | 2 | 6 | 12 | 20 | 30 | 42 | 56 | 72 | 90 | 110 |
| $\displaystyle\sum_{j=i+1}^N \min(b_j,i+1)$ | 9 | $\color{red}{18}$ | $\color{red}{21}$ | $\color{red}{20}$ | $\color{red}{16}$ | $\color{red}{12}$ | $\color{red}{9}$ | $\color{red}{6}$ | $\color{red}{4}$ | $\color{red}{2}$ | 0 |
| $\displaystyle B_i$ | 0 | $\color{red}{1}$ | $\color{red}{0}$ | $\color{red}{0}$ | $\color{red}{-1}$ | $\color{red}{1}$ | $\color{red}{7}$ | $\color{red}{15}$ | $\color{red}{27}$ | $\color{red}{41}$ | $\color{blue}{57}$ |
- $a_0 = 3$
| i | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| $b_i$ | 10 | 9 | 8 | 5 | 5 | 4 | 3 | 3 | $\color{red}{3}$ | 2 | 2 |
| $\displaystyle\sum_{j=0}^i b_j$ | 10 | 19 | 27 | 32 | 37 | 41 | 44 | 47 | $\color{red}{50}$ | $\color{red}{52}$ | $\color{red}{54}$ |
| $\displaystyle i(i+1)$ | 0 | 2 | 6 | 12 | 20 | 30 | 42 | 56 | 72 | 90 | 110 |
| $\displaystyle\sum_{j=i+1}^N \min(b_j,i+1)$ | 9 | 18 | $\color{red}{22}$ | $\color{red}{21}$ | $\color{red}{17}$ | $\color{red}{13}$ | $\color{red}{10}$ | $\color{red}{7}$ | 4 | 2 | 0 |
| $\displaystyle B_i$ | 0 | 1 | $\color{red}{1}$ | $\color{red}{1}$ | $\color{red}{0}$ | $\color{red}{2}$ | $\color{red}{8}$ | $\color{red}{16}$ | $\color{blue}{26}$ | $\color{blue}{40}$ | $\color{blue}{56}$ |
- $a_0 = 4$
| i | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| $b_i$ | 10 | 9 | 8 | 5 | 5 | 4 | $\color{red}{4}$ | 3 | 3 | 2 | 2 |
| $\displaystyle\sum_{j=0}^i b_j$ | 10 | 19 | 27 | 32 | 37 | 41 | $\color{red}{45}$ | $\color{red}{48}$ | $\color{red}{51}$ | $\color{red}{53}$ | $\color{red}{55}$ |
| $\displaystyle i(i+1)$ | 0 | 2 | 6 | 12 | 20 | 30 | 42 | 56 | 72 | 90 | 110 |
| $\displaystyle\sum_{j=i+1}^N \min(b_j,i+1)$ | 9 | 18 | 22 | $\color{red}{22}$ | $\color{red}{18}$ | $\color{red}{14}$ | 10 | 7 | 4 | 2 | 0 |
| $\displaystyle B_i$ | 0 | 1 | 1 | $\color{red}{2}$ | $\color{red}{1}$ | $\color{red}{3}$ | $\color{blue}{7}$ | $\color{blue}{15}$ | $\color{blue}{25}$ | $\color{blue}{39}$ | $\color{blue}{55}$ |
- $a_0 = 5$
| i | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| $b_i$ | 10 | 9 | 8 | 5 | 5 | $\color{red}{5}$ | 4 | 3 | 3 | 2 | 2 |
| $\displaystyle\sum_{j=0}^i b_j$ | 10 | 19 | 27 | 32 | 37 | $\color{red}{42}$ | $\color{red}{46}$ | $\color{red}{49}$ | $\color{red}{52}$ | $\color{red}{54}$ | $\color{red}{56}$ |
| $\displaystyle i(i+1)$ | 0 | 2 | 6 | 12 | 20 | 30 | 42 | 56 | 72 | 90 | 110 |
| $\displaystyle\sum_{j=i+1}^N \min(b_j,i+1)$ | 9 | 18 | 22 | 22 | $\color{red}{19}$ | 14 | 10 | 7 | 4 | 2 | 0 |
| $\displaystyle B_i$ | 0 | 1 | 1 | 2 | $\color{red}{2}$ | $\color{blue}{2}$ | $\color{blue}{6}$ | $\color{blue}{14}$ | $\color{blue}{24}$ | $\color{blue}{38}$ | $\color{blue}{54}$ |
- $a_0 = 6$
| i | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| $b_i$ | 10 | 9 | 8 | $\color{red}{6}$ | 5 | 5 | 4 | 3 | 3 | 2 | 2 |
| $\displaystyle\sum_{j=0}^i b_j$ | 10 | 19 | 27 | $\color{red}{33}$ | $\color{red}{38}$ | $\color{red}{43}$ | $\color{red}{47}$ | $\color{red}{50}$ | $\color{red}{53}$ | $\color{red}{55}$ | $\color{red}{57}$ |
| $\displaystyle i(i+1)$ | 0 | 2 | 6 | 12 | 20 | 30 | 42 | 56 | 72 | 90 | 110 |
| $\displaystyle\sum_{j=i+1}^N \min(b_j,i+1)$ | 9 | 18 | 22 | 22 | 19 | 14 | 10 | 7 | 4 | 2 | 0 |
| $\displaystyle B_i$ | 0 | 1 | 1 | $\color{blue}{1}$ | $\color{blue}{1}$ | $\color{blue}{1}$ | $\color{blue}{5}$ | $\color{blue}{13}$ | $\color{blue}{23}$ | $\color{blue}{37}$ | $\color{blue}{53}$ |
ここから、$a_0$を$j-1\to j$に変えたとき、 $pos =\displaystyle\sum_{k=1}^N \mbox{1}\hspace{-0.25em}\mbox{l}_{\{a_k\geq j\}}$とすると
- $i=pos,\dots,n$では$B_i$が1減る
- $i=j-1,\dots,pos-1$では$B_i$が1増える
ことがわかる。$B$に上記の変更をした後、$\displaystyle\min_{i=0}^N B_i\geq 0$かどうかを確かめる、という操作を$j=1,\dots,N$について行うことで解ける。これは区間add-区間minのSegment Treeを使うと$\mathcal O(N\log N)$。
解答例
#include <bits/stdc++.h> using LL = long long; constexpr LL LINF = 334ll << 53; constexpr int INF = 15 << 26; namespace Problem { using namespace std; struct RangeMin { using T = long long; T operator()(const T &a, const T &b) { return min(a, b); } static constexpr T id() { return T(LINF); }; }; struct RangeAdd { using T = long long; T operator()(const T &a, const T &b) { return a + b; }; static constexpr T id() { return T(0); }; }; struct Min_Add { using V = RangeMin; using A = RangeAdd; V::T operator()(const V::T &a, const A::T &b, const int &h) { return a + b; }; }; template <typename M> struct LazySegmentTree { using Val = typename M::V::T; using Del = typename M::A::T; typename M::V calc; typename M::A composite; M act; vector<Val> val; vector<Del> delay; int n, height; explicit LazySegmentTree(int size) : n(1), height(0) { while (n < size) { n <<= 1; ++height; } val = vector<Val>(n << 1, calc.id()); delay = vector<Del>(n << 1, composite.id()); } explicit LazySegmentTree(vector<Val> &v) : n(1), height(0) { while (n < v.size()) { n <<= 1; ++height; } val = vector<Val>(n << 1, calc.id()); delay = vector<Del>(n << 1, composite.id()); copy(v.begin(), v.end(), val.begin() + n); for (int i = n - 1; i > 0; --i) { val[i] = calc(val[i * 2], val[i * 2 + 1]); } } void propagate(int k, int h) { val[k] = act(val[k], delay[k], h); if (h > 0) { delay[2 * k] = composite(delay[2 * k], delay[k]); delay[2 * k + 1] = composite(delay[2 * k + 1], delay[k]); } delay[k] = composite.id(); } inline void update(int l, int r, Del v) { update(l, r, v, 1, 0, n, height); } void update(int l, int r, Del v, int k, int a, int b, int h) { if (l <= a && b <= r) { delay[k] = composite(delay[k], v); propagate(k, h); return; } propagate(k, h); if (b <= l || r <= a) { return; } // update(l, r, v, 2 * k, a, (a + b) / 2, h - 1); update(l, r, v, 2 * k + 1, (a + b) / 2, b, h - 1); val[k] = calc(val[2 * k], val[2 * k + 1]); } inline Val query(int l, int r) { return query(l, r, 1, 0, n, height); } Val query(int l, int r, int k, int a, int b, int h) { if (b <= l || r <= a) { return calc.id(); } // propagate(k, h); if (l <= a && b <= r) { return val[k]; } Val vall = query(l, r, 2 * k, a, (a + b) / 2, h - 1); Val valr = query(l, r, 2 * k + 1, (a + b) / 2, b, h - 1); return calc(vall, valr); } }; class Solver { public: int n; vector<LL> a, r, sum; Solver(LL n) : n(n), a(n + 1), r(n + 2), sum(n + 1){}; void solve() { for (int i = 0; i < n; ++i) { cin >> a[i + 1]; } sort(a.rbegin(), a.rend()); for (int i = 0; i <= n; ++i) { r[a[i]] = i+1; } for (int i = n; i >= 0; --i) { r[i] = max(r[i], r[i + 1]); } sum[0] = a[0]; for (int i = 1; i <= n; ++i) { sum[i] = sum[i - 1] + a[i]; } //a[0]=0のときのBを計算 vector<LL> v(n + 1); for (LL k = 0; k <= n; ++k) { LL lhs = sum[k]; LL rhs = k * (k + 1); if (r[k + 1] >= k + 1) { rhs += (r[k + 1] - (k + 1) ) * (k + 1); rhs += sum[n] - sum[r[k + 1]-1]; } else { rhs += sum[n] - sum[k]; } v[k] = rhs - lhs; } LazySegmentTree<Min_Add> seg(v); vector<LL> ans; if (sum[n]%2==0 && seg.query(0, n + 1) >= 0) { ans.push_back(0); } //a[0]を増やしながらErdős–Gallai theoremの条件を確認 for (int i = 1; i <= n; ++i) { LL pos = r[i] ; r[i]++; seg.update(pos, n + 1, -1); seg.update(i - 1, pos, 1); if ((sum[n] + i) % 2 == 0 && seg.query(0, n + 1) >= 0) { ans.push_back(i); } } //output if (ans.size() == 0) { cout << -1 << "\n"; return; } for (int i = 0; i < (int)ans.size(); ++i) { cout << ans[i] << ' '; } cout << endl; } }; } // namespace Problem int main() { std::cin.tie(0); std::ios_base::sync_with_stdio(false); long long n = 0; std::cin >> n; Problem::Solver sol(n); sol.solve(); return 0; }
note
wikipediaの定理の式が1-indexedで、それを見ながら0-indexedで解いていたら大混乱した 。しかもmain testで落ちた。横着しないで0-indexedの式に書き直すべきだった。
AtCoder ARC070 D No Need
問題概要
枚のカードがあり，カード
には正の整数
が書かれている．以下の条件を満たす
の個数を求めよ．
条件
カード$i$以外の$N-1$枚のカードの部分集合で，書かれている数の総和Sumが$K-a_i \leq Sum\leq K-{1}$となる部分集合は存在しない．
制約
- $1\leq N,K\leq 5000$
- $1 \leq a_i \leq 10^9$
解法
集合$A\subseteq\{1,\ldots,n\}$に対して
$S_{{A}}$:$A$の部分集合のカードに書かれた総和として作ることができるK未満の数の集合
とする．$S_{\{1,\ldots,n\}-\{i\}}$をすべての$i$について求めれば，判定は$O(NK)$で可能．
一般に$S_{A}$と$S_{B}$から$S_{A\cup B}$を求める操作は$O(K^2)$の計算量のDPとなるが，$S_{A}$と$S_{\{j\}}$から$S_{A\cup\{j\}}$を求めることは$O(K)$で可能である． この性質から，$S_{\{1,\ldots,n\}-\{i\}}$たちは平方分割により$O(N^{1.5}K)$で計算できる．bitsetのおかげでいわゆる$"O(N^{1.5}K/64)"$となり間に合う．
note
Segment Tree + NNT で$O(N\log N K\log K)$にもできそう(最初それで解こうとして600点…??という気持ちになった)．
codeforces Marathon Match 2 Colored Path
問題
http://codeforces.com/contest/1014/problem/A2
無限に続くマス目の列があり，0,1,2,..と番号がつけられている．それぞれのマスには'A','B',...,'J'の10種類の色のうち1つがつけられている． 与えられた初期状態では，プレイヤーは10色の塗料のボトルを合計10本持っていて(1本づつとは限らない)いる．
初期状態ではマス0,1,...9に10匹のカメレオンがいて，これらを塗料のボトルを投げつける操作を10000回行うことでできるだけ番号の大きいマスへ動かしたい．
操作は以下のように行われる． * プレイヤーは手持ちのボトルのうち1本と対象のカメレオン1匹を指定する． * 対象となったカメレオンは，現在自分がいるマスの番号より大きく，投げられたボトルと同じ色を持ち，他のカメレオンに占有されていないようなマスのうち最も番号の小さいマスへ移動する． * ボトルを投げたプレイヤーは，バッグから1本のボトルを取りだし，手持ちのボトルを再び10本とする．バッグから取り出されるボトルの順序はあらかじめテストケースごとに指定されている．
各テストケースの得点は，10000回の操作終了後にカメレオンがいるマスの番号のうち最も小さい数である．
成績
??? 暫定8位
解法・参加記
7/25
全探索を最も愚直に行うと, 各操作で投げるボトルと対象のカメレオンが10種類ずつあるので通り．重複を取り除けても到底無理．
ここで手持ちの10本のボトルの組み合わせは92378 通りに限られるので，これらに0,1,...,92377と番号をつける．
(カメレオンの位置，手持ちのボトル)を「状態」として，状態の良さを測る評価関数があると仮定した上で次のDPを考える．
回の操作後に手持ちのボトルの組が
の時の評価関数の最大値
このDPは状態数が , 遷移の回数が
となり，ちょっと可能そうに思えてくる．
しかしこのDPでも間に合わないので，i回目→i+1回目の遷移を92738通りのすべてから行うのではなく，評価関数の値の大きいwidth個(width~100)の
からのみ遷移を行うというビームサーチ風の手法にした．
次に評価関数について考えてみる．もしDPで全遷移を計算できたとして，その時最適解が求まるためには，「i回目の操作後に手持ちのボトルの組が同じになるような投げ方について，評価関数の値が
ならば，(i+1)回目以降を最適に投げた時の最終的なスコアは
を経る操作の方が大きい」という性質が成り立っていてほしい．これにできるだけ近そうな評価関数を作ることが目標となる．
とりあえず評価関数を
static int eval(vector<int> position){ int ret = 0; sort(position.begin(), position.end()); for (int i = 0; i < 10; ++i) ret += position[i] * (10-i); return ret; }
として， 37795.24点(width = 60)．
7/26
若干の高速化．それと評価関数を
static int eval(vector<int> position){ int ret = 0; sort(position.begin(), position.end()); for (int i = 0; i < 10; ++i) ret += position[i] * (19-i); return ret; }
に変更．39355.21点(width = 72)．
7/27
実は評価関数でpositionをソートしているのがボトルネックになっていた．だが単純にpositionの総和を評価とすると，カメレオン集団が前後に分断されて一部のみが前進してしまいうまくいかない．そこで，総和 - 標準偏差*0.2 を評価関数として分断を防いだ．この変更によりスコアは少し下がったが，高速化の効果が上回り，39754.64 点(width = 101).
7/28
最後の20回に限り，評価関数に(最も後ろのカメレオンの位置)×10を足して目的関数に近づけた．
ビームサーチの多様性が低そう(あるターンで評価の高かったものばかり次に引き継がれる)だったが改善方法が思いつかず．評価関数に0~9の乱数を加算してみたらスコアが上がったが，探索方法の改善にうまく繋がらない．
これらの変更によって40000点を突破した．40091.95 点(width=97)．
note
焼きなましできるのだろうか？どこかの操作を変えた時に更新がO(S)より速くなる気がしないので，10000回の操作を200*50などに分割して，前のブロックからブロックごとに焼きなまして行けばいいのか？
AtCoder codeFlyer2018 本戦 H 三角形と格子点
問題
三角形ABCの3頂点を
- A
- B
- C
となるようにとる．このようにして作られる個の三角形について三角形の内部に存在する格子点の個数の総和を求めよ．(mod
)
制約
残りの制約はA,B,Cの位置関係について
x_a < x_c< x_b, y_a < y_b < y_c
が成り立つことを示している．
解法
ピックの定理から各三角形について
面積 = (内部の格子点数) + (周上の格子点数)/2-1
が成り立つ．これを個の三角形について和をとると
(面積の和) =(答) + (周上の格子点数の和)/2 - .
よって面積の和と周上の格子点数の和が分かれば解ける．
面積の和
A,B,Cの座標を
- A
- B
- C
とすると面積の和は
.
面積の平均がになることからもこれがわかる．
周上の格子点数の和
ある三角形ABCについて，
となる三角形は
個ある．
- そのような三角形ABCの辺AB上(点Bを除く)にある格子点数は
.
よって個の三角形について辺AB上(点Bを除く)にある格子点数の総和は以下のように求められる．
ただしとする．
また，オイラー関数は
を満たすことに注意．
\begin{align} &\sum_{i = x_2-x_1 -w +1}^{x_2-x_1 +w -1}\sum_{j = y_2-y_1 -h +1}^{y_2-y_1 +h -1} hw(w-|i-x_2+x_1|)(h-|j-y_2+y_1|){\text gcd}(i,j)\\ =&hw\sum_{i = x_2-x_1 -w +1}^{x_2-x_1 +w -1}(w-|i-x_2+x_1|)\sum_{k|i}\sum_{p\in S_{ik}}(h-|p-y_2+y_1|)k\\ =&hw\sum_{i = x_2-x_1 -w +1}^{x_2-x_1 +w -1}(w-|i-x_2+x_1|)\sum_{k|i}\sum_{p\in S_{ik}}(h-|p-y_2+y_1|)\sum_{l|k}\phi(l) \\ =&hw\sum_{i = x_2-x_1 -w +1}^{x_2-x_1 +w -1}(w-|i-x_2+x_1|)\sum_{k|i}\sum_{l|k}\phi(l)\sum_{p\in \cup_{l|m}S_{im}}(h-|p-y_2+y_1|) \\ =&hw\sum_{i = x_2-x_1 -w +1}^{x_2-x_1 +w -1}(w-|i-x_2+x_1|)\sum_{k|i}\sum_{l|k}\phi(l)\sum_{l|p}(h-|p-y_2+y_1|) \\ =&hw\sum_{i = x_2-x_1 -w +1}^{x_2-x_1 +w -1}(w-|i-x_2+x_1|)\sum_{l|i}\phi(l)\sum_{l|p}(h-|p-y_2+y_1|) \\ \end{align}
は等差数列の和なので, それぞれの
についてO(1)時間で求められる．
の計算のために，まず
の範囲について，区間ふるいによる素数判定を行うと同時に素因数分解しておく．すると任意の
について
の約数pとオイラー関数の組(
,
)を約数の個数について線形時間で求められる．
同じ計算を辺BC, CAについても行えばよい．
code
#include <bits/stdc++.h> using LL = long long; constexpr LL LINF = 334ll << 53; constexpr int INF = 15 << 26; constexpr LL MOD = 1E9 + 7; namespace Problem { using namespace std; template <long long Mod = MOD> class Modint { public: long long v; Modint(const Modint &x) { v = x.v; } explicit Modint(int x) { if (x < 0 or Mod <= x) x %= Mod; if (x < 0) x += Mod; v = x; } explicit Modint(long long x) { if (x < 0 or Mod <= x) x %= Mod; if (x < 0) x += Mod; v = x; } explicit Modint(unsigned x) { if (Mod <= x) x %= Mod; v = x; } explicit Modint(unsigned long long x) { if (Mod <= x) x %= Mod; v = x; } template <typename T> explicit Modint(T x) { if (x < 0 or Mod <= x) x %= Mod; v = x; } Modint() : v(0) {} long long get() const { return v; } Modint mpow(Modint &n) const { return mpow(n.v); } Modint mpow(long long n) const { long long i = 1, p = v; Modint ret(1); while (i <= n) { if (i & n) ret *= p; i = (i << 1); p = (p * p) % Mod; } return ret; } Modint operator-() const { return (v ? Modint(Mod - v) : Modint(0)); } explicit operator int() const { return v; } explicit operator long long() const { return v; } Modint &operator+=(const Modint &a) { v = (v + a.v); if (v >= Mod) v -= Mod; return *this; } Modint &operator-=(const Modint &a) { v = (v - a.v); if (v < 0) v += Mod; return *this; } Modint &operator*=(const Modint &a) { v = (v * a.v) % Mod; return *this; } Modint &operator/=(const Modint &a) { v = (v * a.mpow(Mod - 2).v) % Mod; return *this; } template <class T> Modint &operator+=(const T &a) { v = (v + Modint(a).v); if (v >= Mod) v -= Mod; return *this; } template <class T> Modint &operator-=(const T &a) { v = (v - Modint(a).v); if (v < 0) v += Mod; return *this; } template <class T> Modint &operator*=(const T &a) { v = (v * Modint(a).v) % Mod; return *this; } template <class T> Modint &operator/=(const T &a) { v = (v * Modint(a).mpow(Mod - 2).v) % Mod; return *this; } friend Modint operator+(const Modint &a, const Modint &b) { return Modint(a) += b; } friend Modint operator-(const Modint &a, const Modint &b) { return Modint(a) -= b; } friend Modint operator*(const Modint &a, const Modint &b) { return Modint(a) *= b; } friend Modint operator/(const Modint &a, const Modint &b) { return Modint(a) /= b; } template <class T, class U> friend Modint operator+(const T &a, const U &b) { return Modint(a) += Modint(b); } template <class T, class U> friend Modint operator-(const T &a, const U &b) { return Modint(a) -= Modint(b); } template <class T, class U> friend Modint operator*(const T &a, const U &b) { return Modint(a) *= Modint(b); } template <class T, class U> friend Modint operator/(const T &a, const U &b) { return Modint(a) /= Modint(b); } }; template <long long M> ostream &operator<<(ostream &os, const Modint<M> m) { return os << m.v; } class SegmentSieve { long long lower, upper, range; vector<long long> primes; // prime number in[2,sqrt(b)] vector<vector<pair<long long, int>>> prime_factor; vector<char> is_prime_seg; void sieve() { vector<char> isprime(range + 1, 1); vector<long long> table(upper - lower); iota(table.begin(), table.end(), lower); isprime[0] = isprime[1] = 0; for (long long i = 2; i <= range; i++) { if (isprime[i]) { primes.push_back(i); for (long long j = i * i; j <= range; j += i) isprime[j] = 0; for (long long j = lower - lower % i; j < upper; j += i) { if (j >= lower) { is_prime_seg[j - lower] = false; int ord = 0; while (table[j - lower] % i == 0) { table[j - lower] /= i; ord++; } prime_factor[j - lower].emplace_back(i, ord); } } } } for (int i = 0; i < (int)prime_factor.size(); ++i) { if (table[i] != 1) prime_factor[i].emplace_back(table[i], 1); } } public: SegmentSieve(long long a, long long b) : lower(a), upper(b), prime_factor(b - a), is_prime_seg(b - a, 1) { //[a,b) range = max(b - a + 1, (long long)sqrt(b) + 1); sieve(); }; bool is_prime(long long x) { // assert(lower <= x && x < upper); return is_prime_seg[x - lower]; } vector<pair<long long, long long>> divisor_euler(long long x) { // assert(lower <= x && x < upper); vector<pair<long long, long long>> ret = {{1ll, 1}}; for (const auto p : prime_factor[x - lower]) { int sz = (int)ret.size(); for (int i = 0; i < sz; ++i) { long long divisor = ret[i].first; long long euler_ = ret[i].second; for (int d = 1; d <= p.second; ++d) { divisor *= p.first; euler_ *= (d == 1 ? (p.first - 1) : p.first); ret.emplace_back(divisor, euler_); } } } return ret; } }; class Solver { public: LL w, h; vector<LL> x, y; Solver() : x(3), y(3){}; void solve() { for (int i = 0; i < 3; ++i) { cin >> x[i] >> y[i]; } cin >> w >> h; Modint<> area(0), on_edge(0); area = Modint<>(h).mpow(3) * Modint<>(w).mpow(3) / 2 * (Modint<>(x[1] - x[0]) * Modint<>(y[2] - y[0]) - Modint<>(x[2] - x[0]) * Modint<>(y[1] - y[0])); on_edge += calc(x[0], y[0], x[1], y[1]); on_edge += calc(x[0], y[0], x[2], y[2]); on_edge += calc(x[1], y[1], x[2], y[2]); cout << area + Modint<>(h).mpow(3) * Modint<>(w).mpow(3) - on_edge / 2 << endl; } Modint<> calc_sum(LL k, LL center, LL width) { Modint<> ret(0); //等差数列の和の部分 long long r = center - center % k; long long l = (center - width + k - 1) / k * k; ret += Modint<>((r - l) / k + 1) * Modint<>(h - abs(center - r) + h - abs(center - l)) / 2; long long r2 = center + width - (center + width) % k; long long l2 = (center + 1 + k - 1) / k * k; ret += Modint<>((r2 - l2) / k + 1) * Modint<>(h - abs(center - r2) + h - abs(center - l2)) / 2; return ret; } Modint<> calc(LL a, LL b, LL c, LL d) { if (a > c) swap(a, c); if (b > d) swap(b, d); SegmentSieve segs(c - a - w + 1, c - a + w); Modint<> ret(0); /*for (LL i = c - a - w + 1; i < c - a + w; ++i) { Modint<> tmp(0); for (LL j = d - b - h + 1; j < d - b + h; ++j) { tmp += Modint<>(__gcd(i, j)) * (h - abs(j - d + b)); } ret += Modint<>(h * w) * tmp * (w - abs(i - c + a)); }*/ for (LL i = c - a - w + 1; i < c - a + w; ++i) { Modint<> tmp(0); for (auto p : segs.divisor_euler(i)) { auto div = p.first; auto euler = p.second; tmp += calc_sum(div, d - b, h) * euler; } ret += Modint<>(h * w) * tmp * (w - abs(i - c + a)); } return ret; } }; } // namespace Problem int main() { std::cin.tie(0); std::ios_base::sync_with_stdio(false); Problem::Solver sol; sol.solve(); return 0; }
note
知らなかったこと
- 長さ 𝑂(𝐻) の区間に含まれる整数の約数の個数の和は 𝑂(𝐻 log 𝐻) (解説より)
特にオイラー関数の性質を知らなかったので，自力でACは無理だった．
