<|instruction|>
Solve the below problem. The solution should start with an abridged problem statement. Then key observations. Then full solution based on the observations. Then C++ and Python implementations with comments.

414. Orthogonal Circles
Time limit per test: 0.5 second(s)
Memory limit: 262144 kilobytes
input: standard
output: standard



Two circles are orthogonal, if they intersect and for any point of intersection their tangent lines at that point are perpendicular.

Consider a set of n circles on the plane. You are to find a circle orthogonal to all of them.

Input
The input file contains an integer n (1 ≤ n ≤ 105), followed by n triples of integers: xi, yi, ri, denoting the center coordinates and the radii of the circles (-106 ≤ xi, yi ≤ 106, 1 ≤ ri ≤ 106).

The circles may coincide.

Output
If there exists exactly one circle orthogonal to all the given ones, output its center coordinates and radius as real numbers separated with single spaces. The numbers will be considered correct if they are within 10-6 relative or absolute error of the exact answers.

In case there's no such circle, output -1. In case there are many, output -2.

Example(s)
sample input
sample output
3
0 0 5
7 16 10
10 0 5
5.0 5.0 5.0



The example is illustrated by the following picture:

<|response|>
## 1. Abridged problem statement

You are given **n** circles \((x_i, y_i, r_i)\) on the plane \((1 \le n \le 10^5)\).
Find a circle \((x, y, r)\) that is **orthogonal** to every given circle.

Output:
- the unique such circle's center and radius if it exists;
- `-1` if no such circle exists;
- `-2` if there are infinitely many.

Numbers are accepted within \(10^{-6}\) absolute/relative error.

---

## 2. Key observations

### Observation A: Orthogonality turns into a simple equation

Two circles with centers \(C=(x,y)\) and \(C_i=(x_i,y_i)\) and radii \(r\), \(r_i\) are orthogonal iff:
\[
|C - C_i|^2 = r^2 + r_i^2
\]
i.e.
\[
(x-x_i)^2 + (y-y_i)^2 = r^2 + r_i^2
\]

### Observation B: Linearize using a substitution

Expand:
\[
x^2 - 2xx_i + x_i^2 + y^2 - 2yy_i + y_i^2 = r^2 + r_i^2
\]
Rearrange:
\[
(-2x_i)x + (-2y_i)y + (x^2+y^2-r^2) = -(x_i^2+y_i^2-r_i^2)
\]

Let:
\[
z = x^2 + y^2 - r^2
\]
Then each given circle yields a **linear** equation in \((x,y,z)\):
\[
(-2x_i)x + (-2y_i)y + 1\cdot z = -(x_i^2+y_i^2-r_i^2)
\]

So we get a linear system with **3 variables** and up to **\(10^5\)** equations.

### Observation C: The system outcome directly maps to the required output

After solving the linear system:
- **Inconsistent** ⇒ no such circle ⇒ `-1`
- **Rank < 3** ⇒ infinitely many \((x,y,z)\) ⇒ infinitely many orthogonal circles ⇒ `-2`
- **Unique \((x,y,z)\)** ⇒ compute:
  \[
  r^2 = x^2 + y^2 - z
  \]
  If \(r^2 \le 0\) ⇒ no valid real circle radius ⇒ `-1`, else output \((x,y,\sqrt{r^2})\).

---

## 3. Full solution approach

1. **Read input** circles \((x_i,y_i,r_i)\).
2. For each circle, create one linear equation in \((x,y,z)\):
   - \(a=-2x_i\)
   - \(b=-2y_i\)
   - \(c=1\)
   - \(d=-(x_i^2 + y_i^2 - r_i^2)\)
   giving \(a x + b y + c z = d\).
3. Solve the overdetermined system using **Gauss-Jordan elimination** on an augmented matrix of size \(n \times 4\), but only **3 variable columns**, so runtime is essentially \(O(n)\) with a small constant.
   - Use partial pivoting for stability.
   - Use an epsilon (e.g. \(10^{-8}\)) to treat tiny values as zero.
4. After elimination:
   - If a row becomes \(0=nonzero\) ⇒ print `-1`.
   - If the rank (number of pivots) is < 3 ⇒ print `-2`.
   - Otherwise extract unique \((x,y,z)\), compute \(r^2=x^2+y^2-z\).
     - If \(r^2 \le \varepsilon\) ⇒ print `-1`
     - Else print \(x, y, \sqrt{r^2}\).

---

## 4. C++ Solution

