Least square approximation stepic my

#include <iostream>
#include <iomanip>
#include <math.h>

using namespace std;

struct Point {
	int x, y;
};

class Matrix {
private:
	int row, col;
	double** matrix;

public:
	Matrix(int n, int m) {
		(*this).row = n;
		(*this).col = m;
		matrix = new double* [n];
		for (int i = 0; i < n; i++) {
			matrix[i] = new double[m];
		}
	}
	void set(int i, int j, double res) {
		matrix[i][j] = res;
	}

	Matrix transpose() {
		Matrix tr((*this).col, (*this).row);
		for (int i = 0; i < (*this).col; i++)
			for (int j = 0; j < (*this).row; j++)
				tr.matrix[i][j] = (*this).matrix[j][i];
		return tr;
	}

	Matrix inverse() {
		Matrix now(this->row, this->col), res(this->row, this->col);
		for (int i = 0; i < row; i++) {
			for (int j = 0; j < col; j++) {
				now.matrix[i][j] = matrix[i][j];
				res.matrix[i][j] = i == j;
			}
		}

		for (int line = 0; line < row; line++) {
			int maxLine = line;
			for (int i = line + 1; i < now.row; i++) {
				if (abs(now.matrix[i][line]) > abs(now.matrix[maxLine][line]))
					maxLine = i;
			}

			if (maxLine != line)
				now.permutation(line, maxLine, res);

			for (int i = line + 1; i < now.row; i++)
				now.elimination(line, i, res);
		}

		for (int line = row - 1; line > 0; line--)
			for (int i = line - 1; i >= 0; i--)
				now.elimination(line, i, res);
		now.normalization(res);
		return res;
	}

	friend Matrix operator*(const Matrix& a, const Matrix& b) {
		Matrix res(a.row, b.col);
		for (int i = 0; i < res.row; i++) {
			for (int j = 0; j < res.col; j++) {
				res.matrix[i][j] = 0;
				for (int k = 0; k < a.col; k++)
					res.matrix[i][j] += a.matrix[i][k] * b.matrix[k][j];
			}
		}
		return res;
	}

	friend ostream& operator<<(ostream& cout, const Matrix& base) {
		for (int i = 0; i < base.row; i++) {
			for (int j = 0; j < base.col; j++) {
				cout << fixed << setprecision(2) << base.matrix[i][j] + 1e-9;

				if (j != base.col - 1)
					cout << " ";
			}
			cout << endl;
		}
		return cout;
	}

	void permutation(int r1, int r2, Matrix& inverse) {
		double* temp = *(matrix + r1);
		*(matrix + r1) = *(matrix + r2);
		*(matrix + r2) = temp;

		temp = *(inverse.matrix + r1);
		*(inverse.matrix + r1) = *(inverse.matrix + r2);
		*(inverse.matrix + r2) = temp;
	}

	void elimination(int upperLine, int bottomLine, Matrix& inverse) {
		double constant = matrix[bottomLine][upperLine] / matrix[upperLine][upperLine];

		if (matrix[bottomLine][upperLine] != 0) {

			for (int j = 0; j < col; j++) {
				matrix[bottomLine][j] -= matrix[upperLine][j] * constant;
				inverse.matrix[bottomLine][j] -= inverse.matrix[upperLine][j] * constant;
			}
		}
	}

	void normalization(Matrix& inverse) {
		for (int i = 0; i < row; i++) {
			double constant = matrix[i][i];
			for (int j = 0; j < col; j++) {
				matrix[i][j] /= constant;
				inverse.matrix[i][j] /= constant;
			}
		}
	}
};

void approximation() {
	int n, degree;
	cin >> n;
	Point* points = new Point[n];
	for (int i = 0; i < n; i++)
		cin >> points[i].x >> points[i].y;
	cin >> degree;
	Matrix a(n, degree + 1), b(n, 1);

	for (int i = 0; i < n; i++) {
		b.set(i, 0, points[i].y);

		for (int j = 0; j < degree + 1; j++)
			a.set(i, j, pow(points[i].x, j));
	}
	cout << "A:\n" << a;

	Matrix at = a.transpose();

	Matrix at_A = at * a;
	cout << "A_T*A:\n" << at_A;

	Matrix at_A_inverse = at_A.inverse();
	cout << "(A_T*A)^-1:\n" << at_A_inverse;

	Matrix at_b = at * b;
	cout << "A_T*b:\n" << at_b;
	Matrix res = at_A_inverse * at_b;
	cout << "x~:\n" << res;
}

int main()
{
	approximation();
}