#include <Eigen/QR>

#include <vector>
#include <iostream>


void polyfit(const std::vector<double> &xv, const std::vector<double> &yv, std::vector<double> &coeff, int order)
{
	Eigen::MatrixXd A(xv.size(), order + 1);
	Eigen::VectorXd yv_mapped = Eigen::VectorXd::Map(&yv.front(), yv.size());
	Eigen::VectorXd result;

	assert(xv.size() == yv.size());
	assert(xv.size() >= order + 1);

	for (size_t i = 0; i < xv.size(); i++)
		for (size_t j = 0; j < order + 1; j++)
			A(i, j) = pow(xv.at(i), j);

	result = A.householderQr().solve(yv_mapped);

	coeff.resize(order + 1);
	for (size_t i = 0; i < order + 1; i++)
		coeff[i] = result[i];
}

int main()
{
	std::vector<double> x_values, y_values, coeff;

	x_values = {0.0, 1.0, 2.0, 3.0, 4,  5,   6, 7    };
	y_values = {1,   -0.5, 2.0, -1, 2, -3,  2 , -3  };

	polyfit(x_values, y_values, coeff, x_values.size()-1);
	for (int i = 0; i <= x_values.size()-1; i++){
		if (coeff[i] > 0)
			std::cout << "+";
		std::cout << coeff[i] << "*x^"<< i << " ";
	}

	system("pause");
	return 0;
}