1~nまでの逆元をO(n)で求める方法
参考にしたのは以下の@rng_58さんによるツイートです
導出が自明ではないのでメモ（少なくとも僕にとっては）
ideone.com/SOygoW これで 1 から N までの逆元が O(N) でもとまります
— rng_2さん (@rng_58) 2013年3月15日@komaki__ i * (MOD/i) = MOD - MOD%i を変形すると出てきます
— rng_2さん (@rng_58) 2013年3月15日逆元列挙のコード
以下のコードで1~nまでの逆元(mod m)をO(n)で求められます
// 1~nの逆元を求める(mod m) vector<ll> list_mod_inverse(ll n, ll m) { vector<ll> inv(n + 1); inv[1] = 1; for (int i = 2; i <= n; ++i) inv[i] = inv[m % i] * (m - m / i) % m; return inv; }
導出
説明では、mを法としています。また、式中の/除算は切り捨て除算です。
上記のコードではとして、逆元を求めています。この式を導出します。
mはと表すことができます。この式中のkは(m / i)に等しくなります。よって、
となります。この式を変形をしていきます。
(※以下の式はです。）
より、
を掛ける
を掛ける
負の数を消す
以上です
AOJ 580 Fish (おそらく非想定解？)
AOJで他の解法投げると見れなくなると思うので，とりあえずコード貼っておきます
解法
直方体とその領域で重なっている数を持っておく．
ソースコード
struct Rect { int v[3][2]; ll calc_area() const { ll res = 1; rep(d, 3) res *= v[d][1] - v[d][0]; return res; } }; typedef pair<Rect, int> Area; // (rect, num of fishes) bool intersect(const Rect& a, const Rect& b, Rect& r) { rep(d, 3) { r.v[d][0] = max(a.v[d][0], b.v[d][0]); r.v[d][1] = min(a.v[d][1], b.v[d][1]); if (r.v[d][0] >= r.v[d][1]) return false; } return true; } int main() { int n, k; cin >> n >> k; const int inf = ten(8); Rect inf_rect; rep(d, 3) { inf_rect.v[d][0] = -inf; inf_rect.v[d][1] = inf; } vector<Area> area; area.pb(make_pair(inf_rect, 0)); rep(fi, n) { Rect fish; rep(i, 2) rep(d, 3) cin >> fish.v[d][i]; vector<Area> narea; foreach (it, area) { Rect and_rect; if (intersect(fish, it->first, and_rect)) { narea.pb(make_pair(and_rect, it->second + 1)); int bound[3][4]; rep(d, 3) { bound[d][0] = -inf; bound[d][1] = and_rect.v[d][0]; bound[d][2] = and_rect.v[d][1]; bound[d][3] = inf; } rep(i, 3) rep(j, 3) rep(k, 3) { if (i == 1 && j == 1 && k == 1) { // excluding: temp_rect == and_rect continue; } Rect temp_rect; temp_rect.v[0][0] = bound[0][i]; temp_rect.v[0][1] = bound[0][i + 1]; temp_rect.v[1][0] = bound[1][j]; temp_rect.v[1][1] = bound[1][j + 1]; temp_rect.v[2][0] = bound[2][k]; temp_rect.v[2][1] = bound[2][k + 1]; Rect and_r; if (intersect(temp_rect, it->first, and_r)) narea.pb(make_pair(and_r, it->second)); } } else narea.pb(*it); } area = narea; } ll res = 0; foreach (it, area) if (it->second >= k) res += it->first.calc_area(); cout << res << endl; }
Facebook Hacker Cup 2013 Qualification Round
まだジャッジが終わっていないので、間違っていたらごめんなさい
Balanced Smileys
解法
適当に探索。メモ化してやると、状態数はO(|s|^2)なので余裕で間に合う
ソースコード
map<string, bool> dp; bool ok_char(char c) { return (isalpha(c) && islower(c)) || c == ' ' || c == ':'; } bool dfs(const string& s) { if (dp.count(s)) return dp[s]; bool res = false; if (ok_char(s[0])) res |= dfs(s.substr(1)); if (s.size() >= 2) { if (ok_char(s[s.size() - 1])) res |= dfs(s.substr(0, s.size() - 1)); if (s.substr(0, 2) == ":)" || s.substr(0, 2) == ":(") res |= dfs(s.substr(2)); if (s[0] == '(') { for (int i = 1; i < s.size(); ++i) { if (s[i] == ')') res |= dfs(s.substr(1, i - 1)) && dfs(s.substr(i + 1)); } } } return dp[s] = res; } int main() { dp[""] = true; string s; getline(cin, s); int T = to_T<int>(s); for (int case_num = 1; case_num <= T; ++case_num) { getline(cin, s); string res = dfs(s) ? "YES" : "NO"; printf("Case #%d: %s\n", case_num, res.c_str()); } }
Find the Min
解法
k以降の連続したk+1個の要素を見ると、0~kの数がそれぞれちょうど1つしかない。周期k+1になっているので、k~(2*k+1)までの値を求めれば、大きいnの場合でもO(1)で求めることができる。
k~(2*k+1)を求める方法は、愚直にやって全体でO(k^2)でも6分あるので間に合いそうだけど、BITを使うと速く求められる。
ある要素を求めるときに、k個前に含まれる場合に1、含まれない場合に0を割り当てる。そうすると、始めに0となる数を探せばいい。これは二分探索でできる。
ソースコード
template <class T> class BIT { public: vector<T> a; int n; BIT(int n) : n(n), a(n + 1) {} BIT() { } void add(int i, T x) { ++i; assert(i > 0); while (i <= n) { a[i] += x; i += i & -i; } } T sum(int i) const { ++i; T res = 0; while (i > 0) { res += a[i]; i -= i & -i; } return res; } T range_sum(int low, int high) const { return sum(high) - sum(low - 1); } T at(int i) const { return sum(i) - sum(i - 1); } void assign(int i, T x) { add(i, x - at(i)); } }; // precondition: all elements are 0 or 1 int find_first_zero(const BIT<int>& bit) { int low = -1, high = bit.n; while (high - low > 1) { int mid = (low + high) / 2; if (bit.sum(mid) == mid + 1) low = mid; else high = mid; } return high; } vector<int> gen(int n, int k, ll a, ll b, ll c, ll r) { vector<int> v(n); v[0] = a; for (int i = 1; i < k; ++i) v[i] = (b * v[i - 1] + c) % r; const int mv = 3 * k + 100; BIT<int> bit(mv + 100); map<int, int> cnt; rep(i, k) { ++cnt[v[i]]; if (v[i] < mv) bit.assign(v[i], 1); } for (int i = k; i < n; ++i) { v[i] = find_first_zero(bit); int out = v[i - k]; int in = v[i]; --cnt[out]; ++cnt[in]; if (out < mv && cnt[out] == 0) bit.assign(out, 0); bit.assign(in, 1); } return v; } int solve(int n, int k, ll a, ll b, ll c, ll r) { const int num = 3 * k + 1000; vector<int> v = gen(num, k, a, b, c, r); --n; const int period = k + 1; int rebased_n = n - k; int peri = rebased_n % period; return v[k + peri]; } int main() { int T; cin >> T; for (int case_num = 1; case_num <= T; ++case_num) { int n, k; ll a, b, c, r; cin >> n >> k; cin >> a >> b >> c >> r; int res = solve(n, k, a, b, c, r); printf("Case #%d: %d\n", case_num, res); } }
SRM 566 Div2 Hard FencingPenguinsEasy
解法
あるフェンスが利用可能である条件は、フェンスを張ったときに左右のどちらか一方のみにペンギンが存在するような場合である。（ここでは、反時計回りにフェンスを張っていくので、左側に全てのペンギンがいる場合に利用可能とする。）
利用可能であるフェンスのみを用いて囲った時の面積を最小化にする方法を考える。
フェンスを張り始める杭を決める。この杭をstartとすると、フェンスを貼り終える杭もstartである。区別のためstart'とする。
i -> jにフェンスを張った時の面積の増加分は三角形(start, i, j)の面積で求められる。
区間[start, i]と区間[i, start']のフェンスの張り方は面積に相互作用しないので、区間[start, i]のフェンスの張り方を決定するときは、面積が最小となるものを選べば良い。よって、以下のdpを解けばいい。
dp[i] = [start, i]にフェンスを張った時の最小の面積として、
dp[j] = min{ dp[i] + triangle_area(start, i, j); フェンスi -> jは利用可能 }
求めたい解: dp[start'] = startから張り始めたときの最小の面積
以上のdpを全ての杭をstartとして計算し、その中で最小のものを全体の解とする。
幾何的な要素は外積を使うと楽
ソースコード
const double PI = acos(-1.0); double cross(double x1, double y1, double x2, double y2) { return x1 * y2 - y1 * x2; } double tri_area(double x1, double y1, double x2, double y2, double x3, double y3) { return abs(cross(x2 - x1, y2 - y1, x3 - x1, y3 - y1)) / 2; } double FencingPenguinsEasy::calculateMinArea(int numPosts, int radius, vector <int> x, vector <int> y) { const double sp = 2 * PI / numPosts; double px[300], py[300]; rep(i, numPosts) { px[i] = radius * cos(i * sp); py[i] = radius * sin(i * sp); } bool can_connect[300][300]; // i -> j rep(i, numPosts) rep(j, numPosts) { // ok if all penguins in left can_connect[i][j] = i != j; rep(k, x.size()) { // (post(i) -> post(j)) * (post(i) -> pen(k)) < 0 if (cross(px[j] - px[i], py[j] - py[i], x[k] - px[i], y[k] - py[i]) < 0) { can_connect[i][j] = false; break; } } } const double inf = 1e60; double res = inf; rep(start, numPosts) { double dp[333]; // [start, i] -> min area rep(i, numPosts) dp[i] = can_connect[start][i] ? 0 : inf; #define next(i) (((i) + 1) % numPosts) for (int i = next(start); i != start; i = next(i)) { for (int j = i + 1; j != i; j = next(j)) if (can_connect[i][j]) chmin(dp[j], dp[i] + tri_area(px[start], py[start], px[i], py[i], px[j], py[j])); } chmin(res, dp[start]); } return res == inf ? -1 : res; }
SRM 566 Div2 Medium PenguinPals
解法
区間dp
dp[i][j] = [i, j]の最大マッチング
ソースコード
int n; string c; int dp[55][55]; int dfs(int l, int r) { if (l >= r) return 0; else if (dp[l][r] != -1) return dp[l][r]; int res = 0; for (int i = l + 1; i <= r; ++i) if (c[l] == c[i]) chmax(res, 1 + dfs(l + 1, i - 1) + dfs(i + 1, r)); chmax(res, dfs(l + 1, r)); return dp[l][r] = res; } int PenguinPals::findMaximumMatching(string colors) { c = colors; n = colors.size(); CL(dp, -1); return dfs(0, n - 1); }
SRM 566 Div2 Easy PenguinTiles
解法
右下にある場合0
一番下の行 or 一番右の列にある場合1
それ以外の場合は1回の操作で↑の状態にできるので2
ソースコード
int PenguinTiles::minMoves(vector <string> tiles) { int h = tiles.size(), w = tiles[0].size(); if (tiles[h - 1][w - 1] == '.') return 0; rep(x, w) if (tiles[h - 1][x] == '.') return 1; rep(y, h) if (tiles[y][w - 1] == '.') return 1; return 2; }
ハル研究所プログラミングコンテスト2012
ハル研究所プログラミングコンテスト2012に参加しました。
応募締め切り時は総合9位でした。（最後に最適解が出ないコードを出してしまって少し不安、最大スコアのコードで順位付けてくれるんだろうか）
どのような方針を取ったか書いていきます。
まず、問題の制約を愚直に列挙すると
- 1 <= W, H <= 24
- 周期210
- アイテムの状態数6
- waitしたか
探索することを考えると、状態空間のサイズはO(W*H*210*6*2)で、普通に競プロでもでそうな制約です。競技やってる人ならすぐに方針が見えたと思います。
初めに取ったアプローチはdjikstraです。状態を(座標, 周期, アイテムの状態, waitしたか)として、コストはソースコードにあったスコア計算のマイナスになる部分を使いました。
ただし、このままではサーバで20secに間に合わないので、A*を使います。ヒューリスティックス関数には針を無視したゴールからの距離を用います。これでジャッジサーバで19secほどでました。
最適解が出るようになったので、stage数を10^6ほどにしてPCを温めます。
その結果わかったことは、最適解の場合はdamageを一度も受けない、waitはしないです。
これから、状態のwaitしたかを省けます。探索時にdamageを受ける場合はそれ以降探索する必要がありません。
最終的な方針は、状態を(座標, 周期, アイテムの状態)。コストは針を無視した最短路から外れた回数としました。具体的なコストは最短路に沿って移動した場合に0, アイテムを使った場合に1, 最短路から外れた場合に2となります。
コストが0~2だけなので、priority_queueを使わずに、queue q[コスト];のように複数のキューを用意しました。小さいコストのキューから見ていき、空になったらコスト+1のキューに進みます。なので、コード的にはbfsっぽい感じになっています。（考え方的にはA*と一緒です）
その他高速化のためにしたこと
キューにつっこむ状態は32bitに詰める
コストなどを格納しておくテーブルはint -> charにする
テーブル配列の次元をできるだけ下げる。例えば、table[y][x] -> table[1024]として、table[(y << 5) | x]でアクセスする
一番速かった実行時間は3.29secでした。自分のチューニング力だとこの辺が限界です。
この方針で一番ネックになっていたのは、周期の状態数が大きいことでしたが、自分で解決案が出せませんでした。
トップ層の人たちは状態を(座標, アイテムの状態, コスト)として、状態の値を各周期で到達可能かを持つ方針みたいです。
