6. метод вращений

 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
#define _CRT_SECURE_NO_WARNINGS
#include <iostream>
#include <iomanip>
#include <cmath>
#include <algorithm>
#include <vector>
#define vd vector<double>
#define vvd vector<vector<double> >
#define ui unsigned int
using namespace std;
inline void print();
inline double sqr(const double& a) {
return a * a;
}
inline double coeff(const double& a, const double& b, const double& c) {
return a / sqrt(sqr(b) + sqr(c));
}
vd operator +(const vd& a, const vd& b) {
vd ret = vd(a.size());
for (int i = 0; i < a.size(); ++i)
ret[i] = a[i] + b[i];
return ret;
}
vd operator *(const vd& a, double d) {
vd ret = vd(a.size());
for (int i = 0; i < a.size(); ++i)
ret[i] = a[i] * d;
return ret;
}
int n;
vvd g;
vd ans;
inline void rotate(const int& ind) {
for (int i = ind + 1; i < n; ++i) {
double c = coeff(g[ind][ind], g[ind][ind], g[i][ind]);
double s = coeff(g[i][ind], g[ind][ind], g[i][ind]);
vd old_ind = g[ind];
vd old_i = g[i];
g[ind] = old_ind * c + old_i * s;
g[i] = old_ind * (-s) + old_i * c;
}
}
inline void triangle_solve() {
int k = n - 1;
ans.resize(n);
for (int i = k; i >= 0; --i) {
ans[i] = g[i][n] / g[i][i];
for (int j = i - 1; j >= 0; --j) {
g[j][n] -= ans[i] * g[j][i];
g[j][i] = 0;
}
}
cout << "solution:" << endl;
for (int i = 0; i < n; ++i)
cout << "x[" << i + 1 << "] " << setw(12) << fixed << setprecision(6) << ans[i] << endl;
}
int main(){
//freopen("input.txt", "r", stdin);
//freopen("output.txt", "w", stdout);
cin >> n;
g.resize(n, vd(n + 1));
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n + 1; ++j) {
cin >> g[i][j];
}
}
print();
for (int i = 0; i < n - 1; ++i) {
rotate(i);
print();
}
triangle_solve();
return 0;
}
inline void print() {
for (int i = 0; i < n; ++i) {
cout << "( ";
for (int j = 0; j < n; ++j)
cout << setw(12) << fixed << setprecision(6) << g[i][j];
cout << " ) " << fixed << setprecision(6) << g[i][n];
cout << endl;
}
cout << endl;
}