#include <array>

template<typename T, size_t N>
std::array<T,N> newton(const std::array<T,N> &xs, const std::array<T,N> &ys)
{
	std::array<T,N> as = ys;
	for (size_t k = 1; k < N; k++) {
		for (size_t j = N-1; j > k-1; j--) {
			as[j] = (as[j] - as[j-1]) / (xs[j] - xs[j-k]);
		}
	}
	return as;
}

int main()
{
	std::array<double, 4> xs = { 0.0, 1.0, -1.0, 2.0  }; 
	std::array<double, 4> ys = { 1.0, 5.0, -1.0, 17.0 }; // 1 + 2x + x^2 + x^3
	std::array<double, 4> as = newton(xs,ys);            // 1 + 4x + 1x(x-1) + 1x(x-1)(x+1)
	printf("%g %g %g %g\n", as[0], as[1], as[2], as[3]);
	return 0;
}

I2luY2x1ZGUgPGFycmF5PgoKdGVtcGxhdGU8dHlwZW5hbWUgVCwgc2l6ZV90IE4+CnN0ZDo6YXJyYXk8VCxOPiBuZXd0b24oY29uc3Qgc3RkOjphcnJheTxULE4+ICZ4cywgY29uc3Qgc3RkOjphcnJheTxULE4+ICZ5cykKewoJc3RkOjphcnJheTxULE4+IGFzID0geXM7Cglmb3IgKHNpemVfdCBrID0gMTsgayA8IE47IGsrKykgewoJCWZvciAoc2l6ZV90IGogPSBOLTE7IGogPiBrLTE7IGotLSkgewoJCQlhc1tqXSA9IChhc1tqXSAtIGFzW2otMV0pIC8gKHhzW2pdIC0geHNbai1rXSk7CgkJfQoJfQoJcmV0dXJuIGFzOwp9CgppbnQgbWFpbigpCnsKCXN0ZDo6YXJyYXk8ZG91YmxlLCA0PiB4cyA9IHsgMC4wLCAxLjAsIC0xLjAsIDIuMCAgfTsgCglzdGQ6OmFycmF5PGRvdWJsZSwgND4geXMgPSB7IDEuMCwgNS4wLCAtMS4wLCAxNy4wIH07IC8vIDEgKyAyeCArIHheMiArIHheMwoJc3RkOjphcnJheTxkb3VibGUsIDQ+IGFzID0gbmV3dG9uKHhzLHlzKTsgICAgICAgICAgICAvLyAxICsgNHggKyAxeCh4LTEpICsgMXgoeC0xKSh4KzEpCglwcmludGYoIiVnICVnICVnICVnXG4iLCBhc1swXSwgYXNbMV0sIGFzWzJdLCBhc1szXSk7CglyZXR1cm4gMDsKfQ==