<|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.

162. Pyramids
time limit per test: 0.25 sec.
memory limit per test: 65536 KB
input: standard input
output: standard output



Recently in Farland, a country in Asia, a famous scientist Mr. Log Archeo has discovered ancient pyramids. But unlike those in Egypt and Central America, they have triangular (not rectangular) foundation. That is, they are tetrahedrons in mathematical sense. In order to find out some important facts about the early society of the country (it is widely believed that the pyramid sizes are in tight connection with Farland ancient calendar), Mr. Archeo needs to know the volume of the pyramids. Unluckily, he has reliable data about their edge lengths only. Please, help him!

Input
The file contains six positive integer numbers not exceeding 1000 separated by spaces, each number is one of the edge lengths of the pyramid ABCD. The order of the edges is the following: AB, AC, AD, BC, BD, CD.

Output
A real number - the volume printed accurate to four digits after decimal point.

Sample test(s)

Input
Sample input #1
1 1 1 1 1 1

Sample input #2
1000 1000 1000 3 4 5

Output
Sample output #1
0.1179

Sample output #2
1999.9937

Note
Author:	Nikolay Durov
Resource:	ACM ICPC 2002-2003 NEERC, Northern Subregion
Date:	November, 2002

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

You are given the six edge lengths of a tetrahedron \(ABCD\) in this order:

\[
AB,\ AC,\ AD,\ BC,\ BD,\ CD
\]

All are positive integers \(\le 1000\). Output the tetrahedron volume accurate to **4 digits after the decimal point**.

---

## 2) Key observations needed to solve the problem

1. A tetrahedron’s volume is uniquely determined by its **six pairwise edge lengths** (assuming the tetrahedron is geometrically valid).
2. Computing the volume via “base area + height” requires reconstructing coordinates or solving for a height—possible but more complicated.
3. The **Cayley–Menger determinant** directly gives the volume of a simplex (triangle, tetrahedron, etc.) using only pairwise distances.
4. For a tetrahedron, if \(M\) is the Cayley–Menger \(5\times 5\) matrix, then:
   \[
   V^2 = \frac{\det(M)}{288},\quad V=\sqrt{\frac{\det(M)}{288}}
   \]
5. \(\det(M)\) can be computed efficiently using **Gaussian elimination with partial pivoting** (matrix size is only 5, so it’s extremely fast).
6. Due to floating-point rounding, \(\det(M)\) might become a tiny negative number (like \(-1e{-10}\)) for a valid tetrahedron; using `abs(det)` before `sqrt` avoids NaNs.

---

## 3) Full solution approach based on the observations

### Step A: Read input and square edge lengths
Let:
- \(a=AB,\ b=AC,\ c=AD,\ d=BC,\ e=BD,\ f=CD\)
- Use squared lengths \(a^2,b^2,c^2,d^2,e^2,f^2\)

### Step B: Build Cayley–Menger matrix
Construct:

\[
M =
\begin{pmatrix}
0 & 1 & 1 & 1 & 1 \\
1 & 0 & a^2 & b^2 & c^2 \\
1 & a^2 & 0 & d^2 & e^2 \\
1 & b^2 & d^2 & 0 & f^2 \\
1 & c^2 & e^2 & f^2 & 0
\end{pmatrix}
\]

### Step C: Compute determinant \(\det(M)\)
Use Gaussian elimination:
- For each column \(i\), choose the pivot row with the largest absolute value in that column (partial pivoting).
- Swap rows if needed (each swap flips the sign of the determinant).
- Multiply determinant accumulator by the pivot.
- Eliminate values below the pivot.

This computes the determinant in \(O(5^3)\), trivial for the constraints.

### Step D: Convert determinant to volume
\[
V = \sqrt{\frac{|\det(M)|}{288}}
\]

### Step E: Print
Print with exactly **4 digits** after the decimal point.

---

## 4) C++ implementation

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

---

## 5) Python implementation with detailed comments

```python
import sys
import math

def determinant_gauss(mat):
    """
    Determinant via Gaussian elimination with partial pivoting.
    mat: list of lists (square matrix), will not be modified outside.

    Returns det(mat) as float.
    """
    a = [row[:] for row in mat]  # deep copy
    n = len(a)
    det = 1.0

    for i in range(n):
        # Find pivot row
        pivot = i
        for r in range(i + 1, n):
            if abs(a[r][i]) > abs(a[pivot][i]):
                pivot = r

        # Near-zero pivot => determinant is 0
        if abs(a[pivot][i]) < 1e-12:
            return 0.0

        # Swap to the top if needed; swap changes sign
        if pivot != i:
            a[i], a[pivot] = a[pivot], a[i]
            det = -det

        piv = a[i][i]
        det *= piv

        # Eliminate below
        for r in range(i + 1, n):
            factor = a[r][i] / piv
            for c in range(i, n):
                a[r][c] -= factor * a[i][c]

    return det

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

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

    # Cayley–Menger matrix
    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 avoid tiny negative due to rounding)
    volume = math.sqrt(abs(det) / 288.0)

    print(f"{volume:.4f}")

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

This method is robust, fast (tiny matrix), and uses only the given edge lengths—perfect for the problem constraints and required precision.