Least squares

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

using namespace std;

double binary_power(double a, int n){
    if(n == 0){
        return 1;
    }

    if(n % 2 == 0){
        double res = binary_power(a, n / 2);

        return res * res;
    }
    return a * binary_power(a, n - 1);
}

void fill(vector < vector <double> > & M, vector <double> & x, vector <double> & y, vector <double> & b){
    int N = x.size();
    int k = M.size() - 1;
    vector <double> Sum(2 * k + 1);

    Sum[0] = N;
    for(int i = 1; i <= 2 * k; ++i){
        for(int j = 0; j < N; ++j){
            Sum[i] += binary_power(x[j], i);
        }
    }

    for(int i = 0; i <= k; ++i){
        for(int j = 0; j <= k; ++j){
            M[i][j] = Sum[2 * k - (i + j)];
        }
    }

    for(int i = 0; i <= k; ++i){
        for(int j = 0; j < N; ++j){
            b[i] += binary_power(x[j], k - i) * y[j];
        }
    }
}

void triangulation(vector < vector <double> > & A, vector <double> & b){
    int N = A.size();

    for(int j = 0; j < N; ++j){
        if(A[j][j] == 0){
            int index = 0;
            for(int i = 1; i < N; ++i){
                if(fabs(A[i][j]) > A[index][j] && A[i][j] != 0){
                    index = i;
                }
            }
            swap(A[j], A[index]);
            swap(b[j], b[index]);
        }
        for(int i = j + 1; i < N; ++i){
            double m = A[i][j] / A[j][j];
            for(int k = j; k < N; ++k){
                A[i][k] -= m * A[j][k];
            }
            b[i] -= m * b[j];
        }
    }
}

void solve(vector < vector <double> > & A, vector <double> & x, vector <double> & b){
    int N = A.size();

    for(int i = N - 1; i >= 0; --i){
        double Sum = 0;
        for(int j = i + 1; j < N; ++j){
            Sum += A[i][j] * x[j];
        }
        x[i] = (b[i] - Sum) / A[i][i];
    }
}

int main(){
    cout << "Enter the number of points:\n";
    int N;
    cin >> N;

    vector <double> x(N), y(N);

    cout << "Enter " << N << " pairs of coordinates (x, y):\n";
    for(int i = 0; i < N; ++i){
        cin >> x[i] >> y[i];
    }

    cout << "Enter polynomial degree:\n";
    int k;
    cin >> k;

    vector < vector <double> > M(k + 1, vector <double> (k + 1));
    vector <double> b(k + 1);
    fill(M, x, y, b);

    vector <double> a(N);

    triangulation(M, b);
    solve(M, a, b);

    cout << "Polynomial coefficients (starting from the highest power):\n";
    for(int i = 0; i <= k; ++i){
        cout << a[i] << " ";
    }
    cout << "\n";

    system("PAUSE");
    return 0;
}