SRM 643 Div1 Easy, TheKingsFactorization
本番は素数リスト作ったけどその必要も無かった.
問題
http://community.topcoder.com/stat?c=problem_statement&pm=13594&rd=16086
1018 以下の整数 N が与えられる. N を素因数分解せよ. ただしヒントとして, 昇順 0, 2, 4, 6, ... 番目の素因数が与えられる. (e.g. N = 60 の場合, 素因数 {2, 2, 3, 5} のうち {2, 3} が与えられる.)
方針
ヒントを考慮しない場合, √1018 = 109 までの試し割りが必要となり, 計算量オーバーする.
ヒントの素因数以外で, 106 を超える素因数は高々 1 つである. よって, 106 までを試し割れば良い.
解答
using ll = long long; const int MAX = 1000001; class TheKingsFactorization { public: vector<long long> getVector (long long N, vector<long long> primes) { for (ll p : primes) N /= p; for (int i = 2; i < MAX; i++) { while (N%i == 0) { primes.push_back(i); N /= i; } } if (N > 1) primes.push_back(N); sort(primes.begin(), primes.end()); return primes; } };
SRM 641 Div1 Easy, TrianglesContainOrigin
問題
http://community.topcoder.com/stat?c=problem_statement&pm=13309&rd=16084
2 次元平面上の点が N (< 2500) 個与えられる. これらの点が構成する三角形のうち, 内側に原点を含むものの個数を求めよ.
方針
偏角でソートし, 全ての点 i, j (i < j) について, 条件を満たす点 k (> j) の個数を二分探索で求める.
解答
using ll = long long; const double PI = atan(1)*4; class TrianglesContainOrigin { public: long long count (vector <int> x, vector <int> y) { int N = x.size(); vector<double> args(N); loop (N, i) args[i] = atan2(y[i], x[i]); sort(args.begin(), args.end()); ll ans = 0; loop (N, j) loop (j, i) { if (args[j] > args[i]+PI) continue; auto s = lower_bound(args.begin(), args.end(), args[i] + PI); auto t = lower_bound(args.begin(), args.end(), args[j] + PI); ans += distance(s, t); } return ans; } };
SRM 642 Div1 Med, TaroCutting
Med は難しい...
問題
http://community.topcoder.com/stat?c=problem_statement&pm=13319&rd=16085
木が N 本ある. n 番目の木は高さ height[n] であり, 1 日に add[n] だけ伸びる. また, カッターが M 個ある. m 番目のカッターをある木に対して使うと, 木の高さは device[m] になる.
太郎は 1 日の終わりにカッターで木を切る. 各々のカッターは 1 日に 1 回, 1 つの木に対してのみ使うことが出来る.
カッターを上手く使って, time 日目の終わり時点(カッターを使った後)での木の高さの和を最小にしたい. その時の和を求めよ.
- N, M, time <= 150
- height[n], add[n], device[m] <= 10000
方針
1 つの木を 2 回以上切っても無意味なので, 1 つの木を切るのは 0 回または 1 回である. また, add が大きい木は出来るだけ遅い日に切った方が良い.
add が大きい木を優先的に切るために, 木を add の降順でソートする.
dp(t, i, j) を「(t+1 日目以降で最適な切り方をした時の)t 日目に, i 番目までの木に対して, j 番目までのカッターをうまく使った時の, i 番目までの木の総長」とすると, 次の漸化式が得られる.
- dp(t, i, -1) = dp(t+1, i, M-1)
- カッターを使わない場合は, t+1 日目以降の最適値と一致する
- dp(t, i, j) = 以下の最小値
- dp(t, i-1, j) + height[i] + time * add[i]
- dp(t, i-1, j-1) + device[i] + (time-t) * add[i]
- dp(t, i, j-1)
上の漸化式をメモ化再帰で実装する.
解答
using P = pair<int, int>; const int INF = 1<<29; class TaroCutting { public: int N, M; vector<P> tree; vector<int> device; int time; int getNumber (vector <int> height, vector <int> add, vector <int> device, int time) { N = height.size(); M = device.size(); this->device = device; this->time = time; loop (N, i) tree.emplace_back(add[i], height[i]); sort(tree.begin(), tree.end(), greater<P>()); fill(memo[0][0], memo[0][0]+160*160*160, -1); return rec(1, N-1, M-1); } int memo[160][160][160]; int rec (int t, int i, int j) { if (t > time) return INF; if (i < 0) return 0; if (j < 0) return rec(t+1, i, M-1); if (memo[t][i][j] != -1) return memo[t][i][j]; int ans = rec(t, i-1, j) + time * tree[i].first + tree[i].second; ans = min(ans, rec(t, i-1, j-1) + device[j] + (time-t) * tree[i].first); ans = min(ans, rec(t, i, j-1)); return memo[t][i][j] = ans; } };
考察
本番で漸化式作って証明までするのは不可能だと思う. おとなしく最小費用流で書くべき*1な気がする.
SRM 642 Div1 Easy, WatingForBus
DP の練習に良い問題だとおもった
問題
http://community.topcoder.com/stat?c=problem_statement&pm=13540&rd=16085
同じ停留所に発着するバスが N (< 100) 台あり, 0 から N-1 までの番号が割り当てられている. i 番目のバスが発車してから戻ってくるまでに, time[i] (< 105) の時間がかかる. 時刻 0 に, N 台のうちランダムに選ばれた 1 台が発車する. バスが戻ってくると同時に, また N 台のうちランダムに選ばれた 1 台が発車する. これを繰り返す. 同じバスが続けて選ばれることもある. i 番目のバスが選ばれる確率は prob[i] / 100 で与えられる.
あなたは 時刻 s に停留所に到着した. 次のバスが発車するまでの待ち時間の期待値を出力せよ.
方針
- 時刻 0 から s のそれぞれでバスが発(着)車する確率 P を DP で求める
ある時刻 i について, P[i] = Σ[j = 0..N-1] P[i-time[j]] * prob[j] / 100 - 上で求めた確率を使って, 「時刻 s から s+i-1 の間にバスが発(着)車せず、時刻 s+i にバスが発(着)車する確率」P' を求める
- Σ P' * i が求める期待値
解答
const int M = 200010; class WaitingForBus { public: double whenWillBusArrive (vector <int> time, vector <int> prob, int s) { int N = time.size(); vector<double> P(M, 0); P[0] = 1; loop (M-1, i) { loop (N, j) { if (i+1 >= time[j]) P[i+1] += P[i+1-time[j]] * prob[j] / 100; } } double ans = 0; loop (100001, i) { loop (N, j) { if (i < time[j] && s+i >= time[j]) ans += i * P[s+i-time[j]] * prob[j] / 100; } } return ans; } };
社会人 2 年目が競技プログラミングを始めた 1 年を振り返る
この記事は Competitive Programming Advent Calendar Div2014 12 日目です。
競技プログラミングを始めて 1 年がたちました。この一年をざっと振り返ってみたいと思います。
【1 月】PrintScreen / Aizu Online Judge
新卒 SIer 1 年目だった僕は、一日 8 時間のうち 7 時間待機、1 時間でプリントスクリーンを Excel に貼り付ける業務を任され、時間を持て余していた。そこで Aizu Online Judge と出会い、ハマった。AOJ サブミット用のスクリプトを作ったりした。
小難しいデータ構造だとか理論を知らなくても、発想の転換で簡単に解けたりするのはおもしろいなーと。例えば 1028: Ideal Coin Payment and Change なんかは、ただの全数探索なんだけど発想の転換の妙を初めて感じた問題で、答えを見て震えたのを覚えている。
【2 月】僻地出張 / TopCoder
僻地に長期出張、ホテル暮らしとなり、かつ相変わらず業務の大半は待機ということで、自由に使える時間が一日 15 時間くらいある生活。ここらへんで TopCoder を始めてみる。初回の戦績としては、90 分かけてなんとか easy を解き、medium を2つ撃墜というまずまずの結果だった。
その後は参加する度にレートを下げて、灰コーダーになる。
【4 月】Web 系
業務の待機時間に飽き飽きし、転職。楽しい Web 系が始まった。自己紹介で競技プログラマーアピールした。灰色なんて言えない。毎日練習に励む。緑組に復帰。
あと AtCoder が毎週あってすごく有難かった。
【6 月】青コーダー
部署内の週次報告回で TopCoder の戦績を発表したりしていた。練習のかいあってか、ある時突然 Div.2 が簡単に感じるようになる。SRM 623 で一気にレート上げて青コーダーになった。
青コーダーになってからはグラフ理論の素振りをよくやっていた気がする。蟻本を買う。
【9 月】Rating 1499
9 戦連続 1 完でジリジリレートを上げて、あと 1 ポイントで黄色というところまできた。初戦からぴったり半年で黄色コーダーや!と思いきや、その後まさかの 2 連続 0 完で大きく下げる。
初めての幾何問題 に手も足も出ず、悔しい思いをする。
【10 月】
社内でライブ SRM とか、社外向けのエンジニアイベントで競技プログラミングについて喋ったりした。
ライブ SRM は、、、スピード感全然無いしインタラクティブじゃ無いし緊張して頭回らないしやるべきじゃないと思った。そもそもレギュレーション違反である。
【11 月】黄コーダー
撃墜回で上げることが多い気がする。correct 率 15 % の SRM 638 Easy を pass して黄コーダーになった。めっちゃはしゃいだ。
レート推移
競技プログラミングは楽しい!
他人と何かで競って、自分の立ち位置の上がったり下がったりが目に見えるのは楽しい。受験とか音ゲーに似ていると思う。
競技プログラミングの良い所は楽しいだけじゃなく実務にも活かせることだと思っていて、競プロを始めてから複雑な処理がシンプルにサクッと書けるようになったし、計算量の勘が効くようになったおかげで、システム設計の方でも注力すべきボトルネックが見えてくるようになったり、実感するところは多い。
学生時代に競プロやってなかったのが惜しいですが、これから競プロの交友関係も作っていけたらなーと思ってます。
SRM 639 Div1 Easy, AliceGame
問題
http://community.topcoder.com/stat?c=problem_statement&pm=13490&rd=16082
Alice, Kirito 2人がゲームをする. ゲームはターン1から始まり有限ターンからなる. 各ターンでそれぞれ勝者が決まり, i ターン目で勝利すると2*i-1ポイントもらえる. long long int x, y が与えられる. ゲーム終了時に Alice が x ポイント, Kirito が y ポイントとなるようなパターンは存在するか. 存在するならば Alice が勝利する必要がある最小ターン数を出力せよ. 存在しないならば -1 を出力せよ.
0 <= x, y <= 1012
方針
ターン数 T は √(x+y) で得られる.
計 i ターンの勝利で x ポイント稼げる必要十分条件:
- x が 1 以上 2*T-1 以下の相異なる i 個の奇数の和で表せる
- ⇔ (x+i)/2 が 1 以上 T 以下の相異なる i 個の整数の和で表せる
この条件を満たすかどうかを i = 1 .. T で確かめれば良い.
解答
class AliceGame { public: long long findMinimumValue(long long x, long long y) { ll turns = sqrt(x+y)+0.01; if (turns*turns != x+y) return -1; for (ll i = 0; i <= turns; i++) { ll min = i*(i+1)/2, max = (2*turns-i+1)*i/2; if ((x+i) % 2 == 0 && (x+i)/2 >= min && (x+i)/2 <= max) return i; } return -1; } };
雑感
大きい奇数から順に足していく方法だと, 2 が残る場合のコーナーケースがあり, ここで大勢落ちていたっぽい. 2 が残る場合だけ除外するとそれでも通るが, 証明がわからない.
Haskellでエラトステネスのふるい (Data.Set編)
最近、Project Euler の問題をぽちぽち解いていってます。どうせやるなら慣れていない言語でやろう!ということで Python とか Haskell を触ってるのですが、イミュータブルなデータ構造が基本の Haskell でアレコレやるのは本当に難しい。
例えば Problem 10 - Project Euler を考えてみます。
The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17.
Find the sum of all the primes below two million.
2*106 未満の素数の和を求めよという問題。エラトステネスのふるい で素数リストを作って足せば良い問題です。
C++
C++ で実装するとこんな感じになります。
vector<bool> isPrime(N, true); isPrime[0] = isPrime[1] = false; for (int i = 2; i*i < N; i++) { if (!isPrime[i]) continue; for (int j = i*i; j < N; j += i) isPrime[j] = false; }
比較のために、単純な試し割り法も載せてみます。
vector<bool> isPrime(N, true); isPrime[0] = isPrime[1] = false; for (int i = 2; i < N; i++) { for (int j = 2; j*j <= i; j++) { if (i%j == 0) isPrime[i] = false; } }
試し割りの場合は二重ループが全て実行されて計算量 O(N√N) であるのに対して、エラトステネスのふるいのループ回数は大体 N/2 + N/3 + N/5 + N/7 + ... で、素数の逆数和 = loglogN より、O(NloglogN) くらいになります。
ポイントは、各値に約数があるかどうかを調べるのではなく、各値の倍数を素数リストから除外 していること。
Haskell
Haskell でナイーブな試し割り法を実装するとすごくシンプルになります*1。
primes :: [Integer] -> [Integer] primes [] = [] primes (x:xs) = x : primes [y | y <- xs, mod y x /= 0]
一方でエラトステネスの実装はシンプルにはいかないようです。Haskell ではイミュータブルなデータ構造が基本なので、「素数リストから除外」を定数時間で実行させるのが難しい。
例えば C++ の例のように、素数リストを Bool の List で表して除外処理を実装した場合、除外位置より手前の全要素をコピーすることになり、除外処理の計算量は O(N) 、全体の計算量は O(N2 loglogN) となってしまいます。これではナイーブな試し割りより性能が悪いです。
そこで、Data.Set を使ってみます。これはいわゆる平衡木で、insert と delete が O(logN) で実行できるデータ構造です。コピーして再構成するのに O(N) かかりそうなものだけど、再構成するのは対象ノードからルートまでのパスだけで良いので O(logN) で済むらしい。haskell - How is insert O(log(n)) in Data.Set? - Stack Overflow
こちらの記事 では Data.Set を使って O(N*log2 N) のふるいが実装されています。下のコードはこれを少し改良して、C++で書いた上述のアルゴリズムに近づけたものです*2。
seive :: Integer -> Set Integer -> Set Integer seive p xs | p*p > findMax xs = xs | otherwise = seive nextP (xs \\ mults) where mults = fromAscList [p*p, (p+1)*p .. findMax xs] nextP = fromJust $ lookupGT p xs primes :: Integer -> [Integer] primes m = toList $ seive 2 $ fromAscList [2..m]
xs (素数候補リスト) から mults (素数の倍数) を除外する処理を繰り返し実行します。一つの値の除外処理に O(logN) かかるので、計算量は O(N*logN*loglogN) となります。Problem 10 - Project Euler で求められる 2*106 の場合、5秒くらいで計算できました。
参考
試し割り: 誇大宣伝ではないエラトステネスの篩 - HaHaHa!(old) - haskell
優先順位つきキュー O(nlog2 n): Haskellでエラトステネスの篩 - 簡潔なQ
Data.Set O(nlog2 n): Tech Tips: Haskellで有名アルゴリズム実装(4)
ミュータブルな実装: Haskellでエラトステネスの篩(STUArray) - 似非学問的な手記