Journey to Jupiter

#include<bits/stdc++.h>
// #define int long long
#define f first
#define s second
using namespace std;

ifstream in("jupiter.in");
#define cin in

typedef long long ll;
const int MOD = 1e9 + 9;
const int inf = 1e9;
const double eps = 1e-9;
const double PI = acos(-1);


struct point3 {
    double x, y, z;
    // basic operations
    point3 operator-(point3 p) {
        return point3{x - p.x, y - p.y, z - p.z};
    }
    point3 operator/(double d) {
        return {x / d, y / d, z / d};
    }
    double operator^(point3 p) {
        return x * p.x + y * p.y + z * p.z;
    }
    double squareLen() {
        return (*this) ^ (*this);
    }
    double len() {
        return sqrt(squareLen());
    }
    point3 unitVec() {
        // the same vector but of length 1 (unit)
        return (*this) / len();
    }
};

int solveSystem(vector<vector<double>> &matrix, vector<double> &ans) {
    int n = matrix.size(), m = matrix[0].size() - 1;
    int row = 0;
    for(int col = 0; col < m; ++col) {
        int bestRow = row;
        for(int i = row; i < n; ++i) {
            if(abs(matrix[bestRow][col]) < abs(matrix[i][col])) {
                bestRow = i;
            }
        }
        if(bestRow >= n || abs(matrix[bestRow][col]) < eps) {
            continue;
        }
        for(int i = 0; i <= m; ++i) {
            swap(matrix[bestRow][i], matrix[row][i]);
        }

        // impart toata linia la matrix[row][col]
        for(int i = col + 1; i <= m; ++i) {
            matrix[row][i] = matrix[row][i] / matrix[row][col];
        }

        matrix[row][col] = 1.0;
        for(int i = row + 1; i < n; ++i) {
            for(int j = col + 1; j <= m; ++j) {
                matrix[i][j] -= matrix[i][col] * matrix[row][j];
            }
            matrix[i][col] = 0.0;
        }
        ++row;
    }

    bool infty = false;

    for(int i = n - 1; i >= 0; --i) {
        bool found = false;
        for(int j = 0; j <= m && !found; ++j) {
            if(abs(matrix[i][j]) > eps) {
                if(j == m) {
                    return 0;
                }
                ans[j] = matrix[i][m];
                for(int k = j + 1; k < m; ++k) {
                    ans[j] -= ans[k] * matrix[i][k];
                }
                found = true;
            }
        }
        if(!found) {
            infty = true;
        }
    }
    if(infty) {
        return 2;
    }
    return 1;
}

double nx, ny, nz, ax, ay, az, L, H;;

point3 solvePoint(point3 &A, point3 &normVec, double angle) {
    double S = L * L / 3.0 * sin(angle);
    double C = L * L / 3.0 * cos(angle);

    vector<vector<double>> systemMatrix(3, vector<double> (4, 0));

    systemMatrix[0][1] = -A.z, systemMatrix[0][2] = A.y, systemMatrix[0][3] = S * normVec.x;
    systemMatrix[1][0] = A.z, systemMatrix[1][2] = -A.x, systemMatrix[1][3] = S * normVec.y;
    systemMatrix[2][0] = A.x, systemMatrix[2][1] = A.y, systemMatrix[2][2] = A.z, systemMatrix[2][3] = C;

    vector<double> ans(3, 0);
    int verdict = solveSystem(systemMatrix, ans);
    assert(verdict == 1); // ecuatia are o singura solutie
    return point3{ans[0], ans[1], ans[2]};
}

void solve()
{
    cin >> nx >> ny >> nz >> ax >> ay >> az >> L >> H;

    point3 normVec = point3{nx, ny, nz}.unitVec();
    point3 A = {ax, ay, az};

    point3 B = solvePoint(A, normVec, 2.0 * PI / 3.0);
    point3 C = solvePoint(A, normVec, -2.0 * PI / 3.0);

    A = A - point3{0, 0, -H};
    B = B - point3{0, 0, -H};
    C = C - point3{0, 0, -H};

    cout << fixed << setprecision(10) << A.len() << ' ' << B.len() << ' ' << C.len() << '\n';
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(0);
    int T;
    cin >> T;
    while(T--) {
        solve();
    }
    return 0;
}