ABC 137 F - Polynomial Construction

解法

想定解はかなり頭がいいが、今回は汎用的なラグランジュ補間で解く。
今 \(n\) 次の多項式 \(P(x)\) が \(P(x_i) = y_i ~(i = 0, \ldots, n)\) を満たすとする。このとき、$$P(x) = \sum_{i = 0}^n y_i \frac{\prod_{k \neq i} (x - x_k)}{\prod_{k \neq i} (x_i - x_k)}$$ が成り立つ。$$Q_i = \frac{y_i}{\prod_{k \neq i} (x_i - x_k)}$$ と置いておく。
各 \(i\) に対して \(Q_i\) を求めるのは \(O(n)\) で可能。
\(\prod_{k \neq i} (x - x_k)\) については、予め $$\prod_{k=0}^n (x - x_k) = x^{n+1} + c_nx^n + \ldots + c_1x + c_0$$ を求めておく。
\(c_i\) は、\(\prod_k^m (x - x_k)\) を計算して係数を求めた後、\((x - x_{m + 1})\) を掛けると係数が簡単な漸化式で表せるので、\(O(n^2)\) で求められる。
各 \(i\) について考えるときは \((x - x_i)\) で割る必要があるが、高次のほうから $$ c_n^\prime = c_{n+1} (= 1) \\ c_j^\prime = c_{j+1} + x_ic_{j + 1}^\prime$$ と計算できる。
漸化式の意味は、\(c_j^\prime\) が \(\{x_0, \ldots, x_{i - 1}, x_{i + 1}, \ldots, x_n\}\) から \(n - j\) 個選んでかけ合わせたものの総和、ということを考えると、理解できると思う。
こうして $$\prod_{k \neq i} (x - x_k) = c_n^\prime x^n + \ldots + c_1^\prime x + c_0^\prime$$ が求まったので、これに \(Q_i\) をかけたものを各 \(i\) について足し合わせればOK。

以上より、\(O(n^2)\) でこの問題が解けた。

ソースコード

#include <bits/stdc++.h>

using ll = long long;

ll modpow(ll x, ll n, const ll mod) {
    ll res = 1;
    while(n > 0) {
        if(n & 1) (res *= x) %= mod;
        (x *= x) %= mod;
        n >>= 1;
    }
    return res;
}
ll inverse(ll x, const ll mod) {
    return modpow(x, mod - 2, mod);
}

std::vector<ll> lagrange_interpolation(std::vector<ll> xs, std::vector<ll> ys, const int m) {
    const int n = xs.size();
    for(int i = 0; i < n; ++i) {
        xs[i] %= m;
        ys[i] %= m;
    }
    std::vector<ll> all_c(n + 1);
    all_c[0] = 1;
    for(int i = 0; i < n; ++i) {
        std::vector<ll> nxt(n + 1);
        for(int j = 0; j < n; ++j) {
            nxt[j + 1] = all_c[j];
        }
        for(int j = 0; j < n; ++j) {
            nxt[j] = (m + nxt[j] - xs[i] * all_c[j] % m) % m;
        }
        all_c = std::move(nxt);
    }

    std::vector<ll> c(n);
    for(int i = 0; i < n; ++i) {
        ll qi = 1;
        for(int j = 0; j < n; ++j) {
            if(i == j) continue;
            qi = qi * (m + xs[i] - xs[j]) % m;
        }
        qi = inverse(qi, m) * ys[i] % m;

        auto tmp_c = all_c;
        for(int j = n - 1; j >= 0; --j) {
            c[j] = (c[j] + qi * tmp_c[j + 1]) % m;
            tmp_c[j] = (tmp_c[j] + tmp_c[j + 1] * xs[i]) % m;
        }
    }
    return c;
}

using namespace std;

int main() {
    int p; cin >> p;
    vector<ll> xs(p), ys(p);
    for(int i = 0; i < p; ++i) {
        xs[i] = i;
        cin >> ys[i];
    }
    auto res = lagrange_interpolation(xs, ys, p);
    for(int i = 0; i < p; ++i) {
        cout << res[i] << " \n"[i + 1 == p];
    }
}

