基于最小二乘法的C++二次抛物线拟合

//
//a0x^2+a1x+a2 = y
//r:相关系数^2
bool Parabola(double* x, double* y, int n, double* a, double* r)
{
const int CURVE_INDEX = 3;
double B[CURVE_INDEX];
double **A =new double*[CURVE_INDEX];
for (int i=0; i<CURVE_INDEX; i++)
{
A[i] = new double[CURVE_INDEX];
}
double sumX=0, sumXX=0, sumXXX=0, sumXXXX=0, sumY=0, sumXY=0, sumXXY=0;
for (int i=0; i<n; i++)
{
sumX += x[i];
sumXX += x[i]*x[i];
sumXXX += x[i]*x[i]*x[i];
sumXXXX += x[i]*x[i]*x[i]*x[i];
sumY += y[i];
sumXY+= x[i]*y[i];
sumXXY += x[i]*x[i]*y[i];
}
A[0][0] = sumXXXX;A[0][1]=sumXXX;A[0][2]=sumXX;
A[1][0] = sumXXX;A[1][1]=sumXX;A[1][2]=sumX;
A[2][0] = sumXX;A[2][1]=sumX;A[2][2]=n;
//
B[0] = sumXXY; B[1] = sumXY; B[2] = sumY;

bool solved = Solve(A, B, a, CURVE_INDEX);
for (int i=0; i<CURVE_INDEX; i++)
{
delete A[i];
}
delete A;
//
double *ny = new double
;
double meanNY=0, meanYNY=0;
for (int i=0; i<n; i++)
{
ny[i] = a[0]*x[i]*x[i]+a[1]*x[i]+a[2];
meanNY += ny[i];
meanYNY += y[i]*ny[i];
}
meanYNY /= n;
double meanY = sumY / n;
meanNY /= n;
double sumNY2=0, sumY2=0, sumYNY2=0;
for (int i=0; i<n; i++)
{
sumNY2 += (ny[i]-meanNY)*(ny[i]-meanNY);
sumY2 += (y[i] - meanY)*(y[i] - meanY);
sumYNY2 += (ny[i]-meanNY)*(y[i] - meanY);
}
*r=0;
if(fabs(sumNY2) > EPSILONG && fabs(sumY2) > EPSILONG)
{
*r = (sumYNY2*sumYNY2)/(sumNY2*sumY2);
}
delete []ny;

return solved;
}