C++实现多项式拟合求系数

#include <iostream>
#include <vector>
#include <cmath>
using namespace std;
void polyfit(const vector<double>& x, const vector<double>& y, int n, vector<double>& a)
{
	int len = x.size();
	int i, j, k;
	vector<double> b(n + 1, 0);
	vector<vector<double>> A(n + 1, vector<double>(n + 1, 0));
	for (i = 0; i < len; i++)
	{
		for (j = 0; j <= n; j++)
		{
			b[j] += pow(x[i], j) * y[i];
			for (k = 0; k <= n; k++) 
		    {
				A[j][k] += pow(x[i], j + k);
			}
		}
	}
	for (i = 0; i <= n; i++)
	{
		for (j = i + 1; j <= n; j++)
		{
			double ratio = A[j][i] / A[i][i];
			for (k = i; k <= n; k++)
			{
				A[j][k] -= ratio * A[i][k];
			}
			b[j] -= ratio * b[i];
		}
	}
	for (i = n; i >= 0; i--)
	{
		double sum = 0;
		for (j = i + 1; j <= n; j++)
		{
			sum += a[j] * A[i][j];
		}
		a[i] = (b[i] - sum) / A[i][i];
	}
}
int main()
{
    vector<double> x;
	x.push_back(0.23);
	x.push_back(0.72);
	x.push_back(0.92);
	vector<double> y;
	y.push_back(0.23);
	y.push_back(0.665);
	y.push_back(0.92);
	vector<double> a(3);
	polyfit(x, y, 2, a);
	cout << "拟合多项式为：y = " << a[2] << "x^2 + " << a[1] << "x + " << a[0] << endl; 
    return 0;
}