感想

係数求めるやつは初めて書いたのでちょっと混乱した。

ICPC 2019 国内予選 参加記

はじめに

7/12 に行われた国内予選に SleepingDragon で出ました。
チームメイトはいつもの kazuma, nakano です。

コンテスト内容

  • 問題を読む前に bashrc に alias の設定だけする
  • A を読むとはいなので書く
    • 入力形式を見間違えていてサンプルが合わなかった
    • すぐ気がついたので直して AC (0:05:11)
  • B を kazuma が書いてる間に C を聞く
    • やるだけなので終わったら書くことに
    • B が AC (0:14:05) したので適当に書いて AC (0:24:37)
    • C の FA らしい。謎すぎる。
  • D を kazuma に投げられてしまう
    • nakano から解法を聞く
    • 理解できなかったが自信ありそうだったので理解しないまま書く
    • 理解しないまま書いたので無事破滅する
    • 理解していないのでデバッグができないため nakano に丸投げする (Fを考える
    • nakano が直す箇所を教えてくれて直して AC (1:32:28)
  • E を kazuma が裏で頑張ってくれていたので AC (1:44:48)
  • F をみんなで考える
    • (使う頂点集合、使わない頂点集合) の二部グラフを考えると、色が変わるのはカット辺だけになるから見通しがいいかなと思った
    • 無限に悩む
    • コンテスト終了(は???)

結果

5完で全体6位、大学内3位で予選通過。

感想

さすがに今年は仕事しなさすぎでびっくりした。
D で虚無ったのはまあ反省で、それ以上に F が解けなかったのがかなりつらい。
僕が提示した指針を引きずりすぎた(しみんなも引きずってた)ので戦犯といえば戦犯、仕方ないといえば仕方ない。
院試が終わったら精進再開したいね。
京大から4チーム通過したのは多分初めてっぽいので、めでたい。アジアでは倒す。

AOJ 1387 - String Puzzle

解法

各分割区間と一致する場所が、「より左にある」という制約が重要で、これにちゃんと気がつけば解ける。
これによって

 文字列のある位置の文字と別の位置の文字が、与えられた条件によって等しい
 ⇔ それぞれの位置について、条件をたどって等しいことが言える位置集合のなかで、最も左の位置が等しい

が成り立つ。
したがって、最初に与えられた文字を最も左端に寄せた時にどこに対応するかを計算すれば、各クエリに対して答えるのは簡単(そのクエリの位置を左端に寄せていって、対応する位置の文字を読むだけ)。

左端への寄せ方だが、雑にやると TLE する。
ある区間に一致する区間が、その区間とかなり overlap していると、対応する位置に細かく移していくと O(n) になるからである。
これの解決は簡単で、overlap の場合は一気に飛ばすようにすれば良いだけである。
愚直に一回のステップで左に動く大きさが \(y[i] - h[i]\) なので、これが今見ている区間を追い越すギリギリを求めれば良い。
すなわち、今の位置を \(p\) として、\(p\) が乗っている区間を \(i\) とすれば $$p = p - \left(\left\lfloor\frac{p - y[i]}{y[i] - h[i]}\right\rfloor + 1\right) \times (y[i] - h[i])$$ と更新すると良い。

これで O((a + q)b) で解ける。

ソースコード

#include <bits/stdc++.h>
using namespace std;

int main() {
    int n, a, b, q; cin >> n >> a >> b >> q;
    vector<int> x(a);
    vector<char> c(a);
    for(int i = 0; i < a; ++i) {
        cin >> x[i] >> c[i];
    }
    vector<int> y(b), h(b);
    for(int i = 0; i < b; ++i) {
        cin >> y[i] >> h[i];
    }

    auto find_left_most_pos = [&] (int p) {
        while(true) {
            const int i = upper_bound(begin(y), end(y), p) - begin(y) - 1;
            if(i < 0 || h[i] == 0) break;
            p -= ((p - y[i]) / (y[i] - h[i]) + 1) * (y[i] - h[i]);
        }
        return p;
    };
    map<int, char> ans;
    for(int i = 0; i < a; ++i) {
        ans[find_left_most_pos(x[i])] = c[i];
    }

    while(q--) {
        int z; cin >> z;
        z = find_left_most_pos(z);
        cout << (ans.count(z) ? ans[z] : '?');
    }
    cout << endl;
}

感想

こんなに簡単なのに本番だと解かれないんだから本当に怖い。
参加してたときも Standings みてかなりビビって手を付けなかった気がする。

ゼッケンドルフの定理と全探索

はじめに

友人と会話してて得られた知見で、共有したいと思います。

ゼッケンドルフの定理

ゼッケンドルフの定理は

任意の正の整数が、連続するフィボナッチ数を含まないような形で、相異なる1つ以上のフィボナッチ数の和として一意に表現できる

というものです。
これを使うと、全探索がオーダーレベルで効率良く行えることがあります。

全探索

たとえば、一般的にこういう問題を考えます。

\(N\) 個のデータが順序ついて並んでいて、その部分集合であって、ある性質を満たすものを求める。

もしこの「ある性質」が、連続する2つのデータを取り出すと成り立たないようなものであれば、ゼッケンドルフの定理からオーダーが改善できます。
(成り立たないは直接的ですが、最小化問題などのときに、2連続は考えなくてもよい、みたいな状況もありえそうです。)

\(N\) 個の点を \(N\) ビットで表現するとすると、連続2点をえらばないので、
これはゼッケンドルフ表現と捉えることができます。
このような選び方が自然数と一対一対応するので、例えば \(N \leq 40\) だとすると、大体 \(F_{40} = 10^8\) 程度の部分集合しかないことがわかります。
なので、全探索が間に合う(かもしれない)ということですね。

最後に

結局考えてた問題ではこの性質は使わずに解きました(というか、こんな性質はなかったので使えなかった、悲しい)。

AOJ 1191 - Rotate and Rewrite

解法

こういう問題に対する典型的な発想として、一文字ずつ構成していく、というものがある。
つまり、A の先頭位置と B の先頭位置を決めうった後、一文字ずつ作っていって、A, B それぞれをちょうど使い切れたら、それが答えの候補になる。

このために dp をする。dp は2段階からなる。

先に2段階目の dp を考える。A, B の先頭位置を s1, s2 と決めうつと、

 dp[i][j] := A[s1..s1+i), B[s2...s2+j) を使い切って一致させるときの最大長

となる。求めたい答えは dp[n][m] である。
遷移式は、各 l1, l2 に対して

 もし A[s1+i...s1+i+l1) と b[s2+j...s2+j+l2) を使い切って同じ文字にできるなら
 dp[i + l1][j + l2] = max(dp[i + l1][j + l2], dp[i][j] + 1)

となる。条件部分を判定するために、別の dp が必要。
これは A, B のどちらについても同じなので、A だけ考えると

 can_make[l][r][k] := A[l..r) を使い切って値 k に一致させられるか?

である。これを求めるために補助テーブル

 r_dp[l][r][i][j] := A[l..r) まで使い切って、書き換えルール i の j 文字目まで一致させられるか?

を同時に作っていく。これによって、遷移は

 ある k があって、r_dp[l][k][i][j - 1] かつ A[k..r) を使い切ってルール i の j 文字目にできるならば
 r_dp[l][r][i][j] = true

とかける。この後、r_dp の更新を can_make に伝搬させるため

 書き換えルール i について、A[l..r) から x1, ..., xk をすべて作れたら
 can_make[l][r][y] = true

と更新する。
更にその後、can_make の更新を r_dp に伝搬させるため、

 書き換えルール i の一文字目が A[l..r) を使い切って作れるなら
 r_dp[l][r][i][1] = true

と更新する。
これを r - l の小さい方からやっていくと求まる。
あとは can_make を使って最初の dp を行えば解ける。

計算量は O(n^3m^3 * 30 + (n^3 + m^3)r*10) ぐらいでヤバ過ぎるが、適当に枝刈りを入れると 2sec ぐらいで通った。
critical な枝刈り部分はソースコードのコメントに書いた。

ソースコード

#include <bits/stdc++.h>
using namespace std;

template<typename T>
std::vector<T> table(int n, T v) { return std::vector<T>(n, v); }

template <class... Args>
auto table(int n, Args... args) {
    auto val = table(args...);
    return std::vector<decltype(val)>(n, std::move(val));
}

constexpr int max_a = 30;
constexpr int max_k = 10;

auto calc_can_make(vector<int> const& a, vector<vector<int>> const& xs, vector<int> const& ys) {
    const int n = a.size() / 2, r_sz = xs.size();
    auto can_make = table(a.size(), a.size(), max_a + 1, false);
    // r_dp[i][j][k][l] := can match a[i..j) with xs[k][0..l) ?
    auto r_dp = table(a.size(), a.size(), r_sz, max_k + 1, false);
    for(int i = 0; i + 1 < n * 2; ++i) {
        can_make[i][i + 1][a[i]] = true;
    }
    for(int i = 0; i < n; ++i) {
        for(int j = 0; j < r_sz; ++j) {
            r_dp[i][i][j][0] = true;
        }
    }
    for(int len = 1; len <= n; ++len) {
        for(int l = 0; l + len < 2 * n; ++l) {
            const int r = l + len;
            for(int i = 0; i < r_sz; ++i) {
                for(int j = 1; j <= (int)xs[i].size(); ++j) {
                    for(int k = l; k < r; ++k) {
                        if(r_dp[l][k][i][j - 1] && can_make[k][r][xs[i][j - 1]]) {
                            r_dp[l][r][i][j] = true;
                        }
                    }
                }
            }
            for(int i = 0; i < r_sz; ++i) {
                if(r_dp[l][r][i][xs[i].size()]) {
                    can_make[l][r][ys[i]] = true;
                }
            }
            for(int i = 0; i < r_sz; ++i) {
                if(can_make[l][r][xs[i][0]]) {
                    r_dp[l][r][i][1] = true;
                }
            }
        }
    }
    return can_make;
}

int main() {
    int n, m, r;
    while(cin >> n >> m >> r, n) {
        vector<int> a(n * 2), b(m * 2);
        vector<vector<int>> xs(r);
        vector<int> ys(r);
        for(int i = 0; i < n; ++i) {
            cin >> a[i];
            a[i + n] = a[i];
        }
        for(int i = 0; i < m; ++i) {
            cin >> b[i];
            b[i + m] = b[i];
        }
        for(int i = 0; i < r; ++i) {
            int k; cin >> k;
            xs[i].resize(k);
            for(auto& x : xs[i]) cin >> x;
            cin >> ys[i];
        }

        auto can_make_a = calc_can_make(a, xs, ys);
        auto can_make_b = calc_can_make(b, xs, ys);
        int ans = -1;
        for(int s1 = 0; s1 < n; ++s1) {
            for(int s2 = 0; s2 < m; ++s2) {
                vector<vector<int>> dp(n + 1, vector<int>(m + 1, -1));
                dp[0][0] = 0;
                for(int i = 0; i < n; ++i) {
                    for(int j = 0; j < m; ++j) {
                        if(dp[i][j] == -1) continue;
                        for(int k = 1; k <= 30; ++k) {
                            for(int l1 = 1; l1 <= n - i; ++l1) {
                                if(!can_make_a[s1 + i][s1 + i + l1][k]) continue; // critical
                                for(int l2 = 1; l2 <= m - j; ++l2) {
                                    if(can_make_a[s1 + i][s1 + i + l1][k] && can_make_b[s2 + j][s2 + j + l2][k]) {
                                        dp[i + l1][j + l2] = max(dp[i + l1][j + l2], dp[i][j] + 1);
                                    }
                                }
                            }
                        }
                    }
                }
                ans = max(ans, dp[n][m]);
            }
        }

        cout << ans << endl;
    }
}

感想

いわれてみれば確かになあという感じなんだけど全く見えなかった。
やりたいことのゴールから逆算していくと見えるのかなあ。一生解ける気がしない。
あとオーダーヤバスギ(国内予選だからまあ良いんだろうけど…)。

AOJ 1158 - ICPC: Intelligent Congruent Partition of Chocolate

解法

\(_{36}C_{18}\) は流石に間に合わない。
しかし、2つの領域を合同かつ連結になるように選んでいけば、探索空間はそんなに大きくないように思える(実際大きくない)。
なので、今探索中の領域1に連結なチョコレートを付け加えたときに、領域2にも対応する位置のチョコレートを使うと確定させるような探索が良さそう。
そのためにまず、領域2が領域1に対してどれだけ回転しているか&裏返しかを決めておく。
このもとで、一番左上のチョコレートを領域1として使うと確定させ、そのチョコレートが領域2においてどこに対応するかを全部試す。
すると、チョコレートを領域1に付け加えたとき、それが領域2のどこに対応するかがわかるようになる。
領域1にチョコレートを付け加えるときに、連結になるようにしておけば自動的に領域2も連結になるのでOK。
途中でバッティングしたり範囲外に出たりせず全部使い切れれば YES とすればよい。

ソースコード

#include <bits/stdc++.h>
using namespace std;

using ll = long long;

constexpr int dy[4] = {0, 1, 0, -1};
constexpr int dx[4] = {1, 0, -1, 0};

int main() {
    int w, h;
    while(cin >> w >> h, w) {
        vector<vector<int>> r(h, vector<int>(w)), idx(h, vector<int>(w, -1));
        vector<int> xs, ys;
        for(int i = 0; i < h; ++i) {
            for(int j = 0; j < w; ++j) {
                cin >> r[i][j];
                if(r[i][j] == 1) {
                    idx[i][j] = ys.size();
                    ys.push_back(i);
                    xs.push_back(j);
                }
            }
        }
        const int n = xs.size();
        if(n & 1) {
            cout << "NO" << endl;
            continue;
        }

        auto in_range = [&] (int y, int x) {
            return 0 <= y && y < h && 0 <= x && x < w;
        };
        auto calc_dy = [] (int mir, int rot, int s_dy, int s_dx) {
            if(rot == 0)      return s_dy;
            else if(rot == 1) return (mir ? -s_dx : s_dx);
            else if(rot == 2) return -s_dy;
            else              return (mir ? s_dx : -s_dx);
        };
        auto calc_dx = [] (int mir, int rot, int s_dy, int s_dx) {
            if(rot == 0)      return (mir ? -s_dx : s_dx);
            else if(rot == 1) return s_dy;
            else if(rot == 2) return (mir ? s_dx : -s_dx);
            else              return -s_dy;
        };

        auto solve = [&] () {
            for(int j = 1; j < n; ++j) {
                for(int rot = 0; rot < 4; ++rot) {
                    for(int mir = 0; mir < 2; ++mir) {
                        unordered_set<ll> vis;
                        stack<pair<ll, ll>> st;
                        vis.insert(1);
                        st.emplace(1, 1LL << j);
                        while(!st.empty()) {
                            const ll used1 = st.top().first, used2 = st.top().second;
                            if(__builtin_popcountll(used1) * 2 == n) {
                                return true;
                            }
                            st.pop();
                            for(int k = 0; k < n; ++k) {
                                if(!(used1 & (1LL << k))) continue;
                                for(int d = 0; d < 4; ++d) {
                                    const int ny1 = ys[k] + dy[d], nx1 = xs[k] + dx[d];
                                    const int ny2 = ys[j] + calc_dy(mir, rot, ny1 - ys[0], nx1 - xs[0]);
                                    const int nx2 = xs[j] + calc_dx(mir, rot, ny1 - ys[0], nx1 - xs[0]);
                                    if(!in_range(ny1, nx1) || idx[ny1][nx1] == -1) continue;
                                    if(!in_range(ny2, nx2) || idx[ny2][nx2] == -1) continue;
                                    const int id1 = idx[ny1][nx1], id2 = idx[ny2][nx2];
                                    if((used1 | used2) & (1LL << id1) || (used1 | used2) & (1LL << id2) || id1 == id2) continue;
                                    if(vis.count(used1 | (1LL << id1))) continue;
                                    vis.insert(used1 | (1LL << id1));
                                    st.emplace(used1 | (1LL << id1), used2 | (1LL << id2));
                                }
                            }
                        }
                    }
                }
            }
            return false;
        };

        cout << (solve() ? "YES" : "NO") << endl;
    }
}

感想

国内予選ならではって感じの問題だった。

AOJ 2682 - Polygon Guards

解法

答えは 9 以下なので全探索で間に合う(えぇ…)。
前処理で、ある点に配置したときにどこが見えるようになるかを bit で管理すると判定が O(1) になる。
視線の線分が多角形内部に全部入っているかどうかの判定は適当にやるとハマる。
まず、その線分と多角形の交点を全部求める。
線分の両端から線分の内側にちょっと進んだところが多角形の内部にあるかをまず判定。
交点については、線分の方向ベクトルの2方向にちょっと進んだところが多角形内部にあるかで判定。
この2つをすべてクリアした場合のみ、その点が見える。

AOJ でも 2sec ぐらいかかったので本番で出たらもうちょっと高速化しないとダメかも。

ソースコード

#include <bits/stdc++.h>
using namespace std;

// ライブラリ部分は省略

int dfs(polygon const& ps, int i, ll cur_vis, ll used, vector<ll> const& vis) {
    const int n = ps.size();
    if(cur_vis == (1LL << n) - 1) return __builtin_popcountll(used);
    if(i == n) return 9;
    int res = dfs(ps, i + 1, cur_vis, used, vis);
    if(__builtin_popcountll(used) < 8) {
        res = min(res, dfs(ps, i + 1, cur_vis | vis[i], used | (1LL << i), vis));
    }
    return res;
}

int main() {
    int n; cin >> n;
    polygon ps(n);
    for(int i = 0; i < n; ++i) {
        int x, y; cin >> x >> y;
        ps[i] = point(x, y);
    }

    vector<ll> vis(n);
    for(int i = 0; i < n; ++i) {
        vis[i] = 1LL << i;
        for(int j = 0; j < n; ++j) {
            if(i == j) continue;
            segment s(ps[i], ps[j]);
            bool check = true;
            for(int k = 0; k < n; ++k) {
                segment t(ps[k], ps[(k + 1) % n]);
                if(!isis_ss(s, t)) continue;
                const auto delta = (ps[j] - ps[i]) / abs(ps[i] - ps[j]) * 0.001;
                for(auto p : is_ss(s, t)) {
                    if(abs(s.a - p) > eps && abs(s.b - p) > eps) {
                        check &= is_in_polygon(ps, p + delta) != 2 && is_in_polygon(ps, p - delta) != 2;
                    }
                }
                check &= is_in_polygon(ps, ps[i] + delta) != 2;
                check &= is_in_polygon(ps, ps[j] - delta) != 2;
            }
            if(check) {
                vis[i] |= 1LL << j;
            }
        }
    }

    cout << dfs(ps, 0, 0, 0, vis) << endl;
}

感想

この問題みたいな多角形だと、n / 4 以下でできることが示せるらしい。
一般の多角形なら n / 3 とのこと。
たしかに多角形を三角形 or 四角形で分割することを考えるとそんな気もする。