Polynomial interpolation

#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_X(vector < vector <double> > & X, vector <double> & x){
    int N = X.size();

    for(int i = 0; i < N; ++i){
        for(int j = 0; j < N; ++j){
            X[i][j] = binary_power(x[i], N - j - 1);
        }
    }
}

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);
    vector <double> y(N);

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

    vector < vector <double> > X(N, vector <double> (N));
    fill_X(X, x);

    vector <double> a(N);

    triangulation(X, y);
    solve(X, a, y);

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

    system("PAUSE");
    return 0;
}