Circle Inversion - Yixiong Gao

Analysis

二维平面上，给定反演中心 \(O\) 和半径 \(r\) ，定义点 \(P\) 基于圆 \((O,r)\) 的反演变换：\(P\mapsto P' : |OP||OP'|=r^2\) 。

圆内点反演完在圆外，圆外点反演完在圆内，圆上点是反演变换的不动点。

进一步的我们可以定义图形基于圆 \(O\) 的反演（即变换后的点集），容易发现一些结论：

不过 \(O\) 的圆 \(\Leftrightarrow\) 不过 \(O\) 的圆（先求出在 \(O\) 与圆心连线上直径两个端点的反演点再恢复，注意圆心反演后不是圆心！！）
过 \(O\) 的直线 \(\Leftrightarrow\) 过 \(O\) 的直线（且是同一条直线，不会有变化）
过 \(O\) 的圆 \(\Leftrightarrow\) 不过 \(O\) 的直线（先求出 \(O\) 对应直径的另一个端点的反演点，再求垂直于 \(O\) 与反演点连线的直线）
- 圆内区域反演完之后是反演直线的不包含 \(O\) 的半平面
反演前两图形的相对关系，反演后不变（反演前相切，反演后也相切；相交等以此类推）
- 特殊的，两切于 \(O\) 的圆反演后是两平行直线（切点 \(O\) 反演后在无穷远处，即两平行直线交于无穷远处）

反演半径一般选取点到中心距离的几何平均。完整计算几何模版见这里。

struct C {
    P c; double r;
    C(letp &c = zero, double r = 0.0) : c(c), r(r) {}

    bool operator == (letc &a) const {return c == a.c && z(r - a.r);}
    
    // 点与圆的关系 : -1 圆上 | 0 圆外 | 1 圆内
    int is_in(letp &p) const {double d = p.dis(c); return z(d - r) ? -1 : d < r - eps;}
    
    P inverse(letp &p) const {
        assert(p != c);
        const P dlt = p - c;
        return c + dlt * (r * r / dlt.norm());
    }
 
    tuple<int, C, L> inverse(letl &l) const {
        letc null_c = C(zero, 0.0);
        letl null_l = L(zero, zero);
        if (l.ori(c) == 0) return mt(2, null_c, l);
        letp v = (l.ori(c) == 1 ? P{l.v.y, -l.v.x} : P{-l.v.y, l.v.x});
        const double d = r * r / l.dis(c);
        letp p = c + unit(v) * d;
        return mt(1, C{(c + p) / 2, d / 2}, null_l);
    }
         
    tuple<int, C, L> inverse(letc &a) const {
        letc null_c = C(zero, 0.0);
        letl null_l = L(zero, zero);
        letp v = a.c - c;
        if (a.is_in(c) == -1) {
            const double d = r * r / (a.r + a.r);
            return mt(2, null_c, L{c + unit(v) * d, {-v.y, v.x}});
        }
        if (c == a.c) return mt(1, C{c, r * r / a.r}, null_l);
        const double dis = abs(v);
        const double k1 = r * r / (dis - a.r), k2 = r * r / (dis + a.r);
        const double d = (k1 + k2) / 2, rad = (k1 - k2) / 2;
        return mt(1, C{c + v * (d / dis), rad}, null_l);
    }
};

CF 77 E - Martian Food

给定两个内切的圆，第三个圆位于两圆心连线且与两圆相切，第四个圆开始和前三个圆相切，以此类推。

求第 \(k\) 个圆的半径。多组 \(t\le 10^4,r,R\le 10^4,k\le 10^4\) 。

以两圆切点反演，前两个圆变成平行直线，后面的圆都与两平行线相切，第三个圆心反演后位于前两圆心连线。

然后朝一侧找到第 \(k\) 个圆，再反演回来即得到了答案。注意不能反演圆心，圆心反演完不是圆心！

