Least square approximation my

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

#define GNUPLOT_NAME "C:\\gnuplot\\bin\\gnuplot.exe -persist"

using namespace std;

struct Point {
    double x, y;
    Point(int x, int y) : x(x), y(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;
    }

    double getElement(int i, int j) {
        return matrix[i][j];
    }

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

double* approximation(Point* points) {
    Matrix a(16, 3), b(16, 1);

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

        for (int j = 0; j < 3; j++)
            a.set(i, j, pow(points[i].x, j));
    }

    Matrix at = a.transpose();
    Matrix at_A = at * a;
    Matrix at_A_inverse = at_A.inverse();
    Matrix at_b = at * b;
    Matrix res = at_A_inverse * at_b;
    double* out = new double[3];
    for (int i = 0; i < 3; i++)
        out[2 - i] = res.getElement(i, 0);
    return out;
}

int main(){
    FILE* pipe = _popen(GNUPLOT_NAME, "w");

    if (pipe != NULL) {
        fprintf(pipe, "set xrange [-3:19]\nset yrange[4:15]\nset multiplot\n");

        Point points[16] = { {1, 6}, {2, 8}, {3, 11},
        {5, 12}, {6, 11}, {7, 11},
        {8, 12}, {9, 12}, {10, 11},
        {11, 12}, {12, 12}, {12, 10},
        {13, 9}, {14, 7}, {15, 10}, {16, 9}
        };

        double* coef = approximation(points);

        fprintf(pipe, "plot '-' using 1:2 title ' %.2f * x^2 %s %.2f %s %.2f ' with lines linecolor 2\n", coef[0],
            coef[1] > 0 ? "+" : "-", abs(coef[1]),
            coef[1] > 0 ? "* x +" : "* x -", coef[2]);

        for (int i = 0; i < 2200 + 1; i++) {
            double x = -3 + i / 100.0;
            double y = coef[0] * pow(x, 2) + coef[1] * x + coef[2];
            fprintf(pipe, "%f\t%f\n", x, y);
        }
        fprintf(pipe, "e\n");

        fprintf(pipe, "set nokey\nplot '-' using 1:2 title ' ' with points pointtype 5 pointsize 2\n");
        for (int i = 0; i < 16; i++) {
            fprintf(pipe, "%f\t%f\n", points[i].x, points[i].y);
        }
        fprintf(pipe, "e\n");

        fflush(pipe);

        _pclose(pipe);
    }
}