6. МСГ

#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(vd g);
inline double sqr(const double& a) {
  return a * a;
}

inline double dot(const vd& a, const vd& b) {
  double ret = 0;
  for (int i = 0; i < a.size(); ++i)
    ret += a[i] * b[i];
  return ret;
}

inline double norm(const vd& a) {
  double ret = 0;
  for (int i = 0; i < a.size(); ++i)
    ret += sqr(a[i]);
  return sqrt(ret);
}

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

vd operator *(const vvd& a, const vd& b) {
  vd ret = vd(a.size(), 0.0);
  for (int i = 0; i < a.size(); ++i) {
    ret[i] = dot(a[i], b);
  }
  return ret;
}

int n;
int k = 0;
vvd A;
vd B;
vvd p;
vvd q;
vvd x;
vvd e;

void solve();

int main(){
  freopen("input.txt", "r", stdin);
  //freopen("output.txt", "w", stdout);
  n;
  cin >> n;
  A.resize(n, vd(n, 0.0));
  B.resize(n);
  p.resize(2, vd(n, 0.0));
  q.resize(2, vd(n, 0.0));
  x.resize(2, vd(n, 0.0));
  e.resize(2, vd(n, 0.0));

  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < n; ++j) {
      cin >> A[i][j];
    }
    cin >> B[i];
  }

  /*
  vd test = A * B;
  print(test);
  return 0;*/

  e[k] = B - A * x[k];
  p[k] = e[k];

  solve();
  cout << "solution:" << endl;
  for (int i = 0; i < n; ++i)
    cout << "x[" << i + 1 << "] " << setw(12) << fixed << setprecision(6) << x[k][i] << endl;
  cout << endl;
}

void solve() {
  /*for (int i = 0; i < n; ++i)
    cout << "x[" << i + 1 << "] " << setw(12) << fixed << setprecision(6) << x[k][i] << endl;
  cout << endl << endl;*/
  q[k] = A * p[k];
  double alpha = dot(e[k], p[k]) / dot(q[k], p[k]);
  x[1 - k] = x[k] + p[k] * alpha;
  e[1 - k] = e[k] - q[k] * alpha;

  double norm_e = norm(e[1 - k]);
  if (norm(e[1 - k]) < 1e-18)
    return;

  double beta = dot(e[1 - k], q[k]) / dot(p[k], q[k]);
  p[1 - k] = e[1 - k] - p[k] * beta;
  k = 1 - k;
  solve();
}

inline void print(vd g) {
  cout << "( ";
  for (int j = 0; j < n; ++j)
    cout << setw(12) << fixed << setprecision(6) << g[j];
  cout << endl;
}