Jacobi method

#include <array>
#include <iomanip>
#include <iostream>
template <std::size_t Size>
static int jacobi(std::array<double, Size> &x,
                  const std::array<std::array<double, Size>, Size> &a,
                  const std::array<double, Size> &b, const double e,
                  int count) {
  ++count;
  std::array<double, Size> next;
  bool last = true;
  for (std::size_t i = 0; i < Size; ++i) {
    next[i] = b[i];
    for (std::size_t j = 0; j < Size; ++j) {
      if (i == j) { continue; }
      next[i] -= a[i][j] * x[j];
    }
    next[i] /= a[i][i];
    if (std::abs(next[i] - x[i]) >= e) { last = false; }
  }
  if (last) { return count; }
  x = next;
  return jacobi(x, a, b, e, count);
}
template <class T, std::size_t Size>
std::ostream &operator<<(std::ostream &out, const std::array<T, Size> &a) {
  out << '[';
  for (auto &&i = a.cbegin(); i != a.cend(); ++i) {
    out << (i == a.cbegin() ? "" : ", ") << *i;
  }
  out << ']';
  return out;
}
int main() {
  const int n = 5;
  std::array<std::array<double, n>, n> a = {{{9., 3., 1., 1., 1.},
                                             {1., 8., 4., 1., 1.},
                                             {1., 1., 4., 1., 1.},
                                             {1., 1., 4., 7., 1.},
                                             {1., 1., 4., 4., 8.}}};
  std::array<double, n> b = {-8., 18., -12., -15., 24.};
  // initialize zero vector
  std::array<double, n> x = {0., 0., 0., 0., 0.};
  const double e = 1.e-6;
  int count = jacobi(x, a, b, e, 0);
  std::cout << std::setprecision(6) << x << std::endl << count << std::endl;
}