推荐理由:

细节较多的计算几何,适合一做。

Describe

在平面上有三个没有公共部分的圆,求平面上一点使得到三个圆的切线的夹角相等。

Solution

容易得出:

稍加转换:

那么我们只要找到一个点 ,满足他到 个点的距离符合给定比例。

我们设出

平方一下,然后乘以分母。可以得到

我们把 看做常数,那么原方程可以化成一个

有点常识的人就可以知道这是圆的一般方程,然后再分 种情况讨论。

那么 的轨迹是一条直线。

只需要把三个圆两两向交,再分开判断就行了。

code

#include <cstdio>
#include <cmath>
#include <iostream>
using namespace std;
const bool LINE = true;
const bool CIRCLE = false;
const double EPS = 1e-8;
struct Point { double x, y;};
struct Node {
    bool type; double r;
    Point O, p[2];
} node[3];
int x[3], y[3], r[3];
double dis(const Point & p1, const Point & p2) {
    return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
}
void getget(int i, int j, int k) {
    if(r[i] == r[j]) { node[k].type = LINE;
        node[k].p[0].x = (x[i] + x[j]) / 2.0;
        node[k].p[0].y = (y[i] + y[j]) / 2.0;
        double dx = x[j] - x[i], dy = y[j] - y[i];
        double temp = dx; dx = -dy, dy = temp;
        node[k].p[1].x = node[k].p[0].x + dx;
        node[k].p[1].y = node[k].p[0].y + dy;
    } else { node[k].type = CIRCLE;
        if(r[i] > r[j]) swap(i, j);
        double dx = x[j] - x[i], dy = y[j] - y[i];
        double d = sqrt(dx * dx + dy * dy);
        node[k].p[0].x = x[i] + dx * r[i] / (r[i] + r[j]);
        node[k].p[0].y = y[i] + dy * r[i] / (r[i] + r[j]);
        node[k].p[1].x = x[i] - dx * r[i] / (r[j] - r[i]);
        node[k].p[1].y = y[i] - dy * r[i] / (r[j] - r[i]);
        node[k].r = dis(node[k].p[0], node[k].p[1]) / 2.0;
        node[k].O.x = (node[k].p[0].x + node[k].p[1].x) / 2.0;
        node[k].O.y = (node[k].p[0].y + node[k].p[1].y) / 2.0;
    }
} bool zero(double x) { return ((x > 0 ? x : -x) < EPS);}
bool work(const Point & u1, const Point & u2, const Point & v1, const Point & v2) {
    return zero((u1.x - u2.x) * (v1.y - v2.y) - (v1.x - v2.x) * (u1.y - u2.y));
}
Point solve(const Point & u1, const Point & u2, const Point & v1, const Point & v2) { Point ret = u1;
    double t = ((u1.x - v1.x) * (v1.y - v2.y) - (u1.y - v1.y) * (v1.x - v2.x)) / ((u1.x - u2.x) * (v1.y - v2.y) - (u1.y - u2.y) * (v1.x - v2.x));
    ret.x += (u2.x - u1.x) * t, ret.y += (u2.y - u1.y) * t;
    return ret;
}
void yuan(const Point & c, double r, const Point & l1, const Point & l2, Point & p1, Point & p2) {
    Point p = c;
    p.x += l1.y - l2.y;
    p.y += l2.x - l1.x;
    p = solve(p, c, l1, l2);
    double t = sqrt(r * r - dis(p, c) * dis(p, c)) / dis(l1, l2);
    p1.x = p.x + (l2.x - l1.x) * t;
    p1.y = p.y + (l2.y - l1.y) * t;
    p2.x = p.x - (l2.x - l1.x) * t;
    p2.y = p.y - (l2.y - l1.y) * t;
}
void find(const Point & c1, double r1, const Point & c2, double r2, Point & p1, Point & p2) {
    Point u, v;
    double t = (1 + (r1 * r1 - r2 * r2) / dis(c1, c2) / dis(c1, c2)) / 2;
    u.x = c1.x + (c2.x - c1.x) * t;
    u.y = c1.y + (c2.y - c1.y) * t;
    v.x = u.x + c1.y - c2.y;
    v.y = u.y - c1.x + c2.x;
    yuan(c1, r1, u, v, p1, p2);
}
double xmult(const Point & p1, const Point & p2, const Point & p0) {
    return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);
}
double Dis(const Point & p, const Point & l1, const Point & l2) {
    return fabs(xmult(p, l1, l2)) / dis(l1, l2);
}
int Get(int i, int j, Point p[2]) {
    if(node[i].type == LINE && node[j].type == LINE) {
        if(work(node[i].p[0], node[i].p[1], node[j].p[0], node[j].p[1])) {
            return 0;
        }
        p[0] = solve(node[i].p[0], node[i].p[1], node[j].p[0], node[j].p[1]);
        return 1;
    }
    if(node[i].type == CIRCLE && node[j].type == CIRCLE) {
        if(dis(node[i].O, node[j].O) - EPS >= node[i].r + node[j].r) return 0;
        if(dis(node[i].O, node[j].O) + EPS <= fabs(node[i].r - node[j].r)) return 0;
        find(node[i].O, node[i].r, node[j].O, node[j].r, p[0], p[1]);
        return 2;
    } if(node[i].type == CIRCLE) swap(i, j);
    if(Dis(node[j].O, node[i].p[0], node[i].p[1]) - EPS >= node[j].r) return 0;
    yuan(node[j].O, node[j].r, node[i].p[0], node[i].p[1], p[0], p[1]);
    return 2;
}
bool checkp(const Point &p) {
    if(node[2].type == LINE)
        return Dis(p, node[2].p[0], node[2].p[1]) <= EPS;
    else return fabs(dis(p, node[2].O) - node[2].r) <= EPS;
}
signed main() { Point pp;
    for (int i = 0; i < 3; ++ i) {
        scanf("%d %d %d", x + i, y + i, r + i);
    } pp.x = x[0], pp.y = y[0];
    for (int i = 0; i < 3; ++ i) {
        getget(i, (i + 1) % 3, i);
    } Point p[2]; int np = Get(0, 1, p);
    if(np == 1) printf("%.5lf %.5lf\n", p[0].x, p[0].y);
    else if(np == 2) {
        if(dis(p[0], pp) < dis(p[1], pp)) 
            printf("%.5lf %.5lf\n", p[0].x, p[0].y);
        else printf("%.5lf %.5lf\n", p[1].x, p[1].y);
    }
    return 0;
}