```cpp
#include <bits/stdc++.h>

using namespace std;

template<typename T1, typename T2>
ostream& operator<<(ostream& out, const pair<T1, T2>& x) {
    return out << x.first << ' ' << x.second;
}

template<typename T1, typename T2>
istream& operator>>(istream& in, pair<T1, T2>& x) {
    return in >> x.first >> x.second;
}

template<typename T>
istream& operator>>(istream& in, vector<T>& a) {
    for(auto& x: a) {
        in >> x;
    }
    return in;
};

template<typename T>
ostream& operator<<(ostream& out, const vector<T>& a) {
    for(auto x: a) {
        out << x << ' ';
    }
    return out;
};

const long double eps = 1e-8l;

int n;
vector<array<int, 3>> circles;

void read() {
    cin >> n;
    circles.resize(n);
    for(int i = 0; i < n; i++) {
        cin >> circles[i][0] >> circles[i][1] >> circles[i][2];
    }
}

void solve() {
    // We only want to find 3 reals as the answer. Let's treat them as
    // variables, and each of the circles as constraints. We effectively want to
    // solve this system with 3 variables and N constraints. There are a few
    // possibilities - there is exactly one solution, no solutions, or
    // infinitely many. The only tricky bit is that the systems are non-linear,
    // so we can't directly use Gauss or another popular approach. But we will
    // start with defining the constraints and show that we can actually convert
    // this to a linear system.
    //
    // We can notice two circles are orthogonal if the square of the distance
    // between the centers equals the sum of the squares of the radii:
    //
    //     (xi-x)^2 + (yi-y)^2 = ri^2 + r^2
    //     -2*xi*x + -2*yi*y + x^2 + y^2 - r^2 = -xi^2 - yi^2 + ri^2
    //     A*x + B*y + (x^2 + y^2 - r^2) = C
    //
    // where (x, y, r) are our variables. Let's treat z = x^2 + y^2 - r^2,
    // making the system linear which we can solve with a simple Gauss in O(n)
    // as we have only 3 variables. We will have one of the below:
    //
    //     1) A unique solution for (x, y, z), in which case we set the radius
    //        r = sqrt(x^2 + y^2 - z).
    //
    //     2) No solutions in which case we can't really set r.
    //
    //     3) Infinitely many solution. Similarly to (1), any value given by
    //        r = sqrt(x^2 + y^2 - z) will satisfy the system.
    //
    // This means it's enough to solve the system after the transformation of
    // introducing z. Complexity as already mentioned is linear.

    vector<array<long double, 4>> matrix;
    for(auto& circ: circles) {
        int64_t xi = circ[0], yi = circ[1], ri = circ[2];
        long double a = -2.0l * xi;
        long double b = -2.0l * yi;
        long double c = 1.0l;
        long double d = -(xi * xi + yi * yi - ri * ri);
        matrix.push_back({a, b, c, d});
    }

    int rows = matrix.size();
    int cols = 3;

    int pivot_row = 0;
    for(int col = 0; col < cols; col++) {
        int best = -1;
        for(int row = pivot_row; row < rows; row++) {
            if(abs(matrix[row][col]) > eps) {
                if(best == -1 ||
                   abs(matrix[row][col]) > abs(matrix[best][col])) {
                    best = row;
                }
            }
        }

        if(best == -1) {
            continue;
        }

        swap(matrix[pivot_row], matrix[best]);

        long double divisor = matrix[pivot_row][col];
        for(int c = 0; c <= cols; c++) {
            matrix[pivot_row][c] /= divisor;
        }

        for(int row = 0; row < rows; row++) {
            if(row != pivot_row && abs(matrix[row][col]) > eps) {
                long double factor = matrix[row][col];
                for(int c = 0; c <= cols; c++) {
                    matrix[row][c] -= factor * matrix[pivot_row][c];
                }
            }
        }

        pivot_row++;
    }

    for(int row = pivot_row; row < rows; row++) {
        if(abs(matrix[row][3]) > eps) {
            cout << -1 << '\n';
            return;
        }
    }

    if(pivot_row < 3) {
        cout << -2 << '\n';
        return;
    }

    vector<long double> solution(3);
    for(int i = 0; i < 3; i++) {
        solution[i] = matrix[i][3];
    }

    long double x = solution[0];
    long double y = solution[1];
    long double z = solution[2];
    long double r_squared = x * x + y * y - z;
    if(r_squared <= eps) {
        cout << -1 << '\n';
        return;
    }

    long double r = sqrt(r_squared);
    cout << fixed << setprecision(10) << x << ' ' << y << ' ' << r << '\n';
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    int T = 1;
    // cin >> T;
    for(int test = 1; test <= T; test++) {
        read();
        solve();
    }

    return 0;
}
```

---

## 5. Python implementation with detailed comments

```python
import sys
import math

EPS = 1e-8

def solve() -> None:
    data = sys.stdin.buffer.read().split()
    if not data:
        return
    it = iter(data)
    n = int(next(it))

    # Build augmented matrix rows: [a, b, c, d] for equation
    # a*x + b*y + c*z = d, with z = x^2 + y^2 - r^2.
    mat = []
    for _ in range(n):
        xi = int(next(it)); yi = int(next(it)); ri = int(next(it))
        a = -2.0 * xi
        b = -2.0 * yi
        c = 1.0
        d = -(xi * xi + yi * yi - ri * ri)
        mat.append([a, b, c, float(d)])

    rows = n
    cols = 3
    pivot_row = 0

    # Gauss-Jordan elimination with partial pivoting
    for col in range(cols):
        best = -1
        best_abs = 0.0
        for r in range(pivot_row, rows):
            v = abs(mat[r][col])
            if v > EPS and v > best_abs:
                best_abs = v
                best = r

        if best == -1:
            continue  # free variable, no pivot in this column

        mat[pivot_row], mat[best] = mat[best], mat[pivot_row]

        div = mat[pivot_row][col]
        for c in range(cols + 1):
            mat[pivot_row][c] /= div

        for r in range(rows):
            if r == pivot_row:
                continue
            factor = mat[r][col]
            if abs(factor) <= EPS:
                continue
            for c in range(cols + 1):
                mat[r][c] -= factor * mat[pivot_row][c]

        pivot_row += 1

    # Inconsistency check: 0 = d (d must be ~0)
    for r in range(pivot_row, rows):
        if abs(mat[r][3]) > EPS:
            sys.stdout.write("-1\n")
            return

    # Rank < 3 => infinitely many solutions
    if pivot_row < 3:
        sys.stdout.write("-2\n")
        return

    # Unique solution (x,y,z)
    x = mat[0][3]
    y = mat[1][3]
    z = mat[2][3]

    r2 = x * x + y * y - z
    if r2 <= EPS:
        sys.stdout.write("-1\n")
        return

    r = math.sqrt(r2)
    sys.stdout.write(f"{x:.10f} {y:.10f} {r:.10f}\n")


if __name__ == "__main__":
    solve()
```
