include iostream include vector include iomanip using namespace std do

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
#include <iostream>
#include <vector>
#include <iomanip>
using namespace std;
double eps = 1e-4;
bool test (vector <vector <double> > &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 <vector <double> > matrix_mult(vector <vector <double> >& x, vector <vector <double> >& y) {
vector <vector <double> > r(x.size(), vector <double>(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 <vector <double> > matrix (n, vector <double> (n));
vector <vector <double> > t(n, vector <double> (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 <vector <double> > cur_t(n, vector <double> (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 <<fixed <<setprecision(5);
for (int i = 0; i < n; ++i) {
cout <<matrix[i][i] <<endl;
}
cout <<endl;
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
cout <<t[i][j] <<" ";
}
cout <<endl;
}
return 0;
}