## 1) Abridged problem statement

You are given the six edge lengths of a tetrahedron \(ABCD\) (in order: \(AB, AC, AD, BC, BD, CD\)), each a positive integer \(\le 1000\). Compute and print the tetrahedron’s volume with exactly 4 digits after the decimal point.

---

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

### Key idea
A tetrahedron is a 3D simplex. If you know all pairwise distances between its 4 vertices (i.e., its 6 edges), you can compute its volume directly using the **Cayley–Menger determinant**, which generalizes Heron’s formula to higher dimensions.

This avoids:
- reconstructing coordinates,
- computing a base triangle area + height,
- solving systems / vector geometry.

### Cayley–Menger determinant for a tetrahedron
Let the vertices be \(A,B,C,D\). Define the squared distances:
- \(d_{AB}^2, d_{AC}^2, d_{AD}^2, d_{BC}^2, d_{BD}^2, d_{CD}^2\)

Build the \(5 \times 5\) matrix:
\[
M =
\begin{pmatrix}
0 & 1 & 1 & 1 & 1 \\
1 & 0 & d_{AB}^2 & d_{AC}^2 & d_{AD}^2 \\
1 & d_{AB}^2 & 0 & d_{BC}^2 & d_{BD}^2 \\
1 & d_{AC}^2 & d_{BC}^2 & 0 & d_{CD}^2 \\
1 & d_{AD}^2 & d_{BD}^2 & d_{CD}^2 & 0
\end{pmatrix}
\]

Then:
\[
V^2 = \frac{\det(M)}{288}
\quad\Longrightarrow\quad
V = \sqrt{\frac{\det(M)}{288}}
\]

The constant \(288\) comes from the general simplex formula:
\[
V_n^2 = \frac{(-1)^{n+1}}{2^n (n!)^2}\det(M)
\]
For a tetrahedron, \(n=3\):
\[
2^3 (3!)^2 = 8 \cdot 36 = 288
\]
and the sign works out so that valid tetrahedra produce a non-negative value (small negative due to floating errors can happen).

### Computing the determinant
The solution computes \(\det(M)\) via **Gaussian elimination with partial pivoting**:
- For each column \(i\), choose a pivot row with the largest absolute value in that column.
- Swap rows if needed (flips sign of determinant).
- Multiply determinant by the pivot value.
- Eliminate entries below the pivot.

This is \(O(5^3)\), trivial in time.

### Numerical notes
- Uses `double`.
- If the determinant is slightly negative due to rounding, `abs` is taken before `sqrt`.
- Output formatted to 4 digits after the decimal.

---

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

double AB, AC, AD, BC, BD, CD;

void read() { cin >> AB >> AC >> AD >> BC >> BD >> CD; }

double gaussian_det(vector<vector<double>> mat) {
    int n = mat.size();
    double det = 1.0;
    for(int i = 0; i < n; ++i) {
        int pivot = i;
        for(int j = i + 1; j < n; ++j) {
            if(abs(mat[j][i]) > abs(mat[pivot][i])) {
                pivot = j;
            }
        }
        if(pivot != i) {
            swap(mat[i], mat[pivot]);
            det = -det;
        }
        if(abs(mat[i][i]) < 1e-10) {
            return 0.0;
        }
        det *= mat[i][i];
        for(int j = i + 1; j < n; ++j) {
            double factor = mat[j][i] / mat[i][i];
            for(int k = i; k < n; ++k) {
                mat[j][k] -= factor * mat[i][k];
            }
        }
    }
    return det;
}

