## 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 (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

### Key geometry fact: orthogonality condition becomes a linear 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
\]

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 with rows:
- \(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)\).

### Determine which of the three outcomes applies

After elimination:

1. **No solution**: if a row becomes \(0x + 0y + 0z = d\) with \(d \ne 0\), output `-1`.

2. **Infinitely many solutions**: if rank \(< 3\) (fewer than 3 pivots), output `-2`.

3. **Unique solution for \((x,y,z)\)**: compute \(r^2 = x^2 + y^2 - z\).
   - if \(r^2 \le 0\), output `-1`
   - else \(r=\sqrt{r^2}\) and print \(x,y,r\)

### 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. 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;
}
```

---

## 4. Python solution 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})\).
