C++ least squares fit of N points, onto quadratic

least-squares-quadratic.cpp

void FitQuadratic(const std::vector<std::pair<double, double>>& xy, double& a, double& b, double& c) {
	// http://mathforum.org/library/drmath/view/72047.html

	/*
	Now all you have to do is take your data (your eight points) and
	evaluate the various sums

			n
	Sj0 = sum x_i^j
			i=1

	for j = 0 through 4, and

			n
	Sj1 = sum x_i^j*y_i
			i=1

	for j = 0 through 2.  Then you substitute into the formulas for a, b,
	and c, and you are done!	
	*/

	double S[5][2];
	for (size_t j = 0; j < 5; j++) {
		S[j][0] = 0;
		S[j][1] = 0;
		for (size_t i = 0; i < xy.size(); i++) {
			S[j][0] += pow(xy[i].first, (double) j);
			S[j][1] += pow(xy[i].first, (double) j) * xy[i].second; // only used when j = 0,1,2
		}
	}

	/*
	Now we can use Cramer's Rule to give a, b, and c as formulas in these
	Sjk values.  They all have the same denominator:

		(S00*S20*S40 - S10^2*S40 - S00*S30^2 + 2*S10*S20*S30 - S20^3)	
	*/

	double denom = S[0][0] * S[2][0] * S[4][0] - pow(S[1][0], 2) * S[4][0] - S[0][0] * pow(S[3][0], 2) + 2 * S[1][0] * S[2][0] * S[3][0] - pow(S[2][0], 3);

	a = S[0][1] * S[1][0] * S[3][0] - S[1][1] * S[0][0] * S[3][0] - S[0][1] * pow(S[2][0], 2) + S[1][1] * S[1][0] * S[2][0] + S[2][1] * S[0][0] * S[2][0] - S[2][1] * pow(S[1][0], 2);
	b = S[1][1] * S[0][0] * S[4][0] - S[0][1] * S[1][0] * S[4][0] + S[0][1] * S[2][0] * S[3][0] - S[2][1] * S[0][0] * S[3][0] - S[1][1] * pow(S[2][0], 2) + S[2][1] * S[1][0] * S[2][0];
	c = S[0][1] * S[2][0] * S[4][0] - S[1][1] * S[1][0] * S[4][0] - S[0][1] * pow(S[3][0], 2) + S[1][1] * S[2][0] * S[3][0] + S[2][1] * S[1][0] * S[3][0] - S[2][1] * pow(S[2][0], 2);

	a /= denom;
	b /= denom;
	c /= denom;
}