Inversion

typedef double LD;
const LD PI = 3.14159265358979323846;
const LD eps = 1E-10;
const LD R2 = 1.0;
int sgn(LD x) { return fabs(x) < eps ? 0 : (x > 0 ? 1 : -1); }
struct P {
    LD x, y;
    P(LD x = 0, LD y = 0): x(x), y(y) {}
    P operator * (LD k) { return P(x * k, y * k); }
    P operator / (LD k) { return P(x / k, y / k); }
    string prt() const {
        char s[100];
        sprintf(s, "(%.2f, %.2f)", x, y);
        return string(s);
    }
};
typedef P V;
P operator - (const P& a, const P& b) { return P(a.x - b.x, a.y - b.y); }
P operator + (const P& a, const P& b) { return P(a.x + b.x, a.y + b.y); }
struct C {
    P p;
    LD r;
    C(LD x = 0, LD y = 0, LD r = 0): p(x, y), r(r) {}
};
LD dist(V v) { return sqrt(v.x * v.x + v.y * v.y); }

C inv(C c, const P& o) {
    LD d = dist(c.p - o);
    assert(sgn(d) != 0);
    LD a = 1 / (d - c.r);
    LD b = 1 / (d + c.r);
    c.r = (a - b) / 2 * R2;
    c.p = o + (c.p - o) * ((a + b) * R2 / 2 / d);
    return c;
}