ブログ名

競技プログラミングやお歌のお話をする高菜です。

整数のまま行う偏角ソート(行列式のあれです。)の実装バリエーションの検討とご紹介です。

ねぼこさんのこちらの記事に触発され、私の方法もご紹介したい! という気持ちになりました。

nebocco.hatenablog.com

この記事の初めの方は Rust でのコーディングを主眼においてご説明しておりますが、最後の方に C++ 向けのお話も追記させていただきましたのでぜひです。

アイデア

ねぼこさんは偏角を  \lbrack 0, 90 ^ \circ \rbrack, \rbrack 90 ^ \circ, 180 ^ \circ \rbrack, \rbrack 180 ^ \circ, 360 ^ \circ \lbrack と分けていらっしゃるのですが、私は  \lbrack 0, 180 ^ \circ \lbrack, \lbrack 180 ^ \circ, 360 ^ \circ \lbrack などのように  180 ^ \circ ずつに分ける方がメリットが大きいと考えておりますので、こちらについてご紹介していこうと思います。この分け方の利点は何といっても原点との辞書順比較で判定できることです。座標の順番が逆であること(← これも後述するように、始線にこだわらないことで解消できます。またこのように末尾から比較する順序を、逆辞書順といいます。)に注意なのですが、 (y, x) \lt (0, 0)false なものが偏角が小さく、true のものが偏角が大きいです。下の図をご覧いただくとよいです。

(y, x) < 0

実装

結局逆辞順比較で原点と比較したあと、行列式でタイブレーキングすればよいわけでしたね。

さて、作りたい関数のシグニチャはこちらです。

fn argcmp(p0: (i64, i64), p1: (i64, i64)) -> Ordering;

アルゴリズムはこうです。

  • まずは、赤い領域と青い領域、どちらに入っているかをチェックです。bool には Ord が実装されていますから、((y0, x0) < (0, 0)).cmp(&((y1, x1) < (0, 0)))でよいですね。
  • 同じ領域に入った場合は determinant でチェックでしたね。これは (x1 * y0).cmp(&(x0 * y1)) でできます。
  • これは Rust の標準ライブラリ特有の事情ですが、タイブレーキングをするときには Ordering::then_with が便利ですから、使っていきましょう。

つまり、このように実装すればよいです。思ったより短いでしょう! 後ほどご紹介するように、始線にこだわらないことでさらに短くすることができます。