inline void work() {
    static C O = C(zero, 20000);
    double r1 = rd(), r2 = rd();
    L l1 = get<2>(O.inverse(C(P(r1, 0), r1)));
    L l2 = get<2>(O.inverse(C(P(r2, 0), r2)));
    C res = C(P((l1.p.x + l2.p.x) / 2, l1.dis(l2) * rd()), l1.dis(l2) / 2);
    printf("%.10lf\n", get<1>(O.inverse(res)).r);
}

HDU 6158 - The Designer

求类似上一题图的前 \(n\ (n\le 10^7)\) 个圆的面积和，多组 \((t\le 10^3)\) 。

判一下如果剩下的圆面积和足够小就不用做了，注意不是面积足够小，而是面积和。

inline void work() {
    static C O = C(zero, 200);
    double r1 = rd(), r2 = rd();
    L l1 = get<2>(O.inverse(C(P(r1, 0), r1)));
    L l2 = get<2>(O.inverse(C(P(r2, 0), r2)));
    double dlt = l1.dis(l2);
    C res = C((l1.p + l2.p) / 2, dlt / 2);
    double ans = get<1>(O.inverse(res)).area();
    for (int i = 1, t = 0, n = rd(); i < n; ++i) {
        (i & 1) ? ++t, res.c.y = t * dlt : res.c.y = -t * dlt; 
        double add = get<1>(O.inverse(res)).area();
        if (add * (n - i) < 1e-7) break; ans += add;
    }
    printf("%.5lf\n", ans);
}

2017 ICPC Nanning Online - Finding the Radius for an Inserted Circle

三个半径为 \(R\) 的圆互相切，在中间的弧形边三角形中塞小圆，然后朝着一个角塞，问第 \(k\) 个的半径。

反演点就是会往那个角里塞小圆的两个圆的切点。

三个圆变成一对平行直线+一个夹中间的圆，然后朝反方向找第 \(k\) 个圆即可。

int ans[11];
int main() {
    int t; double r;
    scanf("%d%lf", &t, &r);
    C O = C(zero, r * 2);
    C c = get<1>(O.inverse(C(P(0, -sqrt(3) * r), r)));
    for (int i = 1; i <= 10; ++i) {c.c.y -= r * 4; ans[i] = floor(get<1>(O.inverse(c)).r + 1e-8);}
    for (int i = 1, n; i <= t; ++i) {scanf("%d", &n); printf("%d %d\n", n, ans[n]);}
    return 0;
}

2013 ICPC Hangzhou - Problem of Apollonius

求过指定点且与给定两个圆相外切的所有的圆。

关于指定点反演，所求答案就是两个圆的公切线。

inline void work() {
    C c1; c1.c.x = rd(); c1.c.y = rd(); c1.r = rd();
    C c2; c2.c.x = rd(); c2.c.y = rd(); c2.r = rd();
    C O; O.c.x = rd(); O.c.y = rd(); O.r = 100;
    vector<L> s = get<1>(O.inverse(c1)).tangent(get<1>(O.inverse(c2)));
    vector<C> ans;
    for (auto x : s) {
        C inv = get<1>(O.inverse(x));
        if (inv.relation(c1) == 1 && inv.relation(c2) == 1) ans.push_back(inv);
    }
    printf("%d\n", (int)ans.size());
    for (auto x : ans) printf("%.10lf %.10lf %.10lf\n", x.c.x, x.c.y, x.r);
}

CF 372 E - Drawing Circles is Fun

给定二维平面内 \(n\) 个点, 其中任意两点所在的直线上都不包含原点 \(O\) 。

求集合 \(\{\left(P_1, P_2\right)\left(P_3, P_4\right) \ldots\left(P_{2 k-1}, P_{2 k}\right)\}\) 的个数（\(P_i\) 取自给定的 \(n\) 个点），满足:

不存在 \(i \neq j, P_i=P_j\).

对于任意 \(i \neq j,\left(P_{2 i-1}, P_{2 i}\right)\left(P_{2 j-1}, P_{2 j}\right)\) 满足：

