Least square approximation Lera

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

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

using namespace std;

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

public:
    Matrix(int row, int column) {
        (*this).row = row;
        (*this).column = column;
        matrix = new double* [row];
        for (int i = 0; i < row; i++)
            matrix[i] = new double[column];
    }

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

    Matrix& operator=(const Matrix& right) {
        Matrix enmat(right.row, right.column);
        for (int i = 0; i < enmat.row; i++)
            for (int j = 0; j < enmat.column; j++)
                enmat.matrix[i][j] = right.matrix[i][j];
        (*this).matrix = enmat.matrix;
        (*this).row = enmat.row;
        (*this).column = enmat.column;
        return *this;
    }

    Matrix transp() {
        Matrix transp(column, row);
        for (int i = 0; i < row; i++)
            for (int j = 0; j < column; j++) {
                transp.matrix[j][i] = matrix[i][j];
            }
        return transp;
    }

    void permut(int row1, int row2, Matrix& a) {
        double* tempo = *(matrix + row1);
        *(matrix + row1) = *(matrix + row2);
        *(matrix + row2) = tempo;

        tempo = *(a.matrix + row1);
        *(a.matrix + row1) = *(a.matrix + row2);
        *(a.matrix + row2) = tempo;
    }

    void elimin(int topLine, int botLine, Matrix& a) {
        double f = matrix[botLine][topLine] / matrix[topLine][topLine];
        if (matrix[botLine][topLine] != 0) {
            for (int j = 0; j < column; j++) {
                matrix[botLine][j] -= matrix[topLine][j] * f;
                a.matrix[botLine][j] -= a.matrix[topLine][j] * f;
            }
        }
    }

    void norm(Matrix& a) {
        for (int i = 0; i < row; i++) {
            double ct = matrix[i][i];
            for (int j = 0; j < column; j++) {
                matrix[i][j] /= ct;
                a.matrix[i][j] /= ct;
            }
        }
    }

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

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

            if (maxLine != line)
                base.permut(line, maxLine, inverse);

            for (int i = line + 1; i < base.row; i++)
                base.elimin(line, i, inverse);
        }

        for (int line = row - 1; line > 0; line--)
            for (int i = line - 1; i >= 0; i--)
                base.elimin(line, i, inverse);
        base.norm(inverse);
        return inverse;
    }
    friend double* findCoef(double points[15][2]);
    friend ostream& operator<<(ostream& cout, Matrix matrix) {
        for (int i = 0; i < matrix.row; i++) {
            for (int j = 0; j < matrix.column; j++) {
                cout << matrix.matrix[i][j];
                if (j != matrix.column - 1)
                    cout << " ";
            }
            cout << endl;
        }
        return cout;
    }
};

double* findCoef(double points[15][2]) {
    Matrix a(15, 3), b(15, 1);

    for (int i = 0; i < 15; i++) {
        b.matrix[i][0] = points[i][1];
        a.matrix[i][0] = 1;
        a.matrix[i][1] = points[i][0];
        a.matrix[i][2] = points[i][0] * points[i][0];
    }
    Matrix at = a.transp();
    Matrix at_A = at * a;
    Matrix at_A_inverse = at_A.inv();
    Matrix at_b = at * b;
    Matrix res = at_A_inverse * at_b;
    double* coef = new double[3];
    for (int i = 0; i < 3; i++)
        coef[2 - i] = res.matrix[i][0];
    return coef;
}

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

    if (pipe != NULL) {
        fprintf(pipe, "set xrange [0:13]\nset yrange[0:17]\nset multiplot\n");

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

        double* coef = findCoef(points);

        fprintf(pipe, "set nokey\nplot '-' using 1:2 title ' ' with lines linecolor 3\n");

        for (int i = 0; i < 1300; i++) {
            double x = 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 7\n");
        for (int i = 0; i < 15; i++) {
            fprintf(pipe, "%f\t%f\n", points[i][0], points[i][1]);
        }
        fprintf(pipe, "e\n");

        fflush(pipe);

        _pclose(pipe);
    }
}