これは何

ここの問題を解きながら幾何ライブラリを整備していて、あまりまとまった解説記事がないなあと思ったので、全21問について、良いリンクやポイントなどを自分用に簡単にメモします。
自分用なんでだいぶテキトーです。解説集というよりは解説が載ったリンクまとめみたいな感じです。あと数2Bの知識マシマシです。
実装は主にうしさんの幾何テンプレートを参考にしています。基本自分で解いて、どうしてもわからなかったらコードを読んで何やってるか理解しようとする感じで作っているので、ほぼ写経している部分と自分で書いた部分の両方あります。なのでたまに効率悪いこともしているかも...その時は指摘されたいです。
後、AOJの幾何コース全21問のネタバレを含んでいるので注意して下さい。
自分の幾何ライブラリ→ https://siro53.github.io/compro_library/geometry/geometry.hpp

前提

前提で使ってる関数や構造体は、以後当たり前のように使っていきます

点(位置ベクトル)を複素数型で扱う
実軸(real)をx軸、虚軸(imag)をy軸として見る

using Point = complex<double>;

内積、外積、点の回転など、よく使うものは関数化(正直単位ベクトルとか関数にする必要ない気がしてきた)
比較するときは、計算誤差を意識してEPSを使って判定(equal関数)

inline bool equal(const double &a, const double &b) {
    return fabs(a - b) < EPS;
}

// 単位ベクトル(unit vector)を求める
Point unitVector(const Point &a) { return a / abs(a); }

// 法線ベクトル(normal vector)を求める
// 90度回転した単位ベクトルをかける
// -90度がよければPoint(0, -1)をかける
Point normalVector(const Point &a) { return a * Point(0, 1); }

// 内積(dot product) : a・b = |a||b|cosΘ
double dot(const Point &a, const Point &b) {
    return (a.real() * b.real() + a.imag() * b.imag());
}

// 外積(cross product) : a×b = |a||b|sinΘ
double cross(const Point &a, const Point &b) {
    return (a.real() * b.imag() - a.imag() * b.real());
}

// 点pを反時計回りにtheta度回転
Point rotate(const Point &p, const double &theta) {
    return Point(cos(theta) * p.real() - sin(theta) * p.imag(),
                 sin(theta) * p.real() + cos(theta) * p.imag());
}

// ラジアン->度
double radianToDegree(const double &radian) { return radian * 180.0 / PI; }

// 度->ラジアン
double degreeToRadian(const double &degree) { return degree * PI / 180.0; }

円、直線、線分は構造体にする(完全にうしさんのパクり)(これかなり設計が天才みある、すごい)
線分は直線を継承してるけど、正直いる？？？って感じはする...(うまみがわかってない)

// Line : 直線を表す構造体
// b - a で直線・線分を表せる
struct Line {
    Point a, b;
    Line() = default;
    Line(Point a, Point b) : a(a), b(b) {}
    // Ax+By=C
    Line(double A, double B, double C) {
        if(equal(A, 0)) {
            a = Point(0, C / B), b = Point(1, C / B);
        } else if(equal(B, 0)) {
            b = Point(C / A, 0), b = Point(C / A, 1);
        } else {
            a = Point(0, C / B), b = Point(C / A, 0);
        }
    }
};

// Segment : 線分を表す構造体
// Lineと同じ
struct Segment : Line {
    Segment() = default;

    Segment(Point a, Point b) : Line(a, b) {}
};

// Circle : 円を表す構造体
// pが中心の位置ベクトル、rは半径
struct Circle {
    Point p;
    double r;

