## 1) Concise 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 (i.e., the circles intersect and their tangents at intersection points are perpendicular).

Output:
- the unique orthogonal circle’s center and radius (as reals) if it exists uniquely;
- `-1` if no such circle exists;
- `-2` if there are infinitely many.

Accuracy: \(10^{-6}\) absolute or relative error.

---

## 2) Detailed editorial (how the solution works)

### Key geometry fact: orthogonality condition becomes an equation
Two circles:
- unknown: center \((x, y)\), radius \(r\)
- given circle \(i\): center \((x_i, y_i)\), radius \(r_i\)

are orthogonal iff:
\[
(x - x_i)^2 + (y - y_i)^2 = r^2 + r_i^2
\]
(standard characterization: squared center distance equals sum of squared radii).

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)
\]

### Linearization trick: introduce a third variable
Let:
\[
z = x^2 + y^2 - r^2
\]
Then each circle produces 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:
- unknowns: \(x, y, z\) (3 variables)
- equations: \(n\) (up to \(10^5\))

### Solve the linear system with Gaussian elimination (3 variables)
We perform Gauss-Jordan elimination on an \(n \times 4\) augmented matrix:
\[
[A\;B\;C\mid D]
\]
where each row is:
- \(A=-2x_i\)
- \(B=-2y_i\)
- \(C=1\)
- \(D=-(x_i^2 + y_i^2 - r_i^2)\)

Because there are only **3 columns of variables**, elimination is effectively \(O(n)\) (constant small factor).

### Determine which of the three outcomes applies
After elimination:

1) **No solution**  
If we get a row like:
\[
0x + 0y + 0z = d,\;\; d\ne 0
\]
then inconsistent ⇒ output `-1`.

2) **Infinitely many solutions**  
If rank \(< 3\) (fewer than 3 pivots), there are infinitely many \((x,y,z)\) satisfying all constraints ⇒ output `-2`.
(And indeed infinitely many circles satisfy orthogonality constraints.)

3) **Unique solution for \((x,y,z)\)**  
Then compute \(r\) from:
\[
z = x^2 + y^2 - r^2 \Rightarrow r^2 = x^2 + y^2 - z
\]
- if \(r^2 \le 0\), no valid circle radius ⇒ output `-1`
- else \(r=\sqrt{r^2}\) and print \(x,y,r\)

### Why checking \(r^2 > 0\) is necessary
The linear system ensures algebraic orthogonality equations, but a real circle requires \(r>0\).  
If \(r^2\) is non-positive (within epsilon), there is no valid circle.

### Numerical notes
- Use `long double` for stability.
- Use an epsilon (`1e-8`) to treat near-zero as zero.
- Output with sufficient precision (10 decimals).

---

## 3) Provided C++ solution with detailed line-by-line comments

