简单思路如下：

1，采用目标函数对多项式系数求偏导，得到最优值条件，组成一个方程组；

2，方程组的解法采用行列式变换（两次变换：普通行列式——三角行列式——对角行列式——求解），行列式的求解算法上优化过一次了，目前还没有更好的思路再优化运算方法，限幅和精度准备再修改修改

目前存在的问题：

1，代码还是太粗糙

2，数学原理可行，但是计算机运算有幅度溢出和精度问题，这方面欠考虑，导致高阶大数据可能拟合失败

下面附简单注释版代码：

#include "stdafx.h"

#include "stdlib.h"

#include "math.h"

//把二维的数组与一维数组的转换，也可以直接用二维数组，只是我的习惯是不用二维数组

#define ParaBuffer(Buffer,Row,Col) (*(Buffer + (Row) * (SizeSrc + 1) + (Col)))

//

/***********************************************************************************

从txt文件里读取double型的X，Y数据

txt文件里的存储格式为

X1Y1

X2Y2

X3Y3

X4Y4

X5Y5

X6Y6

X7Y7

X8Y8

函数返回X，Y，以及数据的数目（以组为单位）

***********************************************************************************/

static int GetXY(const char* FileName, double* X, double* Y, int* Amount)

{

FILE* File = fopen(FileName, "r");

if (!File)

return -1;

for (*Amount = 0; !feof(File); X++, Y++, (*Amount)++)

if (2 != fscanf(File, (const char*)"%lf %lf", X, Y))

break;

fclose(File);

return 0;

}

/***********************************************************************************

打印系数矩阵，只用于调试，不具备运算功能

对于一个N阶拟合，它的系数矩阵大小是（N + 1）行（N + 2）列

double* Para：系数矩阵存储地址

int SizeSrc：系数矩阵大小（SizeSrc）行（SizeSrc + 1）列

***********************************************************************************/

static int PrintPara(double* Para, int SizeSrc)

{

int i, j;

for (i = 0; i < SizeSrc; i++)

{

for (j = 0; j <= SizeSrc; j++)

printf("%10.6lf ", ParaBuffer(Para, i, j));

printf("\r\n");

}

printf("\r\n");

return 0;

}

/***********************************************************************************

系数矩阵的限幅处理，防止它溢出，目前这个函数很不完善，并不能很好地解决这个问题

原理：矩阵解行列式，同一行乘以一个系数，行列式的解不变

当然，相对溢出问题，还有一个精度问题，也是同样的思路，现在对于这两块的处理很不完善，有待优化

以行为单位处理

***********************************************************************************/

static int ParalimitRow(double* Para, int SizeSrc, int Row)

{

int i;

doubleMax, Min, Temp;

for (Max = abs(ParaBuffer(Para, Row, 0)), Min = Max, i = SizeSrc; i; i--)

{

Temp = abs(ParaBuffer(Para, Row, i));

if (Max < Temp)

Max = Temp;

if (Min > Temp)

Min = Temp;

}

Max = (Max + Min) * 0.000005;

for (i = SizeSrc; i >= 0; i--)

ParaBuffer(Para, Row, i) /= Max;

return 0;

}

/***********************************************************************************

同上，以矩阵为单位处理

***********************************************************************************/

static int Paralimit(double* Para, int SizeSrc)

{

int i;

for (i = 0; i < SizeSrc; i++)

if (ParalimitRow(Para, SizeSrc, i))

return -1;

return 0;

}

/***********************************************************************************

系数矩阵行列式变换

***********************************************************************************/

static int ParaPreDealA(double* Para, int SizeSrc, int Size)

{

int i, j;

for (Size -= 1, i = 0; i < Size; i++)

{

for (j = 0; j < Size; j++)

ParaBuffer(Para, i, j) = ParaBuffer(Para, i, j) * ParaBuffer(Para, Size, Size) - ParaBuffer(Para, Size, j) * ParaBuffer(Para, i, Size);

ParaBuffer(Para, i, SizeSrc) = ParaBuffer(Para, i, SizeSrc) * ParaBuffer(Para, Size, Size) - ParaBuffer(Para, Size, SizeSrc) * ParaBuffer(Para, i, Size);

ParaBuffer(Para, i, Size) = 0;

ParalimitRow(Para, SizeSrc, i);

}

return 0;

}

/***********************************************************************************

系数矩阵行列式变换，与ParaPreDealA配合

完成第一次变换，变换成三角矩阵

***********************************************************************************/

static int ParaDealA(double* Para, int SizeSrc)

{

int i;

for (i = SizeSrc; i; i--)

if (ParaPreDealA(Para, SizeSrc, i))

return -1;

return 0;

}

/***********************************************************************************

系数矩阵行列式变换

***********************************************************************************/

