最小二乘法拟合二元n次曲线C++实现

/* 最小二乘法拟合二元多次曲线 */
/* 已知一堆(X, Y)数据的集合 */
class LeastSquareMethod
{
public:
	///<summary>
    ///用最小二乘法拟合二元多次曲线
    ///例如y=ax+b
    ///其中MultiLine将返回a，b两个参数。
    ///a对应MultiLine[1]
    ///b对应MultiLine[0]
    vector<double> MultiLine(const vector<double>& arrX, const vector<double>& arrY, int dimension)
    {
        int n = dimension + 1;//dimension次方程需要求 dimension+1个 系数
		vector<vector<double>> vGuass(n, vector<double>(n + 1));//高斯矩阵 例如：y=a0+a1*x+a2*x*x

		for (int i = 0; i < n; i++)
		{
		    int j;
			for (j = 0; j < n; j++)
			{
			    vGuass[i][j] = SumArr(arrX, j + i);
			}
			vGuass[i][j] = SumArr(arrX, i, arrY, 1);
		}
		return ComputeGuass(vGuass, n);
    }
	//求数组元素的n次方的和
	double SumArr(const vector<double>& arr, int n) const
    {
        double dSum = 0;
		for (int i = 0; i < arr.size(); i++)
		{
		    if (arr[i] != 0 || n != 0)
		    {
		        dSum += pow(arr[i], n);
		    }
			else
			{
			    dSum += 1;
			}
		}
		return dSum;
    }

	double SumArr(const vector<double>& arr1, int n1, const vector<double>& arr2, int n2) const
    {
        double dSum = 0;
		for (int i = 0; i < arr1.size(); i++)
		{
		    if ((arr1[i] != 0 || n1 != 0) && (arr2[i] != 0 || n2 != 0))
		    {
		        dSum += pow(arr1[i], n1) * pow(arr2[i], n2);
		    }
			else
				dSum += 1;
		}
		return dSum;
    }

	vector<double> ComputeGuass(vector<vector<double>>& Guass, int n)
    {
        int i, j;
		int k, m;
		double temp;
		double max;
		double s;
		vector<double> x(n, 0.0);

		for (j = 0; j < n; j++)
		{
		    max = 0;
			k = j;
			for (i = j; i < n; i++)
			{
			    if (abs(Guass[i][j]) > max))
			    {
			        max = Guass[i][j];
					k = i;
			    }
			}
			if (k != j)
			{
			    for (m = j; m < n + 1; m++)
			    {
			        temp = Guass[j][m];
					Guass[j][m] = Guass[k][m];
					Guass[k][m] = temp;
			    }
			}

			if (0 == max)
			{
				// "此线性方程为奇异线性方程"
			    return x;
			}
			for (i = j + 1; i < n; i++)
			{
			    s  = Guass[i][j];
				for (m = j; m < n + 1; m++)
				{
				    Guass[i][m] = Guass[i][m] - Guass[j][m] * s / (Guass[j][j]);
				}
			}
		}

		for (i = n - 1; i >= 0; i--)
		{
		    s = 0;
			for (j = i + 1; j < n; j++)
			{
			    s = s + Guass[i][j] * x[j];
			}
			x[i] = (Guass[i][n] - s) / Guass[i][i];
		}
		return x;//返回值是函数的系数
    }
};

数学原理完全不懂，根据C#代码实现的C++版本。