圆 \(O P_{2 i-1} P_{2 j-1}\) 和圆 \(O P_{2 i} P_{2 j}\) 只有一个交点, 圆 \(O P_{2 i-1} P_{2 j}\) 和圆 \(O P_{2 i} P_{2 j-1}\) 只有一个交点。

所有圆都过 \(O\) 且只有一个交点，所以都相切在 \(O\) 。

反演之后就是平行直线，即要求 \(P_{2 i-1} P_{2 j-1}\parallel P_{2 i} P_{2 j},\ P_{2 i-1} P_{2 j}\parallel P_{2 i} P_{2 j-1}\) ，即 \(P_{2 i-1} P_{2 j-1}P_{2 i} P_{2 j}\) 是平行四边形。

平行四边形的等价条件是对角中点重合且边不重合，换言之，计数选出若干点对，反演后中点重合在同一个点，且不能重合。

处理出来 \(n^2\) 个线段的中点。对每个中点，处理出来同角度的每个方向有多少个（即重合）。

相当于 \(n\) 个集合里挑若干个（多于一个），每个集合里只能拿至多一个的方案数，即 \(\prod (size_i+1) - \sum size_i - 1\) 。

坑1： 实现用 pair(中点，方向向量）保存信息，为了保证排序正确（中点相同的挨在一起，相同中点按照方向极角序），不能出现反向的方向向量（共线的必须同一方向），所以 \(x\) 为负的时候取反，\(x=0\) 且 \(y<0\) 的时候取反 。
坑2： 精度要求高，坐标范围 \([\frac{1}{50},50]\) ，反演半径最好取几何平均 \(1\) 。
坑3： 最后对集合还要 count 一次。

#define mod 1000000007
#define ppp pair<P, P>
int main() {
    int n = rd(); 
    vector<P> p(n);
    vector<ppp> s;
    for (int i = 0; i < n; ++i) {
        double a = rd(), b = rd(), c = rd(), d = rd();
        p[i].x = a / b; p[i].y = c / d; 
        p[i] = p[i] / p[i].norm();
    }
    for (int i = 0 ; i < n; ++i)
        for (int j = i + 1; j < n; ++j) {
            P u = p[i], v = p[j];
            P mid = (u + v) / 2, dir = u - v; 
            if (z(dir.x)) {if (dir.y < 0) dir = -dir;}
            else if (dir.x < -eps) dir = -dir; 
            s.push_back(mp(mid, dir));
        }
    sort(s.begin(), s.end(), [&](ppp &a, ppp &b) {
      return (a.fr != b.fr) ? a.fr < b.fr : argcmp()(a.sc, b.sc);
    });
    int ans = 0;
    vector<int> cnt;
    auto count = [&]() {
        int prod = 1, sum = 0;
        for (auto x : cnt) {
            prod = 1ll * prod * (x + 1) % mod;
            sum = (sum + x) % mod;
        }
        int res = (prod + mod - sum - 1) % mod;
        ans = (ans + res) % mod;
        cnt.clear();
    };
    for (size_t i = 0; i < s.size(); ++i) 
        if (i > 0 && s[i].fr == s[i - 1].fr) {
            if (para(s[i].sc, s[i - 1].sc)) ++cnt.back();
            else cnt.push_back(1);
        } else {count(); cnt.push_back(1);}
    count();
    printf("%d\n", ans);
    return 0;
}

2022 ICPC Hefei F - Rescue

给定 \(n\le 10^6\) 个圆，保证所有圆都经过原点。

问 \(n\) 个圆的交是否为空，若不为空，则求出交区域内的点离原点最远的距离。

关于原点反演，每个圆就变成了一个半平面，求出来半平面交，问题变成了点到凸包的最近距离。

Algorithms

Geometry

本博客所有文章除特别声明外，均采用 CC BY-SA 4.0 协议，转载请注明出处！