fn argcmp((x0, y0): (i64, i64), (x1, y1): (i64, i64)) -> Ordering {
    ((y0, x0) < (0, 0))
        .cmp(&((y1, x1) < (0, 0))
        .then_with(|| (x1 * y0).cmp(&(x0 * y1)))
}

then でもよいという説はありますが、諸説です。こうすると短くなる一方短絡評価のようにならず、 ((y0, x0)< (0, 0)).cmp(&((y1, x1) < (0, 0))Ordering::Equal な場合以外も (x1 * y0).cmp(&(x0 * y1)) が評価されてしまいます。このあたりはみなさまの好みです。

fn argcmp((x0, y0): (i64, i64), (x1, y1): (i64, i64)) -> Ordering {
    ((y0, x0) < (0, 0))
        .cmp(&((y1, x1) < (0, 0))
        .then((x1 * y0).cmp(&(x0 * y1)))
}

なおテストなのですが、こういう感じのテストがお手軽かつそこそこ強いですから、実装のお供におすすめです。

#[test]
fn test_argcmp() {
    use itertools::Itertools;
    let pts = [
        (1, 0),
        (1, 1),
        (0, 1),
        (-1, 1),
        (-1, 0),
        (-1, -1),
        (0, -1),
        (1, -1),
    ];
    for (i, j) in (0..8).cartesian_product(0..8) {
        assert_eq!(argcmp(pts[i], pts[j]), i.cmp(&j));
    }
}

始線がにこだわらない場合

始線の選択として、

  1. 何度から開始するか
  2. 開始角ピッタリのものは最初と最後どちらにいれるか

というものがあります。なお開集合・閉集合とのアナロジーで、最初に入れる場合は closed な始線、最後に入れる場合は open な始線と呼ぶことにしましょう。先程ご紹介したのは、「 0 ^ \circ-closed」です。

さて最も実装が簡単なのは、実は「 270 ^ \circ-open」です。この場合はなんとペアの比較を逆辞書順比較ではなく通常の辞書順で比較するため、覚えやすいです。

fn argcmp((x0, y0): (i64, i64), (x1, y1): (i64, i64)) -> Ordering {
    ((0, 0) < (x0, y0))
        .cmp(&((0, 0)< (x1, y1)))
        .then_with(|| (x1 * y0).cmp(&(x0 * y1)))
}

 270 ^ \circ-open はさすがにイヤですという方のために、 90 ^ \circ-open もご紹介しましょう。今度は逆辞書順でも通常の辞書順でもなく、辞書順の逆で比較します。

fn argcmp((x0, y0): (i64, i64), (x1, y1): (i64, i64)) -> Ordering {
    ((x0, y0)< (0, 0))
        .cmp(&((x1, y1)< (0, 0)))
        .then_with(|| (x1 * y0).cmp(&(x0 * y1)))
}

……ところで、エッ!? 逆辞書順と辞書順の逆は違うのですか!? 違うのです。逆辞書順の国語辞典末尾が「ア」のものが最初に記載される、いわゆる逆引き広辞苑方式ですが、辞書順の逆の国語辞典先頭が「ン」のもの(たとえば ンナガタカナ)が最初に記載される、ただ普通の広辞苑を逆から読んだだけのものになります。

2022-07-19 追記: 行列式によるソートと半空間への分解との二段階に分けることで、比較関数をスッキリさせる方法

上に紹介した方法を、<[_]>::sort_by() の中に直接書くと、次のようになるわけですが、これはラムダ式の中に実質的に条件分岐(Ordering::and_then)が入っており、速度面などであまり嬉しくないという説があります。

a.sort_by(|&(x0, y0), &(x1, y1)| {
    ((x0, y0)< (0, 0))
        .cmp(&((x1, y1)< (0, 0)))
        .then_with(|| (x1 * y0).cmp(&(x0 * y1)))
});

ただこれは単なるタイブレーキングですから、

  1. |&p| p < (0, 0)true のものを先に、false のものをあとになるように、入力配列を2つに分けます。
  2. 分けられた2つの部分それぞれを、|&(x0, y0), &(x1, y1)| (x1 * y0).cmp(&(x0 * y1)) によりソートする

という風に2段階に分けることで、スッキリ高速になります。ところで 1 を行うときに sort_unstable_by_key の代わりにいわゆる C++ の std::partition のようなものを使うこともできますね。また 2 を先に行ったあとで、キーが false だった部分と true だった部分それぞれについて 1. でソートするというのもありえます。Rust の標準ライブラリにはそれがないため、チチンプイプイ & ヒラケェーイ・ゴマ!!! で、あることにしましょう。Rust に <[T]>::partition: &mut [T] -> usize は、アル!!(暗示)(ないのでご自分で実装するか、<[_]>::sort() して Iterator::position() しましょう。ちなみにIterator::position() する作業は <_>::partition_point() を使うとやや楽ですが、残念ながら Rust 1.52.0 ですから、2022-07-19 時点の AtCoder 環境では使えません。)

let pos = partition(&mut a, |p| (0, 0) < p);
a[..pos].sort_by(|&(x0, y0), &(x1, y1)| x0 * y1 > x1 * y0);
a[pos..].sort_by(|&(x0, y0), &(x1, y1)| x0 * y1 > x1 * y0);

関数をお借りするだけお借りして、何も書かないのでは C++ さんが可愛そうですから。

auto pos = std::partition(std::begin(a), std::end(a), [](auto p) { return std::pair { 0, 0 } < p; }); 
std::sort(std::begin(a), pos, [](auto p0, auto p1) { return p0.first * p1.second > p0.second * p1.first; });
std::sort(pos, std::end(a), [](auto p0, auto p1) { return p0.first * p1.second > p0.second * p1.first; });

2022-07-19 追記: 行列式ではなく傾きを使う方法

こちらのツイートで教えていただきました。使用に注意を要しますし、パフォーマンス如何は私は未確認ですが、とてもおもしろいアルゴリズムです。 以下にご紹介するアルゴリズムは、私のアレンジを含みますから、元ネタの気になる方はリンク先へどうぞです。

まず偏角  \lbrack -\pi / 4, \pi / 4 \rbrack(境界付き、もちろん原点抜き)の部分のみをソートする方法を考えましょう。それは傾き  y / x をキーとしてソートすればよいです。さらに入力が制すであり上界  L を持つ、つまり  \lbrack 0 \le \lvert x \rvert \le y \le L \rbrack を満たすとしましょう。このとき十分大きな  M を用いて  \lfloor y M / x \rfloor をキーとしてソートすればよいです。

 M として例えば  L ^ 2 を取ればよいことは、相異なる点の傾きの差の最小値を見るとすぐにわかります。そうすると傾きが  1 以下であることからキーの最大値は  L ^ 2 以下にになります。また、計算途中の値  y M L ^ 3 以下です。

また細かいですが、負数を除算することがあるのですが、「負数の場合だけ floor ではなく ceil」というよくある仕様であったとしても実は全く同じ  M で動きます。なぜなら仕様の変わり目がちょうど floor/ceil したい有理数の整数点(具体的には  0)であることから、相異なる点の最小値で考えた条件でそのままうまくいくことがわかります。説明が面倒ですからみなさまで確認していただいてですね。(あの?)

なお負数除算についてはこちらです: Arithmetic operators - cppreference.com

The quotient is truncated towards zero (fractional part is discarded).

// L = 1_i128 << 32;
const M: i128 = 1_i128 << 64;
a.sort_by_key(|&(x, y)| y * M / x);

さて、他の象限は次のようにソートします。

  1. 偏角  \lbrack -\pi / 4, \pi / 4 \rbrack:  y / x
  2. 偏角  \lbrack \pi / 4, 3\pi / 4 \rbrack:  -x / y
  3. 偏角  \lbrack 3\pi / 4, 5\pi / 4 \rbrack:  -y / x
  4. 偏角  \lbrack 5\pi / 4, 7\pi / 4 \rbrack:  x / y

なお象限の境界はどちらに入っても構いません。ただ始線の closed/open には影響しますから、注意です。

まとめて書くと、このような感じです。なお引き続き partition 関数は所与のものとします。また分母側には負数が来ないように工夫しています。

const M: i128 = 1_i128 << 64;
let j = partition(&mut a, |&(_, y)| 0 <= y);
let i = partition(&mut a[..j], |&(x, _)| 0 < x);
let k = partition(&mut a[j..], |&(x, _)| x < 0);
a[..i].sort_by_key(|&(x, y)| y as i128 * M / x as i128);
a[i..j].sort_by_key(|&(x, y)| -x as i128 * M / y as i128);
a[j..k].sort_by_key(|&(x, y)| y as i128 * M / -x as i128);
a[k..].sort_by_key(|&(x, y)| -x as i128 * M / -y as i128);

繰り返しになりますが、入力制約に注意です。

なおこの方法ですと  -\pi/4 が始線になり、 0 を始線にするのは一手間がかかってしまいます。 斜めに分けずに象限で分けて、 y / x などの代わりに  (y - x) / (x + y) などを使うという手もあります。

2022-07-19 追記: 偏角ソート、同じ角度のものを区別するともうちょいややこしくなる説について

こちらです。

同じ偏角のものについては例えば、例えば原点に近いほうが先にしたいとします。さて、どうしましょうか。

答えは簡単でして、最初に原点からの距離(の二乗)でソートして、その後で安定ソートのみを用いた偏角ソートを行えばよいだけだったりします。

2022-07-18 追記: C++ など、less 関数を作る必要がある場合

こちらの Vi24E さんのツイートに触発されて追記してしまいました。

引き続き始線  90 ^ \circ (open) が短く書けますからそうしましょう。先ほどお話しした通り、最初にペアを辞書順(X 座標優先)の降順で このように、条件演算子を使って std::pair { 0, 0 } との大小比較をするとよいです。(軽めのテスト → https://wandbox.org/permlink/EPM4QPHfbaT6oKBA

bool argcmp_less(std::pair<int, int> p0, std::pair<int, int> p1) {
    return
        (std::pair {0, 0} < p0) == (std::pair {0, 0} < p1)
        ? p0.first * p1.second > p0.second * p1.first
        : p0 < p1;
}

先ほどお話した通り、始線  0 ^ \circ (closed) にこだわるならばペアの比較を逆辞書順(Y 座標優先)昇順に変更すると良いです。とやや長くなってしまいますが、こうです。(軽めのテスト → [C++] gcc 12.1.0 - Wandbox

bool argcmp_less(std::pair<int, int> p0, std::pair<int, int> p1) {
    auto [x0, y0] = p0;
    auto [x1, y1] = p1;
    return
        (std::pair {0, 0} < std::pair {y0, x0}) == (std::pair {0, 0} < std::pair {y1, x1})
        ? x0 * y1 > x1 * y0
        : std::pair {y0, x0} > std::pair {y1, x1};
}

なおですが、C++ ですとなおさら、タイブレーキングを比較関数の中で行わず、そもそもソートを2段階にするほうがおすすめではあります。

まとめ

要するに偏角ソートのバリエーションは2つの半空間への分け方で決まります。

  1. 始線  4 種類( 0 ^ \circ-,  90 ^ \circ-,  180 ^ \circ-,  270 ^ \circ-)
  2. 始線ピッタリを最初に持ってくるかどうか(-closed, -open)

です。今回ご紹介した限りですと、実装的には

  1. 辞書順か逆辞書順か
  2. 昇順か降順か

の 4 通りになってしまい、足りないわけですが、足りない分は「第一成分は昇順・第二成分は降順」のように混ぜるとちょうどすべて作れますのでぜひです。

また、最後は行列式の正負で判定するわけですが、この正負を変えると当然ソート全体の方向が変わります。

ちなみに浮動小数点数の  \mathop { \mathrm { atan2 } } などでソートする場合、浮動小数点数の誤差が怖いのは言わずもがなですが、実は  \mathop { \mathrm { atan2 } } ですと始線を  0 ^ \circ - closed にしないと壊れるという扱いづらさもあります。せっかく切った半空間の中に偏角の切れ目が来てしまったは大変です。整数の行列式を用いると、それを気にせず任意の場所で切って良いため実装が気楽というメリットもあるわけですね!