AIOJ

Introduction

Inversion transformations are useful for problems involving tangencies among multiple circles or lines. Using the properties of inversion transformations to solve a problem in the inverted plane can greatly simplify the calculations.

Definition

Given an inversion center $O$ and an inversion radius $R$ . If points $P$ and $P'$ in the plane satisfy:

Point $P'$ lies on the ray $\overrightarrow{OP}$ .
$|OP| \cdot |OP'| = R^2$

then $P$ and $P'$ are called inverse points of each other.

Explanation

The inversion of a point $P$ in the plane is shown below:

Inv1

Properties

The inverse of a point outside circle $O$ lies inside circle $O$ , and vice versa; the inverse of a point on circle $O$ is the point itself.
The inverse of a circle $A$ that does not pass through point $O$ is also a circle that does not pass through point $O$ .

Inv2
- Let circle $A$ have radius $r_1$ , and let its inverse circle $B$ have radius $r_2$ . Then:
  $r_2 = \frac{1}{2}\left(\frac{1}{|OA| - r_1} - \frac{1}{|OA| + r_1}\right) R^2$
Proof.

Inv3

By the definition of inversion:

$\begin{aligned} |OC|\cdot|OC'| &= (|OA|+r_1)\cdot(|OB|-r_2) = R^2 \\ |OD|\cdot|OD'| &= (|OA|-r_1)\cdot(|OB|+r_2) = R^2 \end{aligned}$

Eliminate $|OB|$ and solve the equation.
- Let the coordinates of $O$ be $(x_0, y_0)$ , those of $A$ be $x_1, y_1$ , and those of $B$ be $x_2, y_2$ . Then:
  $\begin{aligned} x_2 &= x_0 + \frac{|OB|}{|OA|} (x_1 - x_0) \\ y_2 &= y_0 + \frac{|OB|}{|OA|} (y_1 - y_0) \end{aligned}$
  $|OB|$ can be calculated while finding $r_2$ above.
The inverse of a circle $A$ passing through point $O$ is a straight line not passing through $O$ . Indeed, a point on circle $A$ arbitrarily close to $O$ has an inverse point arbitrarily far from $O$ .

Inv4
If two figures are tangent at a point other than $O$ , their inverses are also tangent.

Examples

ICPC 2013 Hangzhou Regional Contest: Problem of Apollonius

Problem statement

Find all circles that pass through a point outside two given circles and are tangent to both given circles.

Solution

The analytic-geometry solution seems rather difficult at first.

Take the point that the desired circle must pass through as the inversion center, with an arbitrary inversion radius. The inverse of the desired circle is a straight line by Property $3$ , and it is tangent to the inverses of the two given circles by Properties $2$ and $4$ .

Thus, after inversion, the problem becomes finding all common tangents of two circles.

After finding the common tangents, invert them back to the original plane.

Sample code

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <vector>
using namespace std;

constexpr double EPS = 1e-8;   // Precision coefficient
const double PI = acos(-1.0);  // π
constexpr int N = 4;

// Definition of a point.
struct Point {
  double x, y;

  Point(double x = 0, double y = 0) : x(x), y(y) {}

  bool operator<(Point A) const { return x == A.x ? y < A.y : x < A.x; }
};

// Definition of a vector.
using Vector = Point;

// Vector addition.
Vector operator+(Vector A, Vector B) { return Vector(A.x + B.x, A.y + B.y); }

// Vector subtraction.
Vector operator-(Vector A, Vector B) { return Vector(A.x - B.x, A.y - B.y); }

// Scalar multiplication of a vector.
Vector operator*(Vector A, double p) { return Vector(A.x * p, A.y * p); }

// Scalar division of a vector.
Vector operator/(Vector A, double p) { return Vector(A.x / p, A.y / p); }

// Compare with zero.
int dcmp(double x) {
  if (fabs(x) < EPS) return 0;
  return x < 0 ? -1 : 1;
}

// Dot product of two vectors.
double Dot(Vector A, Vector B) { return A.x * B.x + A.y * B.y; }

// Length of a vector.
double Length(Vector A) { return sqrt(Dot(A, A)); }

// Cross product of two vectors.
double Cross(Vector A, Vector B) { return A.x * B.y - A.y * B.x; }

// Projection of a point onto a line.
Point GetLineProjection(Point P, Point A, Point B) {
  Vector v = B - A;
  return A + v * (Dot(v, P - A) / Dot(v, v));
}

// A circle.
struct Circle {
  Point c;
  double r;

  Circle() : c(Point(0, 0)), r(0) {}

  Circle(Point c, double r = 0) : c(c), r(r) {}

