#include <cmath>
#include <algorithm>
#include <iostream>
#include <iomanip>
#include <vector>
#include <stdexcept>

using namespace std;

using DVector = vector<double>;
using DMatrix = vector<vector<double>>;

class Cholesky {
  private:
    DMatrix L;

  public:
    Cholesky (const DMatrix& matrix, double eps = 1e-10)
      : L(matrix.size(), vector<double>(matrix.size())) {
      unsigned int n = L.size();
      for (unsigned int i = 0; i < n; i++)
        if (matrix[i].size() != n) {
          ostringstream sstr;
          sstr << "Matrix nicht quadratisch an Zeile " << i;
          throw invalid_argument{sstr.str()};
        }

      double trace = 0;
      for (unsigned int i = 0; i < n; i++)
        trace += matrix[i][i];

      for (unsigned int j = 0; j < n; j++) {
        L[j][j] = matrix[j][j];
        for (unsigned int k = 0; k < j; k++)
          L[j][j] -= L[j][k] * L[j][k];

        if (abs(L[j][j]) <= abs(eps * trace)) {
          ostringstream sstr;
          sstr << "Diagonalelement zu klein bei " << j;
          throw runtime_error{sstr.str()};
        }
        L[j][j] = sqrt(L[j][j]);

        for (unsigned int i = j + 1; i < n; i++) {
          L[i][j] = matrix[i][j];
          for (unsigned int k = 0; k < j; k++)
            L[i][j] -= L[j][k]*L[i][k];
          L[i][j] /= L[j][j];
        }
      }
    }

    double getEntry(unsigned int i, unsigned int j) const {
      double res = 0;
      for (unsigned int k = 0; k <= min(i, j); k++) {
        res += L[j][k] * L[i][k];
      }
      return res;
    }

    friend ostream& operator<<(ostream& stream, const Cholesky& c) {
      for (unsigned int i = 0; i < c.L.size(); i++) {
        stream << "(";
        for (unsigned int j = 0; j < c.L.size(); j++)
          stream << setw(9) << c.getEntry(i, j);
        stream << ")";
        if (i + 1 < c.L.size())
          stream << endl;
      }
      return stream;
    }

    DVector solve(const DVector& b) {
      DVector y(L.size()), x(L.size());
      for (unsigned int i = 0; i < L.size(); i++) {
        y[i] = b[i];
        for (unsigned int k = 0; k < i; k++)
          y[i] -= L[i][k] * y[k];
        y[i] /= L[i][i];
      }
      for (unsigned int m = 0; m < L.size(); m++) {
        const unsigned int i = L.size() - m - 1;
        x[i] = y[i];
        for (unsigned int k = i + 1; k < L.size(); k++)
          x[i] -= L[k][i] * x[k];
        x[i] /= L[i][i];
      }
      return x;
    }
};

ostream& operator<<(ostream& stream, const DVector& v) {
  stream << "(";
  for (unsigned int i = 0; i < v.size(); i++)
    stream << setw(9) << v[i];
  stream << ")";
  return stream;
}

int main() {
  unsigned int n;
  cout << "n: ";
  cin >> n;

  DMatrix A(n, vector<double>(n));
  for (unsigned i = 0; i < n; i++) {
    cout << "A[" << i << "][0.." << n - 1 << "]: ";
    for (double& x: A[i])
      cin >> x;
  }

  DVector b(n);
  cout << "b[0.." << n - 1 << "]: ";
  for (double& x: b)
    cin >> x;

  cout << endl;

  Cholesky L{A};
  cout << L << endl
       << endl;
  cout << b << endl;

  cout << endl;

  cout << L.solve(b) << endl;
}