static int ParaPreDealB(double* Para, int SizeSrc, int OffSet)

{

int i, j;

for (i = OffSet + 1; i < SizeSrc; i++)

{

for (j = OffSet + 1; j <= i; j++)

ParaBuffer(Para, i, j) *= ParaBuffer(Para, OffSet, OffSet);

ParaBuffer(Para, i, SizeSrc) = ParaBuffer(Para, i, SizeSrc) * ParaBuffer(Para, OffSet, OffSet) - ParaBuffer(Para, i, OffSet) * ParaBuffer(Para, OffSet, SizeSrc);

ParaBuffer(Para, i, OffSet) = 0;

ParalimitRow(Para, SizeSrc, i);

}

return 0;

}

/***********************************************************************************

系数矩阵行列式变换，与ParaPreDealB配合

完成第一次变换，变换成对角矩阵，变换完毕

***********************************************************************************/

static int ParaDealB(double* Para, int SizeSrc)

{

int i;

for (i = 0; i < SizeSrc; i++)

if (ParaPreDealB(Para, SizeSrc, i))

return -1;

for (i = 0; i < SizeSrc; i++)

{

if (ParaBuffer(Para, i, i))

{

ParaBuffer(Para, i, SizeSrc) /= ParaBuffer(Para, i, i);

ParaBuffer(Para, i, i) = 1.0;

}

}

return 0;

}

/***********************************************************************************

系数矩阵变换

***********************************************************************************/

static int ParaDeal(double* Para, int SizeSrc)

{

PrintPara(Para, SizeSrc);

Paralimit(Para, SizeSrc);

PrintPara(Para, SizeSrc);

if (ParaDealA(Para, SizeSrc))

return -1;

PrintPara(Para, SizeSrc);

if (ParaDealB(Para, SizeSrc))

return -1;

return 0;

}

/***********************************************************************************

最小二乘法的第一步就是从XY数据里面获取系数矩阵

double* Para：系数矩阵地址

const double* X：X数据地址

const double* Y：Y数据地址

int Amount：XY数据组数

int SizeSrc：系数矩阵大小（SizeSrc）行（SizeSrc + 1）列

***********************************************************************************/

static int GetParaBuffer(double* Para, const double* X, const double* Y, int Amount, int SizeSrc)

{

int i, j;

for (i = 0; i < SizeSrc; i++)

for (ParaBuffer(Para, 0, i) = 0, j = 0; j < Amount; j++)

ParaBuffer(Para, 0, i) += pow(*(X + j), 2 * (SizeSrc - 1) - i);

for (i = 1; i < SizeSrc; i++)

for (ParaBuffer(Para, i, SizeSrc - 1) = 0, j = 0; j < Amount; j++)

ParaBuffer(Para, i, SizeSrc - 1) += pow(*(X + j), SizeSrc - 1 - i);

for (i = 0; i < SizeSrc; i++)

for (ParaBuffer(Para, i, SizeSrc) = 0, j = 0; j < Amount; j++)

ParaBuffer(Para, i, SizeSrc) += (*(Y + j)) * pow(*(X + j), SizeSrc - 1 - i);

for (i = 1; i < SizeSrc; i++)

for (j = 0; j < SizeSrc - 1; j++)

ParaBuffer(Para, i, j) = ParaBuffer(Para, i - 1, j + 1);

return 0;

}

/***********************************************************************************

整个计算过程

***********************************************************************************/

int Cal(const double* BufferX, const double* BufferY, int Amount, int SizeSrc, double* ParaResK)

{

double* ParaK = (double*)malloc(SizeSrc * (SizeSrc + 1) * sizeof(double));

GetParaBuffer(ParaK, BufferX, BufferY, Amount, SizeSrc);

ParaDeal(ParaK, SizeSrc);

for (Amount = 0; Amount < SizeSrc; Amount++, ParaResK++)

*ParaResK = ParaBuffer(ParaK, Amount, SizeSrc);

free(ParaK);

return 0;

}

/***********************************************************************************

测试main函数

数据组数：20

阶数：5

***********************************************************************************/

int main(int argc, char* argv[])

{

//数据组数

int Amount;

//XY缓存，系数矩阵缓存

double BufferX[1024], BufferY[1024], ParaK[6];

if (GetXY((const char*)"test.txt", (double*)BufferX, (double*)BufferY, &Amount))

return 0;

//运算

Cal((const double*)BufferX, (const double*)BufferY, Amount, sizeof(ParaK) / sizeof(double), (double*)ParaK);

//输出系数

for (Amount = 0; Amount < sizeof(ParaK) / sizeof(double); Amount++)

printf("ParaK[%d] = %lf!\r\n", Amount, ParaK[Amount]);

//屏幕暂留

system("pause");

return 0;

}