  // Return the point corresponding to a polar angle.
  Point point(double a) { return Point(c.x + cos(a) * r, c.y + sin(a) * r); }
};

// Return the number of common tangents of two circles; -1 means infinitely many.
// a[i] and b[i] are the tangency points of the i-th tangent on circles A and B.
int getTangents(Circle A, Circle B, Point* a, Point* b) {
  int cnt = 0;
  if (A.r < B.r) {
    swap(A, B);
    swap(a, b);
  }
  double d2 =
      (A.c.x - B.c.x) * (A.c.x - B.c.x) + (A.c.y - B.c.y) * (A.c.y - B.c.y);
  double rdiff = A.r - B.r;
  double rsum = A.r + B.r;
  if (dcmp(d2 - rdiff * rdiff) < 0) return 0;  // One circle contains the other.

  double base = atan2(B.c.y - A.c.y, B.c.x - A.c.x);
  if (dcmp(d2) == 0 && dcmp(A.r - B.r) == 0) return -1;  // Infinitely many common tangents.
  if (dcmp(d2 - rdiff * rdiff) == 0) {  // Internal tangency: one tangent.
    a[cnt] = A.point(base);
    b[cnt] = B.point(base);
    ++cnt;
    return 1;
  }
  // External common tangents.
  double ang = acos(rdiff / sqrt(d2));
  a[cnt] = A.point(base + ang);
  b[cnt] = B.point(base + ang);
  ++cnt;
  a[cnt] = A.point(base - ang);
  b[cnt] = B.point(base - ang);
  ++cnt;
  if (dcmp(d2 - rsum * rsum) == 0) {  // One internal common tangent.
    a[cnt] = A.point(base);
    b[cnt] = B.point(PI + base);
    ++cnt;
  } else if (dcmp(d2 - rsum * rsum) > 0) {  // Two internal common tangents.
    double ang = acos(rsum / sqrt(d2));
    a[cnt] = A.point(base + ang);
    b[cnt] = B.point(PI + base + ang);
    ++cnt;
    a[cnt] = A.point(base - ang);
    b[cnt] = B.point(PI + base - ang);
    ++cnt;
  }
  return cnt;
}

// Invert circle A, which lies outside point O, into circle B; R is the inversion radius.
Circle Inversion_C2C(Point O, double R, Circle A) {
  double OA = Length(A.c - O);
  double RB = 0.5 * ((1 / (OA - A.r)) - (1 / (OA + A.r))) * R * R;
  double OB = OA * RB / A.r;
  double Bx = O.x + (A.c.x - O.x) * OB / OA;
  double By = O.y + (A.c.y - O.y) * OB / OA;
  return Circle(Point(Bx, By), RB);
}

// Invert a line into circle B passing through O; R is the inversion radius.
Circle Inversion_L2C(Point O, double R, Point A, Vector v) {
  Point P = GetLineProjection(O, A, A + v);
  double d = Length(O - P);
  double RB = R * R / (2 * d);
  Vector VB = (P - O) / d * RB;
  return Circle(O + VB, RB);
}

// Return true if points A and B lie on the same side of the line.
bool theSameSideOfLine(Point A, Point B, Point S, Vector v) {
  return dcmp(Cross(A - S, v)) * dcmp(Cross(B - S, v)) > 0;
}

int main() {
  int T;
  scanf("%d", &T);
  while (T--) {
    Circle A, B;
    Point P;
    scanf("%lf%lf%lf", &A.c.x, &A.c.y, &A.r);
    scanf("%lf%lf%lf", &B.c.x, &B.c.y, &B.r);
    scanf("%lf%lf", &P.x, &P.y);
    Circle NA = Inversion_C2C(P, 10, A);
    Circle NB = Inversion_C2C(P, 10, B);
    Point LA[N], LB[N];
    Circle ansC[N];
    int q = getTangents(NA, NB, LA, LB), ans = 0;
    for (int i = 0; i < q; ++i)
      if (theSameSideOfLine(NA.c, NB.c, LA[i], LB[i] - LA[i])) {
        if (!theSameSideOfLine(P, NA.c, LA[i], LB[i] - LA[i])) continue;
        ansC[ans++] = Inversion_L2C(P, 10, LA[i], LB[i] - LA[i]);
      }
    printf("%d\n", ans);
    for (int i = 0; i < ans; ++i) {
      printf("%.8f %.8f %.8f\n", ansC[i].c.x, ansC[i].c.y, ansC[i].r);
    }
  }

  return 0;
}

Exercises

ICPC 2017 Nanning Regional Online Contest: Finding the Radius for an Inserted Circle

CCPC 2017 Online Contest: The Designer

Inversion Transformation