void solve() {
    // We know that the volume of the pyramid is 1/3 * h * B, where B is the
    // area of the base and h is the height. The area of the base can easily be
    // calculated with Heron's formula, while the height we can bash / compute
    // using complex numbers or by defining a vector base and solving a system.
    // This approach is a bit tedious so we might want to consider alternatives.
    //
    // There is also an alternative approach using more theory. We can use
    // Cayley-Menger's determinant, which can be thought of as a generalization
    // of Heron's formula. It essentially gives us a formula for the volume of
    // the n-dimensional simplex, given by the pairwise distances between all
    // points. In particular, the formula says that the determinant of the
    // pairwise distance matrix with 1s added as the final row and column, times
    // (-1)^(n+1)/((n!)^2*2^n) is the square of the volume of the n-dimensinoal
    // simplex. This is precisely what we have here for n=4 and so we implement
    // this. More details can be found here:
    //
    //     https://en.wikipedia.org/wiki/Cayley%E2%80%93Menger_determinant

    double a2 = AB * AB, b2 = AC * AC, c2 = AD * AD;
    double d2 = BC * BC, e2 = BD * BD, f2 = CD * CD;
    vector<vector<double>> m(5, vector<double>(5));
    m[0] = {0, 1, 1, 1, 1};
    m[1] = {1, 0, a2, b2, c2};
    m[2] = {1, a2, 0, d2, e2};
    m[3] = {1, b2, d2, 0, f2};
    m[4] = {1, c2, e2, f2, 0};
    double determinant = gaussian_det(m);
    double volume = sqrt(abs(determinant) / 288.0);
    cout << fixed << setprecision(4) << volume << endl;
}

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

    int T = 1;
    // cin >> T;
    for(int test = 1; test <= T; test++) {
        read();
        // cout << "Case #" << test << ": ";
        solve();
    }

    return 0;
}
```

---

## 4) Python solution (same approach) with detailed comments

```python
import sys
import math

def determinant_gauss(a):
    """
    Compute determinant of a square matrix a (list of lists) using
    Gaussian elimination with partial pivoting.
    Works in O(n^3), which is tiny here (n=5).
    """
    n = len(a)
    # Make a deep copy so we don't mutate the original matrix
    mat = [row[:] for row in a]

    det = 1.0

    for i in range(n):
        # Choose pivot row with maximum absolute value in column i
        pivot = i
        for r in range(i + 1, n):
            if abs(mat[r][i]) > abs(mat[pivot][i]):
                pivot = r

        # If pivot is (near) zero, determinant is zero
        if abs(mat[pivot][i]) < 1e-12:
            return 0.0

        # Swap pivot row into place; each swap flips determinant sign
        if pivot != i:
            mat[i], mat[pivot] = mat[pivot], mat[i]
            det = -det

        pivot_val = mat[i][i]
        det *= pivot_val

        # Eliminate below pivot
        for r in range(i + 1, n):
            factor = mat[r][i] / pivot_val
            # Start from column i to avoid unnecessary work
            for c in range(i, n):
                mat[r][c] -= factor * mat[i][c]

    return det


def main():
    # Read six edge lengths
    data = sys.stdin.read().strip().split()
    AB, AC, AD, BC, BD, CD = map(float, data)

    # Squared lengths
    a2 = AB * AB
    b2 = AC * AC
    c2 = AD * AD
    d2 = BC * BC
    e2 = BD * BD
    f2 = CD * CD

    # Cayley–Menger 5x5 matrix for tetrahedron volume
    M = [
        [0.0, 1.0, 1.0, 1.0, 1.0],
        [1.0, 0.0, a2,  b2,  c2 ],
        [1.0, a2,  0.0, d2,  e2 ],
        [1.0, b2,  d2,  0.0, f2 ],
        [1.0, c2,  e2,  f2,  0.0],
    ]

    det = determinant_gauss(M)

    # V^2 = det/288, take abs to protect from tiny negative due to rounding
    volume = math.sqrt(abs(det) / 288.0)

    # Print with 4 digits after decimal
    print(f"{volume:.4f}")


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

---

## 5) Compressed editorial

Use the Cayley–Menger determinant to compute tetrahedron volume from its six edge lengths. Form the \(5\times5\) matrix with 0 on the diagonal, squared edge lengths in the distance positions, and a top row/left column of ones (except the top-left 0). Compute \(\det(M)\) via Gaussian elimination. Then:
\[
V=\sqrt{\frac{\det(M)}{288}}
\]
Print \(V\) to 4 decimals (use `abs(det)` to avoid negative zero due to floating errors).