#include <iostream>
#include <cmath>
#define N 995
#define A1 9
#define A23 -1
using namespace std;

void jacobi(double** A, double* b, double* x0, double* x1, double* residuum) {
	for (int i = 0; i < N; i++) {
		double sum1 = 0, sum2 = 0;
		for (int j = 0; j <= i - 1; j++) sum1 += A[i][j] * x0[j];
		for (int j = i + 1; j < N; j++) sum2 += A[i][j] * x0[j];
		x1[i] = (b[i] - sum1 - sum2) / A[i][i];
		x0[i] = x1[i];
	}
	for (int i = 0; i < N; i++) {
		for (int j = 0; j < N; j++) residuum[i] += A[i][j] * x1[j];
		residuum[i] -= b[i];
	}
}

int main() {
	double** A = new double*[N];
	for (int i = 0; i < N; i++) A[i] = new double[N];
	for (int i = 0; i < N; i++)
		for (int j = 0; j < N; j++) A[i][j] = 0.0;

	static double b[N], x0[N], x1[N], residuum[N], norm = 1; //x0 macierz xow poprzednich, x1 macierz xow nowych
	for (int i = 0; i < N; i++) {
		b[i] = sin(i * 10);
		for (int j = 0; j < N; j++) {
			if (i == j) A[i][j] = A1;
			else if (i == j + 1 || i == j - 1) A[i][j] = A23;
			else if (i == j + 2 || i == j - 2) A[i][j] = A23;
		}
	}
	int iter = 0;
	while (norm > pow(10, -6) && norm != 0.0) {
		double result = 0;
		jacobi(A, b, x0, x1, residuum);
		for (int i = 0; i < N; i++) result += pow(residuum[i], 2);
		norm = sqrt(result);
		iter++;
	}
	return 0;
}