#include #include #include using namespace std; double eps = 1e-4; bool test (vector > &matrix) { for (int i = 0; i < matrix.size(); ++i) { for (int j = i + 1; j < matrix.size(); ++j) { if (fabs(matrix[i][j]) > eps) { return true; } } } return false; } vector > matrix_mult(vector >& x, vector >& y) { vector > r(x.size(), vector (y[0].size(), 0.0)); for (int i = 0; i < x.size(); i++) { for (int j = 0; j < y[0].size(); j++) { for (int k = 0; k < x[0].size(); k++) { r[i][j] += x[i][k] * y[k][j]; } } } return r; } int main() { freopen ("input.txt", "r", stdin); int n; cin >>n; vector > matrix (n, vector (n)); vector > t(n, vector (n, 0.0)); for (int i = 0; i < n; ++i) { t[i][i] = 1.0; } for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { cin >>matrix[i][j]; } } while (test(matrix)) { int maxi = 0, maxj = 1; double max_element = 0.0; for (int i = 0; i < n; ++i) { for (int j = i - 1; j >= 0; --j) { if (max_element < fabs(matrix[i][j])) { maxi = i; maxj = j; max_element = fabs(matrix[i][j]); } } } double p = 2.0 * matrix[maxi][maxj], q = matrix[maxi][maxi] - matrix[maxj][maxj]; double d = sqrt (p * p + q * q); double c, s; double tiny_eps = 1e-15; if (fabs(q) < tiny_eps) { c = 1.0 / sqrt(2.0); s = 1.0 / sqrt(2.0); } else { double r = fabs(q) / (2.0 * d); c = sqrt(0.5 + r); if (fabs(p) * 1000.0 < fabs(q)) { s = fabs(p) * (p * q > 0.0 ? 1.0 : -1.0) / (2.0 * c * d); } else { s = sqrt(0.5 - r) * (p * q > 0.0 ? 1.0 : -1.0); } } vector > cur_t(n, vector (n, 0.0)); for (int i = 0; i < n; ++i) { cur_t[i][i] = 1.0; } cur_t[maxi][maxi] = c; cur_t[maxj][maxj] = c; cur_t[maxi][maxj] = -s; cur_t[maxj][maxi] = s; t = matrix_mult(t, cur_t); double a1 = matrix[maxi][maxi], a2 = matrix[maxj][maxj]; matrix[maxi][maxi] = c * c * a1 + s * s * a2 + 2.0 * c * s * matrix[maxi][maxj]; matrix[maxj][maxj] = s * s * a1 + c * c * a2 - 2.0 * c * s * matrix[maxi][maxj]; matrix[maxi][maxj] = 0.0; matrix[maxj][maxi] = 0.0; for (int i = 0; i < n; ++i) { if (i == maxi || i == maxj) { continue; } a1 = matrix[i][maxi], a2 = matrix[i][maxj]; matrix[i][maxi] = c * a1 + s * a2; matrix[maxi][i] = c * a1 + s * a2; matrix[i][maxj] = -s * a1 + c * a2; matrix[maxj][i] = -s * a1 + c * a2; } } cout <