```cpp
#include <bits/stdc++.h>
using namespace std;

// Pretty-print a pair (not used in the final logic, but handy utilities)
template<typename T1, typename T2>
ostream& operator<<(ostream& out, const pair<T1, T2>& x) {
    return out << x.first << ' ' << x.second;
}

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

// Read a whole vector by reading each element
template<typename T>
istream& operator>>(istream& in, vector<T>& a) {
    for (auto& x: a) {
        in >> x;
    }
    return in;
};

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

// Epsilon for floating-point comparisons (treat values with abs <= eps as zero)
const long double eps = 1e-8l;

int n;
// Each circle stored as {xi, yi, ri} (integers as in input)
vector<array<int, 3>> circles;

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

void solve() {
    // We convert the circle orthogonality constraints into a linear system
    // in variables (x, y, z) where z = x^2 + y^2 - r^2.
    //
    // Each input circle gives:
    //   (-2*xi)*x + (-2*yi)*y + 1*z = -(xi^2 + yi^2 - ri^2)

    // Augmented matrix with rows [a b c d] representing:
    // a*x + b*y + c*z = d
    vector<array<long double, 4>> matrix;
    for (auto& circ: circles) {
        long long xi = circ[0], yi = circ[1], ri = circ[2];

        // Coefficients as derived above
        long double a = -2.0l * xi;
        long double b = -2.0l * yi;
        long double c = 1.0l;

        // Right-hand side: -(xi^2 + yi^2 - ri^2)
        long double d = -(xi * xi + yi * yi - ri * ri);

        matrix.push_back({a, b, c, d});
    }

    int rows = (int)matrix.size(); // n
    int cols = 3;                  // variables: x, y, z

    // Gauss-Jordan elimination:
    // pivot_row indicates where the next pivot should be placed
    int pivot_row = 0;

    // For each variable column, try to find a pivot and eliminate
    for (int col = 0; col < cols; col++) {
        int best = -1;

        // Find the row with the largest absolute coefficient in this column
        // among rows >= pivot_row (partial pivoting for stability).
        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;
                }
            }
        }

        // No pivot in this column => this variable is free (may lead to infinite solutions)
        if (best == -1) {
            continue;
        }

        // Move pivot row into position
        swap(matrix[pivot_row], matrix[best]);

        // Normalize pivot row so that pivot coefficient becomes 1
        long double divisor = matrix[pivot_row][col];
        for (int c = 0; c <= cols; c++) { // include RHS (col index 3)
            matrix[pivot_row][c] /= divisor;
        }

        // Eliminate this column from all other rows
        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];
                }
            }
        }

        // We successfully fixed one variable (one pivot)
        pivot_row++;
    }

    // Check for inconsistency:
    // rows below pivot_row should represent 0 = d; if d != 0 => impossible
    for (int row = pivot_row; row < rows; row++) {
        if (abs(matrix[row][3]) > eps) {
            cout << -1 << '\n';
            return;
        }
    }

    // If rank < 3 (fewer than 3 pivots), there are infinitely many solutions
    // in (x,y,z) => infinitely many orthogonal circles => output -2.
    if (pivot_row < 3) {
        cout << -2 << '\n';
        return;
    }

    // Now we have a unique solution; in reduced form, first 3 rows contain x,y,z
    vector<long double> solution(3);
    for (int i = 0; i < 3; i++) {
        solution[i] = matrix[i][3]; // RHS after elimination is the value of variable i
    }

    long double x = solution[0];
    long double y = solution[1];
    long double z = solution[2];

    // Recover r^2 from z = x^2 + y^2 - r^2  => r^2 = x^2 + y^2 - z
    long double r_squared = x * x + y * y - z;

    // r must be real and positive; if r^2 <= 0, no valid circle
    if (r_squared <= eps) {
        cout << -1 << '\n';
        return;
    }

    long double r = sqrt(r_squared);

    // Output with good precision
    cout << fixed << setprecision(10) << x << ' ' << y << ' ' << r << '\n';
}

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

    int T = 1;
    // cin >> T; // problem provides a single test
    for (int test = 1; test <= T; test++) {
        read();
        solve();
    }
    return 0;
}
```

---

## 4) Python solution (same approach) 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  # variables x,y,z
    pivot_row = 0

    # Gauss-Jordan elimination with partial pivoting.
    for col in range(cols):
        best = -1
        best_abs = 0.0

        # Choose pivot row with maximum absolute coefficient in this column.
        for r in range(pivot_row, rows):
            v = abs(mat[r][col])
            if v > EPS and v > best_abs:
                best_abs = v
                best = r

        # No pivot in this column -> variable remains free.
        if best == -1:
            continue

        # Swap pivot into place.
        mat[pivot_row], mat[best] = mat[best], mat[pivot_row]

        # Normalize pivot row so pivot coefficient becomes 1.
        div = mat[pivot_row][col]
        for c in range(cols + 1):  # include RHS
            mat[pivot_row][c] /= div

        # Eliminate this column from all other rows.
        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

    # Check inconsistency: 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

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

    # Extract solution. After Gauss-Jordan, first 3 pivot rows correspond to x,y,z.
    x = mat[0][3]
    y = mat[1][3]
    z = mat[2][3]

    # Recover radius from r^2 = x^2 + y^2 - z.
    r2 = x * x + y * y - z
    if r2 <= EPS:
        sys.stdout.write("-1\n")
        return

    r = math.sqrt(r2)

    # Print with sufficient precision.
    sys.stdout.write(f"{x:.10f} {y:.10f} {r:.10f}\n")


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

---

## 5) Compressed editorial

Orthogonality to circle \(i\) means:
\[
(x-x_i)^2 + (y-y_i)^2 = r^2 + r_i^2
\]
Expand and 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 circle gives a linear equation in \((x,y,z)\):
\[
(-2x_i)x + (-2y_i)y + z = -(x_i^2+y_i^2-r_i^2)
\]
Solve the \(n\)-equation, 3-variable linear system with Gauss-Jordan (only 3 columns ⇒ \(O(n)\)).  
- Inconsistent row ⇒ `-1`.  
- Rank < 3 ⇒ infinitely many solutions ⇒ `-2`.  
- Unique \((x,y,z)\) ⇒ compute \(r^2=x^2+y^2-z\); if \(r^2\le 0\) ⇒ `-1`, else output \((x,y,\sqrt{r^2})\).