プログラミング雑記帳

online-judge-tooolsでAtCoderのテストケースを取得できない

verification-helper でAtCoderを使ってVerifyしようと思ったのですが、案内されるDropboxのアプリがおそらく追加のユーザーを受け付けていないため、うまくできませんでした。

同じ問題で Issue も立てられていますが、それに対する反応はありません。

そこで、Dropboxのアプリを自分で作成してアクセストークンを入手することにしました。

Dropboxアカウントの作成

省略します。

Dropboxアプリの作成

Dropbox Developers からアプリを作成します。

  1. Choose an API: Scoped access
  2. Choose the type of access you need: Full Dropbox
  3. Name your app: 適当に決めてください

App key と App secret の取得

アプリを作成するとアプリの設定ページが表示されます。 その中に App key と App secret という項目があると思うのでこれをメモしておきます。

追記

設定ページの下の方に Generated access token があります。 これを利用する場合は、Permission の変更だけすればよいです。 もしかしたら Permission の変更をしなくてもできるかもしれませんが、試していません。

Permission の変更

このままだと oj からファイルの読み込み等ができないので権限を変更します。 設定ページの Permissions タブから files.metadata.read と sharing.read にチェックを入れて保存します。 oj では Dropbox API v2 の /list_folder と /get_shared_link_file を使用しているらしいので、これらにチェックを入れています。

アプリへのログイン

https://www.dropbox.com/oauth2/authorize?client_id=${APP_KEY}&response_type=code にアクセスします(${APP_KEY} の部分は各自変更してください)。 適当に許可していくと、アクセスコードが生成されるのでこれもメモしておきます。

アクセストークンの取得

curl https://api.dropbox.com/oauth2/token --user ${APP_KEY}:${APP_SECRET} --data grant_type=authorization_code --data code=${ACCESS_CODE}

これを上でメモしたものに置き換えて実行します。 そうするとJSON形式で帰ってくるので、その中の access_token をメモしておきます。

たまに期限切れなどとでることがあるので、その時はアクセスコードをもう一度生成します。リロードで再生成されると思います。

テストケースを取得できるか確認

oj download https://atcoder.jp/contests/agc001/tasks/agc001_a --system --dropbox-token=${ACCESS_TOKEN}

これをメモしたものに置き換えて実行して、成功すれば問題ないです。 実行するとファイルが作られるので問題ない場所で実行してください。 実行する前に oj で AtCoder にログインしておいてください。

SteamのA=Bというゲームを全実績解除した

SteamでA=Bというプログラミングゲームを見つけたので、
やっていたところかなりハマっていつの間にか全実績解除していました。

store.steampowered.com

このゲームは、A=BでAをBに置換するという命令だけしかない言語(実際はもう少し命令がありますが…)を使って与えられた問題を解くだけのゲームです。
教授と学生の面白いやりとりもあります。

かなり難しいゲームなのですが、実直な解答が載ったサイトなどがあまり見つからないため、
GitHubに私の解答を上げました。
github.com

間違っている点や、良い解答などがあれば、プルリクかこのブログにコメントして頂けると嬉しいです。
解説のリクエストでも構いません。

VSCode + WSL + online-judge-tools で競プロ環境構築

競プロ

自分用のメモです。 とりあえずAtCoder用です。

WSLでg++が使えるようになっていることを前提としています。

ディレクトリ構成

とりあえず以下のような構成でファイルを作成します。

contest
├── .vscode
│   └── tasks.json
├── Pipfile
├── Pipfile.lock
└── atcoder
    ├── abc
    │   ├── abc001
    │   ├── abc002
    │   ・・・
    │   ├── abc236
    │   │   ├── a.cpp
    │   │   ├── b.cpp
    │   │   ├── c.cpp
    │   │   ├── d.cpp
    │   │   ├── e.cpp
    │   │   ├── f.cpp
    │   │   ├── g.cpp
    │   │   └── h.cpp
    │   ・・・
    │
    ├── agc
    ├── arc
    ├── dltest.sh
    └── othercontests

online-judge-tools のインストールと設定

pipenv を使ってインストールするので、そうでない人は適宜読み替えてください。

pipenv install online-judge-tools

インストールが終わったらAtCoderにログインします。 ユーザー名とパスワードが求められるので入れます。

pipenv run oj login https://atcoder.jp/

ログインできない場合(2022/03/20 追記)

うまくいかず止まることがあります。 そのときは、Seleniumをインストールするとうまくいくかもしれません。

pipenv install selenium
pipenv run oj login https://atcoder.jp/

サンプルケース取得用のシェルスクリプト作成

dltest.sh というファイル名で適当な場所(contest/atcoder/dltest.sh など)に作成します。

#!/bin/bash