    Circle() = default;

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

問題

トピック1：点とベクトル

1. Projection(射影)

下の式変形より、内積を使えば求められる。
std::normは絶対値の2乗を求める関数。

// 射影(projection)
// 直線(線分)lに点pから引いた垂線の足を求める
Point projection(const Line &l, const Point &p) {
    double t = dot(p - l.a, l.a - l.b) / norm(l.a - l.b);
    return l.a + (l.a - l.b) * t;
}

Point projection(const Segment &l, const Point &p) {
    double t = dot(p - l.a, l.a - l.b) / norm(l.a - l.b);
    return l.a + (l.a - l.b) * t;
}

2. Reflection(反射)

まず射影を求めて、その点を$X$とする。
すると求める点の位置ベクトルは、$\overrightarrow{p_1P} + 2 \overrightarrow{p_{1}X}$。

// 反射(reflection)
// 直線lを対称軸として点pと線対称の位置にある点を求める
Point reflection(const Line &l, const Point &p) {
    return p + (projection(l, p) - p) * 2.0;
}

3. Counter-Clockwise(点の回転方向)

むずい
自分で書くよりSen君のこの記事の方が丁寧でわかりやすい(丸投げ)

// 点の回転方向
// 点a, b, cの位置関係について(aが基準点)
int ccw(const Point &a, Point b, Point c) {
    b -= a, c -= a;
    // 点a, b, c が
    // 反時計回りの時、
    if(cross(b, c) > EPS) {
        return 1;
    }
    // 時計回りの時、
    if(cross(b, c) < -EPS) {
        return -1;
    }
    // c, a, bがこの順番で同一直線上にある時、
    if(dot(b, c) < 0) {
        return 2;
    }
    // a, b, cがこの順番で同一直線上にある場合、
    if(norm(b) < norm(c)) {
        return -2;
    }
    // cが線分ab上にある場合、
    return 0;
}

トピック2：線分と直線

1. 平行・垂直

$\vec{a}$と$\vec{b}$が平行 $\leftrightarrow$ $\vec{a}$と$\vec{b}$の外積が0 ($θ$が180度なので)
$\vec{a}$と$\vec{b}$が垂直 $\leftrightarrow$ $\vec{a}$と$\vec{b}$の内積が0 ($θ$が90度なので)

// 2直線の直交判定 : a⊥b <=> dot(a, b) = 0
bool isOrthogonal(const Line &a, const Line &b) {
    return equal(dot(a.b - a.a, b.b - b.a), 0);
}
// 2直線の平行判定 : a//b <=> cross(a, b) = 0
bool isParallel(const Line &a, const Line &b) {
    return equal(cross(a.b - a.a, b.b - b.a), 0);
}

2. 交差判定

線分の交差判定はccwを使う
要するに、線分abに対して、点cと点dが反対側にあればよい
線分が重なってたり端点で交わっているときもtrueになってしまうことに注意
この不等式で全て網羅できてるのすごい

// 線分sと線分tが交差しているかどうか
bool isIntersect(const Segment &s, const Segment &t) {
    return ccw(s.a, s.b, t.a) * ccw(s.a, s.b, t.b) <= 0 &&
           ccw(t.a, t.b, s.a) * ccw(t.a, t.b, s.b) <= 0;
}

3.交点

そもそも交差してない場合は考えていないので、交差していない場合があるときは先に交差判定しなきゃだめ
直線が一致しているときは1点を返す(if文の中身)
それ以外は、下の画像のような導出により求められる。

// 直線s, tの交点の計算
Point crossPoint(const Line &s, const Line &t) {
    double d1 = cross(s.b - s.a, t.b - t.a);
    double d2 = cross(s.b - s.a, s.b - t.a);
    if(equal(abs(d1), 0) && equal(abs(d2), 0)) {
        return t.a;
    }
    return t.a + (t.b - t.a) * (d2 / d1);
}

// 線分s, tの交点の計算
Point crossPoint(const Segment &s, const Segment &t) {
    return crossPoint(Line(s), Line(t));
}

4.距離

線分と線分の距離の前に、線分と点の距離について考える
線分と点の距離の定義：点pから「線分lのどこか」への最短距離
線分lの端点を点a, bとすると、点pが明らかに点a側にある時は点aと点pの距離が最短。点b側にあるときも同じ。
それ以外は点と直線の距離と求め方は同じ。
よって線分と線分の距離は、4通りの線分と点の距離のうち一番小さいものである

// 線分lと点pの距離を求める
// 定義：点pから「線分lのどこか」への最短距離
double distanceBetweenSegmentAndPoint(const Segment &l, const Point &p) {
    if(dot(l.b - l.a, p - l.a) < EPS) {
        return abs(p - l.a);
    }
    if(dot(l.a - l.b, p - l.b) < EPS) {
        return abs(p - l.b);
    }
    return abs(cross(l.b - l.a, p - l.a)) / abs(l.b - l.a);
}

// 線分sとtの距離
double distanceBetweenSegments(const Segment &s, const Segment &t) {
    if(isIntersect(s, t)) {
        return (double)(0);
    }
    double ans = distanceBetweenSegmentAndPoint(s, t.a);
    chmin(ans, distanceBetweenSegmentAndPoint(s, t.b));
    chmin(ans, distanceBetweenSegmentAndPoint(t, s.a));
    chmin(ans, distanceBetweenSegmentAndPoint(t, s.b));
    return ans;
}

トピック3：多角形

1.面積

多角形の面積は、隣り合う点の外積の総和の1/2
外積は、絶対値を外せば符号付き面積であることを考えるとわかりやすい
参考：n点からなる多角形の面積を求める | イメージングソリューション

// 多角形の面積を求める
double PolygonArea(const vector<Point> &p) {
    double res = 0;
    int n = p.size();
    for(int i = 0; i < n - 1; i++) {
        res += cross(p[i], p[i + 1]);
    }
    res += cross(p[n - 1], p[0]);
    return res * 0.5;
}

2.凸性判定

多角形の点を1周して、3点ずつみる。左、真ん中、右の点を順にpre, now, nxtとする。
この時、pre, now, nxtの順で同一直線上に並んでいる || preから見て、now->nxtと時計回りになっていたらfalse(凹になっているので)

// 凸多角形かどうか
bool isConvex(const vector<Point> &p) {
    int n = p.size();
    int now, pre, nxt;
    for(int i = 0; i < n; i++) {
        pre = (i - 1 + n) % n;
        nxt = (i + 1) % n;
        now = i;
        if(ccw(p[pre], p[now], p[nxt]) == -1) {
            return false;
        }
    }
    return true;
}

3.多角形-点の包含

【第2回】点の多角形に対する内外判定｜【技業LOG】技術者が紹介するNTTPCのテクノロジー｜【公式】NTTPCのCrossing number algorithmを実装した

// 多角形gに点pが含まれているか?
// 含まれる:2, 辺上にある:1, 含まれない:0
int isContained(const vector<Point> &g, const Point &p) {
    bool in = false;
    int n = (int)g.size();
    for(int i = 0; i < n; i++) {
        Point a = g[i] - p, b = g[(i + 1) % n] - p;
        if(imag(a) > imag(b)) {
            swap(a, b);
        }
        if(imag(a) <= EPS && EPS < imag(b) && cross(a, b) < -EPS) {
            in = !in;
        }
        if(cross(a, b) == 0 && dot(a, b) <= 0) {
            return 1;
        }
    }
    return (in ? 2 : 0);
}

トピック4：凸多角形

1. 凸包（convex hull）

色々ネットで資料をググり倒したが、結局図が豊富で蟻本が一番分かりやすかったので蟻本を読みましょう
AOJのやつは問題文をちゃんと読まないと罠にハマるので注意

AOJの凸包やっと通った
まず同一直線上の点を漏らして1敗、下かつ左から出力するという条件を見落として2敗、大量に時間を溶かし申し上げた
— siro53 (@shiro537) 2020年3月5日

// 凸包 O(NlogN)
vector<Point> ConvexHull(vector<Point> &p) {
    int n = (int)p.size(), k = 0;
    sort(ALL(p), [](const Point &a, const Point &b) {
        return (real(a) != real(b) ? real(a) < real(b) : imag(a) < imag(b));
    });
    vector<Point> ch(2 * n);
    // 一直線上の3点を含める -> (< -EPS)
    // 含め無い -> (< EPS)
    for(int i = 0; i < n; ch[k++] = p[i++]) { // lower
        while(k >= 2 && cross(ch[k - 1] - ch[k - 2], p[i] - ch[k - 1]) < EPS)
            --k;
    }
    for(int i = n - 2, t = k + 1; i >= 0; ch[k++] = p[i--]) { // upper
        while(k >= t && cross(ch[k - 1] - ch[k - 2], p[i] - ch[k - 1]) < EPS)
            --k;
    }
    ch.resize(k - 1);
    return ch;
}

2.凸多角形の直径

まだ解いていません

3.凸多角形の切断

まだ解いていません

トピック5：点集合

1. 最近点対

まだ解いていません

トピック6：線分集合

1. 線分交差（マンハッタン幾何）

まだ解いていません

トピック7：円

1. 円の交差判定

【高校数学Ⅱ】2つの円の位置関係5分類、直交する2円 | 受験の月に書いてある通り、5パターンあるので5パターン全部場合分けして実装します

// 2つの円の交差判定
// 返り値は共通接線の数
int isIntersect(const Circle &c1, const Circle &c2) {
    double d = abs(c1.p - c2.p);
    // 2つの円が離れている場合
    if(d > c1.r + c2.r + EPS) {
        return 4;
    }
    // 外接している場合
    if(equal(d, c1.r + c2.r)) {
        return 3;
    }
    // 内接している場合
    if(equal(d, abs(c1.r - c2.r))) {
        return 1;
    }
    // 内包している場合
    if(d < abs(c1.r - c2.r) - EPS) {
        return 0;
    }
    return 2;
}

2. 三角形の内接円

このページを読めばok 三角形の内接円と傍接円 - Wikipedia
他にも、角の二等分線を2本求めてその交点を内心とする求め方もあるが、めっちゃバグらせたのでググってみたら答えがでた。こっちの方が実装軽そうだし、採用...
角の二等分線はこんな感じで、隣り合う2辺について垂直二等分線を求めてそれらの交点と、頂点を結べば求められる

// 三角形の内心
Circle inCircle(const Point &a, const Point &b, const Point &c) {
     D A = abs(b - c), B = abs(a - c), C = abs(a - b);
     Point p(A * real(a) + B * real(b) + C * real(c),
                 A * imag(a) + B * imag(b) + C * imag(c));
     p /= (A + B + C);
     D r = distanceBetweenLineAndPoint(Line(a, b), p);
     return Circle(p, r);
}

3. 三角形の外接円

4. 円と直線の交点

交点を持たない、接する、2点で交わるの3通りで場合分け
接する時は、接点は射影そのもの
2点で交わる時は、射影を使えば後は法線ベクトルで実装は楽

// 円cと直線lの交点
vector<Point> crossPoint(const Circle &c, const Line &l) {
    vector<Point> res;
    double d = distanceBetweenLineAndPoint(l, c.p);
    // 交点を持たない
    if(d > c.r + EPS) {
        return res;
    }
    // 接する
    Point h = projection(l, c.p);
    if(equal(d, c.r)) {
        res.emplace_back(h);
        return res;
    }
    Point e = unitVector(l.b - l.a);
    double ph = sqrt(c.r * c.r - d * d);
    res.emplace_back(h - e * ph);
    res.emplace_back(h + e * ph);
    return res;
}

5. 円の交点

1と同じで同じで5通り場合分け
2点で交わる場合は、以下のようにして求めている

// 2つの円の交点
vector<Point> crossPoint(const Circle &c1, const Circle &c2) {
    vector<Point> res;
    int mode = isIntersect(c1, c2);
    // 2つの中心の距離
    double d = abs(c1.p - c2.p);
    // 2円が離れている場合
    if(mode == 4) {
        return res;
    }
    // 1つの円がもう1つの円に内包されている場合
    if(mode == 0) {
        return res;
    }
    // 2円が外接する場合
    if(mode == 3) {
        double t = c1.r / (c1.r + c2.r);
        res.emplace_back(c1.p + (c2.p - c1.p) * t);
        return res;
    }
    // 内接している場合
    if(mode == 1) {
        if(c2.r < c1.r - EPS) {
            res.emplace_back(c1.p + (c2.p - c1.p) * (c1.r / d));
        } else {
            res.emplace_back(c2.p + (c1.p - c2.p) * (c2.r / d));
        }
        return res;
    }
    // 2円が重なる場合
    double rc1 = (c1.r * c1.r + d * d - c2.r * c2.r) / (2 * d);
    double rs1 = sqrt(c1.r * c1.r - rc1 * rc1);
    if(c1.r - abs(rc1) < EPS) {
        rs1 = 0;
    }
    Point e12 = (c2.p - c1.p) / abs(c2.p - c1.p);
    res.emplace_back(c1.p + rc1 * e12 + rs1 * e12 * Point(0, 1));
    res.emplace_back(c1.p + rc1 * e12 + rs1 * e12 * Point(0, -1));
    return res;
}

6. 円の接線

点pが円の外にあること前提
下の図より、青い部分の長さが等しいので、実は点Pを中心とした円との交点が接点になり、実装をサボれる

// 点pを通る円cの接線
// 2本あるので、接点のみを返す
vector<Point> tangentToCircle(const Point &p, const Circle &c) {
    return crossPoint(c, Circle(p, sqrt(norm(c.p - p) - c.r * c.r)));
}

7.円の共通接線

イメージ図　hはcosθのこと

vector<Line> tangent(const Circle &a, const Circle &b) {
    vector<Line> ret;
    // 2円の中心間の距離
    double g = abs(a.p - b.p);
    // 円が内包されている場合
    if(equal(g, 0)) {
        return ret;
    }
    Point u = unitVector(b.p - a.p);
    Point v = rotate(u, PI / 2);
    for(int s : {-1, 1}) {
        double h = (a.r + b.r * s) / g;
        if(equal(h * h, 1)) {
            ret.emplace_back(a.p + (h > 0 ? u : -u) * a.r,
                             a.p + (h > 0 ? u : -u) * a.r + v);

        } else if(1 - h * h > 0) {
            Point U = u * h, V = v * sqrt(1 - h * h);
            ret.emplace_back(a.p + (U + V) * a.r, b.p - (U + V) * (b.r * s));
            ret.emplace_back(a.p + (U - V) * a.r, b.p - (U - V) * (b.r * s));
        }
    }
    return ret;
}

8. 円と多角形の共通部分

まだ解いていません

有益リンク集

感想

幾何はしんどい　解いてない5問も解き次第随時追加していきます

memo46

競プロの精進記録その他。

幾何ライブラリについてのメモ

これは何

前提

問題

トピック1：点とベクトル

1. Projection(射影)

2. Reflection(反射)

3. Counter-Clockwise(点の回転方向)

トピック2：線分と直線

1. 平行・垂直

2. 交差判定

3.交点

4.距離

トピック3：多角形

1.面積

2.凸性判定

3.多角形-点の包含

トピック4：凸多角形

1. 凸包（convex hull）

2.凸多角形の直径

3.凸多角形の切断

トピック5：点集合

1. 最近点対

トピック6：線分集合

1. 線分交差（マンハッタン幾何）

トピック7：円

1. 円の交差判定

2. 三角形の内接円

3. 三角形の外接円

4. 円と直線の交点

5. 円の交点

6. 円の接線

7.円の共通接線

8. 円と多角形の共通部分

有益リンク集

感想