contest_dir=$1
contest_id=${contest_dir##*/}
problem_id=$2
test_dir=${contest_dir}/test/${contest_id}_${problem_id}


if [ ! -e ${test_dir} ]; then
    pipenv run oj d -d ${test_dir} https://atcoder.jp/contests/${contest_id}/tasks/${contest_id}_${problem_id}
fi

作成したら、パーミッションがきちんと設定されているかを確認します。

tasks.json の設定

.vscode ディレクトリに tasks.json を作成して以下のように書きます。 pipenv で oj をインストールしていない場合は command や args の部分を自分の環境に合うように書き換えてください(いちいち pipenv run をせずに済む方法をご存知の方がいたらコメントください)。

{
    "version": "2.0.0",
    "tasks": [
        {
            "type": "cppbuild",
            "label": "build",
            "command": "/usr/bin/g++",
            "args": [
                "-std=c++17",
                "-O2",
                "-pedantic-errors",
                "-Wall",
                "-D_GLIBCXX_DEBUG",
                "-g",
                "${file}",
                "-o",
                "${fileDirname}/a.out"
            ],
            "problemMatcher": [
                "$gcc"
            ],
            "group": {
                "kind": "build",
                "isDefault": true
            },
            "detail": "compiler: /usr/bin/g++",
            "presentation": {
                "focus": true
            }
        },
        {
            "type": "shell",
            "label": "download test cases",
            "command": "${workspaceFolder}/atcoder/dltest.sh",
            "args": [
                "${fileDirname}",
                "${fileBasenameNoExtension}"
            ]
        },
        {
            "type": "shell",
            "label": "do oj test",
            "command": "pipenv",
            "args": [
                "run",
                "oj",
                "t",
                "-c",
                "${fileDirname}/a.out",
                "-d",
                "${fileDirname}/test/${fileDirnameBasename}_${fileBasenameNoExtension}"
            ],
            "dependsOn": [
                "build",
                "download test cases"
            ],
            "group": {
                "kind": "test",
                "isDefault": true
            },
            "presentation": {
                "focus": true
            }
        },
        {
            "type": "shell",
            "label": "submit",
            "command": "pipenv",
            "args": [
                "run",
                "oj",
                "s",
                "-y",
                "https://atcoder.jp/contests/${fileDirnameBasename}/tasks/${fileDirnameBasename}_${fileBasenameNoExtension}",
                "${file}"
            ]
        }
    ]
}

ショートカットキーの設定

デフォルトで Ctrl+Shift+B でbuildはできるようになっていますが、 testにはショートカットが設定されていません。 Ctrl+Shift+T を割り当てることにしますが、 View: Reopen Closed Editor に Ctrl+Shift+T が既に割り当てられているので、キーバインドを削除してから workbench.action.tasks.test に Ctrl+Shift+T を割り当てます。

tasks.json で設定した他のものは Ctrl+Shift+P から Tasks: Run Task から実行できます。

参考

qiita.com

qiita.com

ABC121 D - XOR World

競プロ

解いて解説を見たら自分の解法の百億倍スマートだったのですが、調べた感じ私と同じ解法が見つからなかったので何かに役立つことがあれば…

問題

atcoder.jp

解法

排他的論理和の性質から
$$f(A, B) = f(0, A - 1) \oplus f(0, B)$$
となります。
なので $f(0, N)$ を求めればよいことが分かります。そこで、簡単のため $g(N) := f(0, N)$ と書くことにします。

ここで、具体的に $N$ に値を入れて $0$ から考えると、$N = 2^n (n \geq 2)$ のときに $g(N) = N$ であることに気が付きます。
$g(N - 1) = 0$ であることを利用すると簡単に証明ができますが、証明は省略します。
※ ここで $N = 4n$ のときに $g(N) = N$ であることに気が付くと解説での解法にたどり着きます。

上記の性質を利用すると、
$N' = 2^n$ かつ $N' \leq N$ を満たすような最大の $N'$ での値は $g(N') = N'$ になるので、 $g(N' - 1) = 0$ となり
$$g(N) = N' \oplus (N' + 1) \oplus \dots \oplus N$$
が導かれます。
上式の項数の偶奇で $N$ の立っているビットの中で最上位のビットが $g(N)$ で立つかどうかが分かるので、最上位ビット以外の残りの計算は
\[
\begin{eqnarray}
&& (N' - N') \oplus (N' - N' + 1) \oplus \dots \oplus (N - N') \\
&=& 0 \oplus 1 \oplus \dots \oplus (N - N') \\
&=& g(N - N')
\end{eqnarray}
\]
となります。
したがって、再帰的に $g(N)$ を解くことができます。

注意点として、利用した性質が成り立つのは $N = 2^n (n \geq 2)$ のときのみなので $N \leq 3$ のときは場合分けする必要があります。

コード

下記のコードでは $N'$ をmsbという関数で求めています。
また、項数の偶奇を見るところでは $N'$ が偶数であることから $N$ の偶奇を見るだけでよいことを利用しています。
他にも、 $N - N'$ の部分が $N \oplus N'$ になっています。
atcoder.jp

抽象化遅延セグメント木の使い方

競プロ

自分用に、抽象化した遅延セグメント木の使い方をとりあえず書いておきます。 ここでは実際に幾つかの問題を解くことで使い方を示します。

抽象化遅延セグメント木の実装

まずはじめに、抽象化遅延セグメント木のコードを置いておきます。

template <typename T, typename E>
struct LazySegmentTree { // 0-indexed
    using F = function<T(T, T)>;
    using G = function<T(T, E)>;
    using H = function<E(E, E)>;
    
    int n, height;
    vector<T> tree;
    vector<E> lazy;
    const F f;
    const G g;
    const H h;
    const T ti;
    const E ei;

    LazySegmentTree() {}
    LazySegmentTree(F f, G g, H h, T ti, E ei) :f(f), g(g), h(h), ti(ti), ei(ei) {}

    void init(int n_) {
        n = 1, height = 0;
        while (n < n_) n *= 2, ++height;
        tree.assign(2 * n, ti);
        lazy.assign(2 * n, ei);
    }

    void build(const vector<T>& v) {
        int n_ = v.size();
        init(n_);
        for (int i = 0; i < n_; ++i) tree[n + i] = v[i];
        for (int i = n - 1; i > 0; --i) tree[i] = f(tree[2 * i], tree[2 * i + 1]);
    }

    inline T reflect(int k) { return (lazy[k] == ei?tree[k]:g(tree[k], lazy[k])); }

    inline void eval(int k) {
        if (lazy[k] == ei) return;
        lazy[2 * k] = h(lazy[2 * k], lazy[k]);
        lazy[2 * k + 1] = h(lazy[2 * k + 1], lazy[k]);
        tree[k] = reflect(k);
        lazy[k] = ei;
    }

    inline void thrust(int k) { for (int i = height; i > 0; --i) eval(k >> i); }

    inline void recalc(int k) {
        while (k >>= 1)
            tree[k] = f(reflect(2 * k), reflect(2 * k + 1));
    }
    
    void update(int s, int t, const E& x) { // [l, r)
        s += n, t += n;
        thrust(s), thrust(t - 1);
        int l = s, r = t;
        while (l < r) {
            if (l & 1) lazy[l] = h(lazy[l], x), ++l;
            if (r & 1) --r, lazy[r] = h(lazy[r], x);
            l >>= 1, r >>= 1;
        }
        recalc(s), recalc(t - 1);
    }

    T find(int s, int t) { // [l, r)
        s += n, t += n;
        thrust(s), thrust(t - 1);
        int l = s, r = t;
        T ll = ti, rr = ti;
        while (l < r) {
            if (l & 1) ll = f(ll, reflect(l++));
            if (r & 1) rr = f(rr, reflect(--r));
            l >>= 1, r >>= 1;
        }
        return f(ll, rr);
    }

    T at(int i) { return find(i, i + 1); }
};

詳しくは説明しませんが、 $T$ は取得クエリに必要な情報を表す型、 $E$ は更新クエリに必要な情報を表す型、 $f$ は $T$ と $T$ をマージする関数、 $g$ は $T$ に $E$ を作用させる関数、 $h$ は $E$ と $E$ をマージする関数、 $ti$ は $T$ の $f$ に関する単位元、 $ei$ は $E$ の $h$ に関する単位元です。

Range Minimum Query and Range Update Query(RMQ and RUQ)

これはAOJのこの問題です。 一応、概要を書いておくと、

数列 $a_0, a_1, a_2, \cdots, a_{n-1}$ に対し、

  • $a_s, a_{s+1}, \cdots, a_t$ を $x$ に変更する。
  • $a_s, a_{s+1}, \cdots, a_t$ の最小値を出力する。

この問題において取得クエリは $find(s, t)$ なので、取得クエリに必要な情報は最小値のみです。したがって $T$ はintとなります。 また、更新クエリは $update(s, t, x)$ なので、更新クエリに必要な情報は $x$ のみです。したがって $E$ もintとなります。

次に、関数 $f$ を考えます。 $A_{i,j} = min(a_i, a_{i+1}, \cdots, a_j)$ とすると、 $$A_{s,t} = min(A_{s,m}, A_{m+1,t}) \ (s \leq m \leq t)$$ と表せるので、 $$f(t1:T, t2:T) = min(t1, t2)$$ で良いことが分かります。

次に、関数 $g$ を考えます。 $update(s, t, x)$ のクエリのあとに、 $A_{s,t}$ がどのようになるかを考えると、 明らかに $A_{s,t} = x$ なので、 $$g(t1:T, e1:E) = e1$$ で良いことが分かります。 \begin{eqnarray} g(t1:T, e1:E) = \begin{cases} t1 \ (e1 = ei) \\ e1 \ (otherwise) \end{cases} \end{eqnarray} とする方が良いかもしれません。

次に、関数 $h$ を考えます。 $update(s, t, x_1)$ のクエリのあとに、 $update(s, t, x_2)$ のクエリが来たとき、 持つ必要がある情報は明らかに $x_2$ だけで十分なので、 $$h(e1:E, e2:E) = e2$$ で良いことが分かります(これは実装において、第二引数をあとに適用するようにしているためです)。 \begin{eqnarray} h(e1:E, e2:E) = \begin{cases} e1 \ (e2 = ei) \\ e2 \ (otherwise) \end{cases} \end{eqnarray} とする方が良いかもしれません。

最後に単位元 $ti$ と $ei$ を考えます。 intの $f$ に関する単位元は $INF$ としておくと都合が良いので $ti = INF$ などとします。しかし、今回の場合は $2^{31}-1$ で初期化されているので $$ti = 2^{31}-1$$ とします。 $h$ は代入演算とみなすことができ、更新クエリで与えられないような値に設定しておくとよいので $$ei = -1$$ などとします。

コード
int n, q; cin >> n >> q;
auto f = [](int a, int b) { return min(a, b); };
auto g = [](int a, int b) { return b; };
LazySegmentTree<int, int> seg(f, g, g, INT_MAX, -1);
seg.init(n);
while (q--) {
    int com, s, t, x; cin >> com >> s >> t;
    if (com) cout << seg.find(s, t + 1) << endl;
    else cin >> x, seg.update(s, t + 1, x);
}

https://onlinejudge.u-aizu.ac.jp/solutions/problem/DSL_2_F/review/4027109/S_ZK/C++14

Range Affine Range Sum

これはAtCoderのこの問題です。 一応、概要を書いておくと、

数列 $a_0, a_1, a_2, \cdots, a_{n-1}$ に対し、

  • $a_l, a_{l+1}, \cdots, a_{r-1}$ のそれぞれを $a_i = ba_i+c$ にする。
  • $\displaystyle \left(\sum_{i=l}^{r-1} a_i\right) mod \ 998244353$ を出力する。

この問題では $T$ や $E$ に区間の長さを持たせることで遅延セグメント木を使って解くことができます。また、クエリの区間が半開区間なのに注意します。

以下では $mod \ 998244353$ と明記しませんが、実装時には適切に書く必要があります。

$T$ を区間の合計値と区間の長さを保持する型とします。 $E$ を $a_i = ba_i+c$ の $b$ と $c$ を保持する型とします。

ここで、関数 $f$ を考えます。 $A_{[i,j)}^{sum}$ を $a_i, a_{i+1}, \cdots, a_{j-1}$ の合計、 $A_{[i,j)}^{len}$ を区間の長さとし、 $A_{i,j} = (A_{[i,j)}^{sum}, \ A_{[i,j)}^{len})$ とすると、 $$f(A_{l,m}:T, A_{m,r}:T) = (A_{[l,m)}^{sum}+A_{[m,r)}^{sum}, \ A_{[l,m)}^{len}+A_{[m,r)}^{len})$$ となります。

次に、関数 $g$ を考えます。 $a_l, a_{l+1}, \cdots, a_{r-1}$ は更新クエリのあと、 $$ba_l+c, \ ba_{l+1}+c, \cdots, \ ba_{r-1}+c$$ となるため、 $e_1 = (e_1^b, e_1^c) = (b, c)$ とすると、 $$g(A_{l,r}:T, e_1:E) = (e_1^bA_{[l,r)}^{sum}+e_1^cA_{[l,r)}^{len}, \ A_{[l,r)}^{len})$$ となります。

次に、関数 $h$ を考えます。 $e_1 = (b_1, c_1), e_2 = (b_2, c_2)$ とすると、 $a_l, a_{l+1}, \cdots, a_{r-1}$ は、更新クエリ $e_1$ のあとに $$b_1a_l+c_1, \ b_1a_{l+1}+c_1, \cdots, \ b_1a_{r-1}+c_1$$ となり、この次に更新クエリ $e_2$ が来ると、 $$b_2b_1a_l+b_2c_1+c_2, \ b_2b_1a_{l+1}+b_2c_1+c_2, \cdots, \ b_2b_1a_{r-1}+b_2c_1+c_2$$ となります。したがって、 $$h(e_1:E, e_2:E) = (e_2^be_1^b, \ e_2^be_1^c+e_2^c)$$ となります。

最後に単位元 $ti$ と $ei$ を考えます。 $T$ の $f$ に関する単位元は、 $T$ が合計値と長さを持つため、 $ti = (0, 0)$ となります。 $E$ の $h$ に関する単位元は、 $h(e_1:E, e_2:E) = e_1$ を満たすような $e_2$ となるため、 $ei = (1, 0)$ となります。

コード
using mint = ModInt<int, 998244353>;
int N, Q; cin >> N >> Q;
vector<pair<mint, mint> > A(N);
for (auto& [sum, len] : A) {
    cin >> sum;
    len = 1;
}
auto f = [](pair<mint, mint> a, pair<mint, mint> b) {
    return pair(a.first + b.first, a.second + b.second);
};
auto g = [](pair<mint, mint> a, pair<mint, mint> b) {
    return pair(b.first * a.first + b.second * a.second, a.second);
};
auto h = [](pair<mint, mint> a, pair<mint, mint> b) {
    return pair(b.first * a.first, b.first * a.second + b.second);
};
LazySegmentTree<pair<mint, mint>, pair<mint, mint> > seg(f, g, h, pair(0, 0), pair(1, 0));
seg.build(A);
while (Q--) {
    int f, l, r; cin >> f >> l >> r;
    if (f) cout << seg.find(l, r).first << '\n';
    else {
        mint b, c; cin >> b >> c;
        seg.update(l, r, pair(b, c));
    }
}

https://atcoder.jp/contests/practice2/submissions/18946367

ModInt は自分で用意してください...

幾何ライブラリ 解説2

競プロ

前回

幾何ライブラリ 解説1 - プログラミング雑記帳

距離

点と直線の距離

高校数学なのですが、射影を使うよりも、外積を使って値を出す方が楽なので一応説明します。
まず、点を P 、直線上の異なる二点を L_1, L_2 とし、また求める解を d とします。
\vec{a} = \overrightarrow{L_1L_2}, \vec{b} = \overrightarrow{L_1P} とし、 \vec{a}, \vec{b} のなす角を \theta とすると、
\vec{a} \times \vec{b} = \|\vec{a}\|\|\vec{b}\|sin\theta より、
\displaystyle \frac{\vec{a} \times \vec{b}}{\|\vec{a}\|} = \|\vec{b}\|sin\theta = \|\vec{b}\|\frac{d}{\|\vec{b}\|} = d となります。

double getDistanceLP(Line l, Point p) {
	return abs(cross(l.p2 - l.p1, p - l.p1) / abs(l.p2 - l.p1));
}
点と線分の距離

まず、点を P 、線分の二端点を S_1, S_2 とします。
射影を求めた時にその点が線分上にあれば、点と直線の距離と同じようにできます。
しかし、その点が線分上にない場合、点と線分の二端点の近い方の距離が解となります。
この判定は内積を使い正負を見ることで分かります。
\overrightarrow{S_1S_2} \cdot \overrightarrow{S_1P} の値が負ならば P と S_1 との距離を求めればよいです。
同様に、 \overrightarrow{S_2S_1} \cdot \overrightarrow{S_2P} の値が負ならば P と S_2 との距離を求めればよいです。

double getDistanceSP(Segment s, Point p) {
	if (dot(s.p2 - s.p1, p - s.p1) < 0.0) return abs(p - s.p1);
	if (dot(s.p1 - s.p2, p - s.p2) < 0.0) return abs(p - s.p2);
	return getDistanceLP(s, p);
}
線分と線分の距離

https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/2/CGL_2_D
線分が交差していれば勿論0です。
そうでないときは線分と片側の線分の端点との距離をそれぞれ4通り求め、それの最小値を出せばよいです。

double getDistanceSS(Segment s1, Segment s2) {
	if (intersectSS(s1, s2)) return 0.0;
	return min({ getDistanceSP(s1, s2.p1), getDistanceSP(s1, s2.p2),
			   getDistanceSP(s2, s1.p1), getDistanceSP(s2, s1.p2) });
}

交点

直線と直線の交点

下の図を用いて説明しますが、直角や平行などを書き忘れたのでなんとなく分かってもらえると助かります。
一つ目の直線の任意の二点を L_1, L_2 とし、二つ目の直線の任意の二点を L_3, L_4 としています。
\displaystyle \vec{L_1} + \frac{d_2}{d_1}\overrightarrow{L_1L_2} が求めたい答えになることが図より分かります。
これは相似から簡単に証明することができます。
したがって、 \displaystyle \frac{d_2}{d_1} を求めます。
\displaystyle \frac{\overrightarrow{L_1L_3} \times \overrightarrow{L_3L_4}}{\overrightarrow{L_1L_2} \times \overrightarrow{L_3L_4}} = \frac{\|\overrightarrow{L_1L_3}\|\|\overrightarrow{L_3L_4}\|sin\varphi}{\|\overrightarrow{L_1L_2}\|\|\overrightarrow{L_3L_4}\|sin\theta} = \frac{\|\overrightarrow{L_1L_3}\|sin\varphi}{\|\overrightarrow{L_1L_2}\|sin\theta} = \frac{d_2}{d_1}
と求めることができます。
直線が一致する場合の処理は問題に依存するので適宜変えますが、別の問題で使うためその時にまた説明します。

上記の説明とコードがほんの少しだけ違いますが問題ないのでこのままにします。

Point getCrossPointLL(Line l1, Line l2) {
	double a = cross(l1.p2 - l1.p1, l2.p2 - l2.p1);
	double b = cross(l1.p2 - l1.p1, l1.p2 - l2.p1);
	if (abs(a) < EPS && abs(b) < EPS) return l2.p1;
	return l2.p1 + (l2.p2 - l2.p1) * (b / a);
}
線分と線分の交点

https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/2/CGL_2_C
交点を持つことを仮定している場合、直線と直線の交点で代用できます。
そうでない場合でも、線分が交差しているという仮定の下だと簡単に考えることができるので、それを書きます。
下の図を用います。
一つ目の線分の二端点を S_1, S_2 とし、二つ目の線分の二端点を S_3, S_4 としています。
PS_1 : PS_2 = t : 1 - t とすると、求める交点は \vec{S_1} + t\overrightarrow{S_1S_2} となります。
また相似より、 t : 1 - t = d_1 : d_2 なので、 \displaystyle t = \frac{d_1}{d_1 + d_2} と d_1, d_2 が分かれば求められます。
d_1, d_2 は既に説明した点と直線の距離を使うことで求められますが、除算の回数を減らすために関数を呼び出さずに似たことをしています。


Point getCrossPointSS(Segment s1, Segment s2) {
	Vector base = s2.p2 - s2.p1;
	double d1 = abs(cross(base, s1.p1 - s2.p1));
	double d2 = abs(cross(base, s1.p2 - s2.p1));
	return s1.p1 + (s1.p2 - s1.p1) * (d1 / (d1 + d2));
}
円と直線の交点

https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/7/CGL_7_D
高校数学なので説明は省略します。円の弦の長さを求めてちょっとやるだけです。

vector<Point> getCrossPointCL(Circle c, Line l) {
	vector<Point> ps;
	Vector pr = project(l, c.c);
	Vector e = (l.p2 - l.p1) / abs(l.p2 - l.p1);
	if (equals(getDistanceLP(l, c.c), c.r)) return vector<Point>{pr, pr};
	double base = sqrt(c.r * c.r - norm(pr - c.c));
	ps.push_back(pr + e * base); ps.push_back(pr - e * base);
	return ps;
}
円と線分の交点

線分と円の交点が2つある場合や接する場合、円と直線の交点を出力すればよいです。
そうでなければ、つまり交点が1つのとき、円と直線の2交点のうちのどちらかを出力します。これは、2交点のうちどちらの交点が線分上に
あるかを見ればよいです。

vector<Point> getCrossPointCS(Circle c, Segment s) {
	Line l(s);
	vector<Point> ps = getCrossPointCL(c, l);
	if (intersectCS(c, s) == 2) return ps;
	if (dot(l.p1 - ps[0], l.p2 - ps[0]) < 0) ps[1] = ps[0];
	else ps[0] = ps[1];
	return ps;
}
円と円の交点

https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/7/CGL_7_E
一つ目の円の中心点と半径を C_1, r_1、 二つ目の円の中心点と半径を C_2, r_2 とします。
また、求める2交点を P_1, P_2 とし、 a = \angle P_1C_1C_2 = \angle P_2C_1C_2 （これは合同な三角形であることから導けます)とします。
C_1 と C_2 の距離を d とすると、
余弦定理より、 \displaystyle cos\angle P_1C_1C_2 = \frac{d^2 + r_1^2 - r_2^2}{2dr_1} とできるので、これを逆三角関数に入れ、角度を出したものを a とします。
それと、 \overrightarrow{C_1C_2} と x 軸のなす角を t とします。
上記で求めた角度を用いることで、
\displaystyle \vec{P_1} = \vec{C_1} + \begin{bmatrix} r_1cos(t + a) \\ r_1sin(t + a) \\ \end{bmatrix},
\displaystyle \vec{P_2} = \vec{C_1} + \begin{bmatrix} r_1cos(t - a) \\ r_1sin(t - a) \\ \end{bmatrix}
と2交点を求めることができます。


vector<Point> getCrossPointCC(Circle c1, Circle c2) {
	double d = abs(c1.c - c2.c);
	double a = acos((c1.r * c1.r + d * d - c2.r * c2.r) / (2 * c1.r * d));
	double t = arg(c2.c - c1.c);
	vector<Point> ps;
	ps.push_back(c1.c + polar(c1.r, t + a)); ps.push_back(c1.c + polar(c1.r, t - a));
	return ps;
}

幾何ライブラリ 解説1

競プロ

幾何ライブラリの説明をしていきます。
ここで扱うものはAOJの計算幾何学に入っている問題です。
参照している問題に書いてある変数等の宣言はしないので、URLが書いてあるものについては一度見てください。
おそらく幾つかに分けることになると思います。
高校数学を前提知識としていますのでその部分の説明はあまりしません。
また、include文などを逐次記述していると膨れてしまうので、記述しません。
最低限理解できればいいと思って書いているので、式変形や変数の定義など書かない時が多々あります。

ここでは、一度記述したものについては再度記述することなく使用しますので、標準のヘッダから提供される関数等と混同しないように注意してください。

分かりにくい場合や説明が足りない場合、コメントを残してもらえれば、追記するかもしれません。

高校数学範囲外の話

・ asin(x), acos(x), atan(x)
下記のサイトの説明がすごく分かりやすいです。
Arctan x についてのみ atan2(y, x) という関数があります。
この関数は対辺と隣辺の比ではなく符号付の長さを引数にすることで定義域を広くしています。
定義域は [-\pi, \pi] となっています。
mathtrain.jp

・外積
二次元幾何以外を扱わない為、二次元幾何における外積の演算内容と簡単な意味のみ説明します。
\vec{a} = (a_x, a_y), \vec{b} = (b_x, b_y) とします。
この時外積 \vec{a} \times \vec{b} は、
\vec{a} \times \vec{b} = a_xb_y - a_yb_x = \|\vec{a}\|\|\vec{b}\|sin\theta
となります。プログラム上では二項目を使いますが、三項目も説明でよく使います。
外積の値は \vec{a}, \vec{b} によって作られる平行四辺形の面積となっています。これは三項目より明らかです。

テンプレート

共通して使うものを幾つか並べます。
コメントアウトで軽い説明をしています。

マクロ
const double EPS = 1e-10;
const double PI = acos(-1);
#define equals(a, b) (fabs((a) - (b)) < EPS) // 誤差を考慮した同値判定
構造体とそれに伴う基本演算

vector(動的配列) と Vector(数学のベクトル) 混乱しそうだな...

struct Point { // 点の構造体
	double x, y; // (x, y)

	Point() {}
	Point(double x, double y) :x(x), y(y) {}

	Point operator+(const Point& p) const { return Point(x + p.x, y + p.y); }
	Point operator-(const Point& p) const { return Point(x - p.x, y - p.y); }
	Point operator*(const double& k) const { return Point(x * k, y * k); }
	Point operator/(const double& k) const { return Point(x / k, y / k); }

	friend istream& operator>>(istream& is, Point& p) {
		is >> p.x >> p.y;
		return is;
	}

	bool operator==(const Point& p) const { return (fabs(x - p.x) < EPS && fabs(y - p.y) < EPS); }
	bool operator<(const Point& p) const { return (x != p.x ? x < p.x : y < p.y); }

	double norm() { return x * x + y * y; }
	double abs() { return sqrt(norm()); }
};

typedef Point Vector; // 点とベクトルは同じように表せる

double norm(Vector a) { return a.x * a.x + a.y * a.y; }
double abs(Vector a) { return sqrt(norm(a)); }
double dot(Vector a, Vector b) { return a.x * b.x + a.y * b.y; } // 内積
double cross(Vector a, Vector b) { return a.x * b.y - a.y * b.x; } // 外積

double arg(Vector p) { return atan2(p.y, p.x); } // 偏角
Point polar(double r, double a) { return Point(cos(a) * r, sin(a) * r); } // 極座標から直交座標への変換

bool isParallel(Vector a, Vector b) { return equals(cross(a, b), 0.0); } // 平行判定
bool isOrthogonal(Vector a, Vector b) { return equals(dot(a, b), 0.0); } // 直交判定

struct EndPoint { // 端点の構造体(マンハッタン幾何)
	Point p;
	int seg, st; // 点のid, 点の種別

	EndPoint() {}
	EndPoint(Point p, int seg, int st) :p(p), seg(seg), st(st) {}

	bool operator<(const EndPoint& ep) const {
		if (p.y == ep.p.y) return st < ep.st;
		return p.y < ep.p.y;
	}
};

struct Segment { // 線分の構造体
	Point p1, p2; // 2端点

	Segment() {}
	Segment(Point p1, Point p2) :p1(p1), p2(p2) {}

	friend istream& operator>>(istream& is, Segment& s) {
		is >> s.p1 >> s.p2;
		return is;
	}
};

typedef Segment Line; // 直線は線分に長さを無くしただけのもの

Point project(Segment s, Point p) { // 点の射影
	Vector base = s.p2 - s.p1;
	double r = dot(p - s.p1, base) / base.norm();
	return s.p1 + base * r;
}

Point reflect(Segment s, Point p) { // 点の反射
	return p + (project(s, p) - p) * 2.0;
}

struct Circle { // 円の構造体
	Point c; // 中心点
	double r; // 半径

	Circle() {}
	Circle(Point c, double r) :c(c), r(r) {}
};

typedef vector<Point> Polygon; // 多角形(順序付きの点の集合)

時計・反時計回り

https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/1/CGL_1_C

まず、
\vec{a} = \vec{p_1} -\vec{p_0}, \vec{b} = \vec{p_2} - \vec{p_0},
\vec{a}, \vec{b} のなす角を \theta とします。
この時、外積 \vec{a} \times \vec{b} = \|\vec{a}\|\|\vec{b}\|sin\theta より時計・反時計回りの判定ができます。
外積の値が正であれば、sin\theta が正であるので反時計回りと判定できます。
同様に、外積の値が負であれば、時計回りと判定できます。
これ以降どちらでもない場合、つまり外積が 0 の時について考えます。
この場合、三点が一直線上にあることが分かるのでもう少し詳しく分類していきます。
まず、内積 \vec{a} \cdot \vec{b} = \|\vec{a}\|\|\vec{b}\|cos\theta が負の時、 \theta = \pi であるので ONLINE_BACK となります。
そうでないときは \theta = 0 です。
この時 ONLINE_FRONT となるには \|\vec{a}\| < \|\vec{b}\| である必要があり、そうでなければ ON_SEGMENT となります。

static const int COUNTER_CLOCKWISE = 1;
static const int CLOCKWISE = -1;
static const int ONLINE_BACK = 2;
static const int ONLINE_FRONT = -2;
static const int ON_SEGMENT = 0;

int ccw(Point p0, Point p1, Point p2) {
	Vector a = p1 - p0, b = p2 - p0;
	if (cross(a, b) > EPS) return COUNTER_CLOCKWISE;
	if (cross(a, b) < -EPS) return CLOCKWISE;
	if (dot(a, b) < -EPS) return ONLINE_BACK;
	if (a.norm() < b.norm()) return ONLINE_FRONT;
	return ON_SEGMENT;
}

交差判定

線分同士の交差判定

https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/2/CGL_2_B

線分 s_1 を直線とみなすと、2つの領域に分けることができます。このとき、点 p_2, p_3 がこの2つの領域の別側に属している、つまりその2点が同じ領域に属していないことが必要条件となります。ただし、p_2, p_3 のうちどちらかが直線上にあるものも条件に含まれます。
同様に、線分 s_2 と 点 p_0, p_1 でも同じことが言えます。
この二つの条件を満たしているとき線分が交差していると判定できます。

bool intersectSS(Point p1, Point p2, Point p3, Point p4) {
	return (ccw(p1, p2, p3) * ccw(p1, p2, p4) <= 0 &&
			ccw(p3, p4, p1) * ccw(p3, p4, p2) <= 0);
}
bool intersectSS(Segment s1, Segment s2) {
	return intersectSS(s1.p1, s1.p2, s2.p1, s2.p2);
}
円と線分の交差判定

線分が円に接する場合は交点が2つあると考えることにします。
まず交点を持たない条件について考えます。
線分を直線とみなしたとき、その直線と円の中心点との距離が円の半径より大きい場合確実に交点を持ちません。これは射影より求められます。この射影で求めた点を C します。
また、円の内部に線分が内包されるときも勿論交点を持ちませんこれは線分の両端点と円の中心点との距離が円の半径より小さい場合と言い換えることができます。
上記の条件だけでは十分ではありません、円の外に線分が位置し、最初の条件を満たさない場合が存在するからです。
しかし、これ以上考えるのは難しいので、上記の条件を満たさなければ、交点を1つもつ条件と2つもつ条件を考え、そのどちらにも属さなかった場合に交点を持たないと判断することにします。
したがって、まず交点を1つもつ条件を考えます。
線分のどちらかの端点が円の内部に、もう片方の端点が円の外部にあるとき、交点をただ1つもつと言えます。これは線分の端点と円の中心点との距離を求めることで分かります。
次に交点を2つもつ条件を考えます。
上記で求めた点 C が線分上にある時に交点を2つもつといえます。
よって上記のどれにも当てはまらなかったときには交点をもちません。

int intersectCS(Circle c, Segment s) {
	if (norm(project(s, c.c) - c.c) - c.r * c.r > EPS) return 0;
	double d1 = abs(c.c - s.p1), d2 = abs(c.c - s.p2);
	if (d1 < c.r + EPS && d2 < c.r + EPS) return 0;
	if ((d1 < c.r - EPS && d2 > c.r + EPS) || (d1 > c.r + EPS && d2 < c.r - EPS)) return 1;
	Point h = project(s, c.c);
	if (dot(s.p1 - h, s.p2 - h) < 0) return 2;
	return 0;
}
円同士の交差判定

https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/7/CGL_7_A

まんま高校数学なので省略します。
戻り値は共通接線の数です。

int intersectCC(Circle c1, Circle c2) {
	if (c1.r < c2.r) swap(c1, c2);
	double d = abs(c1.c - c2.c);
	double r = c1.r + c2.r;
	if (equals(d, r)) return 3; // 外接する
	if (d > r) return 4; // 離れている
	if (equals(d + c2.r, c1.r)) return 1; // 内接する
	if (d + c2.r < c1.r) return 0; // 内包する
	return 2; // 交わる
}

ICPC 2019 Asia Taiwan 参加記

競プロ 参加記

あのICPCさんが問題文をパクるなんてことはないですよね？？？？？？？？？？

０日目

起床ACして待ち合わせ場所に行って空港へGO
空港に弊大学の2チームで集まって機内へ

僕と tekihei さんだけ他４人とシートが離れていた
飛行機でゲームができたので一通りやる
ソリティアやってWindows7に思いを馳せていたところ謎のチャットが届きチャット機能に気付く
先輩がゆるキャン△をおすすめしてきたので横にいる tekihei さんが視聴
ゆるキャン△はいいぞ。

台湾に到着して電車でホテルまで移動
ok さんと同室でした
風呂場に脱衣所的な場所がなく人権を得られないことに気付く、他の二つの部屋を見に行くと人権があるのでクソ、死亡

１日目

ホテルのアレで吉野家で朝飯が食えるみたいなので吉野家で朝飯を食べる
紅茶か味噌汁どっちか選べとかいう謎の選択肢がでた
紅茶を頼んだはずなのにレモネードみたいな味がするものが出てきた、謎

昼飯に油そばを食べ practice へ
油そばなので勿論おいしい。
僕は名古屋県民なので歌志軒って油そば屋を推しておきます。
www.kajiken.biz

practice

Tシャツとかばんを貰う。普段使いしにくい黄色Tシャツらしいがしりません使います。
学生証をカバンに入れたままクロークに預けてしまう、雑魚！ｗ
会場へ行く、時計が二つあるな！？なんか十秒くらいずれてるな？パソコンの時計とは一分以上ずれてるな？は？？？？？？？？？？？？？？

いつの間にかpracticeが始まる
先輩たちがいろいろ試す。試し終えてA,B問題を通す、問題を任せられたので通そうと思ったら頭がおかしくなって弱実装ができなくなって終わる、明日大丈夫か僕みたいな気持ちになった

practiceが終わって夕飯が渡されたので帰る
ThinkMETが餃子を買ってきていて一つ貰う、うめーーーーー
渡された夕飯を食べる、異常に甘いものが多いのですが？？？？夕飯を食べたあと ok さんがおかしくなっていたのでやばいですね！

DDCCに参加して3完する、Cを通せないの実装力不足すぎないですか

例のごとく人権のない風呂に入って就寝。

２日目

朝飯は勿論吉野家です。日本人は味噌汁を食います。

本番

予定時刻近くになっても始まる様子がない（ずれてる２つの時計はすでに開始時刻をすぎている）
しばらくするとアナウンスがあって十五分延期
本番のサイトを見ると三十分延期した開始時刻が書いてあり困惑する。
ちょっとして再読み込みしたら十五分延期した時刻になっていたのでよし
どうこうしてる間に始まる
僕と mo3tthi さんで問題文を読む間 ok さんがテンプレートを写経
写し終えたところで mo3thhi さんから全探索のCがとんでくる
やる
WA
問題文を何も理解していないのに書いてだすのは頭弱すぎるし痛すぎるきちんと問題概要を聞きましょう
流石に次の提出で通す、申し訳なさでいっぱいのなか問題文を読み進める。Lの幾何はうまく考察すればできそうだなぁとか思いながらとりあえず読めるところまで読む
いろいろしてる間に先輩たちがD, K を通す
その後問題共有をしてLを考察する
しばらくして先輩たちが H, J を通す
Lがそれっぽい解法が浮かんだのでライブラリを写経する。その間に全員書ける感じになったのでダメそうなら交代って感じでやっていく
僕が書いたコードは写経ミスを直した後もバグる
終了30分前くらいになって考察抜けに気付いたのでそれを直す。バグる。
殆ど終わりの時にバグがとれた感じになったのでだす
TLE
終わった...

46位でした。僕以外の人類強すぎないか

banquet の時間になる。飯うまいのでうまい。ABCを無視して、ThinkMET に付いていき夜市に行く。サトウキビミルクとかいう飲み物を飲む勿論甘い。うまい。
持病の頻尿が突如襲ってきてやばくなる。あと少しで失禁するところだった、やばば...

ホテルに戻ると ok さんと mo3tthi さんがABCをしているので僕は非人権シャワーを浴びる。
就寝。

４日目

空港に向かい帰る。
その途中で昨日のコンテストのパクり疑惑のツイートが話に上がる


まんまこれじゃん...サンプルも似通ってるし...
終わりです。

帰宅。

ACPC2019 参加記

競プロ 参加記

ACPC2019(会津合宿)に参加したので

Day 1

自己紹介でTwitterでよく見る名前が挙がっていたのでおののいていた。とてもこわい。
スポンサーセッションとかなんかやって、うん
@mao_mao_2000 さんと @Haru_92_10 さんと組んだ

僕はBを解いて、その後Cの解法を相談してるときにクソ解法を僕が投げて実装させようとしたので無理ってなってしまった
咎人だ...
結局インターネッツを使ってそれっぽいアルゴリズム(SCC)調べて張り付けて解いた。


解説くん「ワ―シャルフロイドで解きます。」

その後、D考えたりE考えたりしてたら時間になった。
解けそうだったので悔しい

そういえばお菓子の所に胡瓜とカイワレがあったのおかしくないか？？？？？

Day 2

雛鶴あいちゃんとメイドさんが居ました

@DAyamaCTF さんと @wakuwinmail_C さんと組んだ

Day 1 は知り合いとでこの日に初めて他大学の人と組んだのもあってとても緊張してしまい、Aでバグらせてしまったので先にBを解いてもらった
その後Cの解法を聞いて任せてもらったが、数学がよわよわだったため、教えてもらったことを書いただけになってしまった。つらい
あとは二人を心の中で応援していたら8完していた。こわい

懇親会

すごいことになっていた。
会津勢ばかりつぶれていた気がする

Day 3

@pazzle1230 さんと @kzyKT_M さんと @albusSamurai さんと組んだ

Bを想定TLE解法でずっと考えていたうえにチームメイトと一緒になってバグらせたのでやばかったが、助けてもらって想定TLE解法で通してしまった
INF時間バグらせている一方で殆ど解いていた、やばいですね！！！

終わりに

ACPC初参加でとても楽しかったので来年までにもっと精進してまた参加したいですね