//=====================================================================================
//函數說明
//函數名稱：PolyFit
//函數功能：最小二乘法曲線擬合
//使用方法：double *x ---- 存放n個數據點的X坐標
//          double *y ---- 存放n個數據點的Y坐標
//			int n -------- 給定數據點個數
//			double *a ---- 返回m-1次擬合多項式的m個系數
//          int m -------- 擬合多項式的項數，即擬合多項式的最高次為m-1。要求m<=n,且
//						   m<=20。若m>n或m>20，則本函數自動按m=min{n,20}處理
//			double *dt --- dt[0]返回擬合多項式與各數據點誤差的平方和；dt[1]返回擬合多
//						   項式與各數據點的誤差絕對值之和；dt[2]返回擬合多項式與各數據
//						   點誤差絕對值的最大值
//注意事項：擬合多項式的形式為 y = b0 + b1*(x-Xavr)...
//=====================================================================================
#include "math.h"
void PolyFit(double *x, double *y, int n, double *a, int m, double *dt)
{
	int i, j, k;
	double z, p, c, g, q, d1, d2, s[20], t[20], b[20];
	for (i = 0; i <= m-1; i++)
	{
		a[i] = 0.0;
	}
	if (m > n)
	{
		m = n;
	}
	if (m > 20)
	{
		m = 20;
	}
	z = 0.0;
	for (i = 0; i <= n-1; i++)
	{
		z = z+x[i]/(1.0 *n);
	}
	b[0] = 1.0;
	d1 = 1.0 * n;
	p = 0.0;
	c = 0.0;
	for (i = 0; i <= n-1; i++)
	{
		p = p+(x[i]-z);
		c = c+y[i];
	}
	c = c/d1;
	p = p/d1;
	a[0] = c * b[0];
	if (m > 1)
	{
		t[1] = 1.0;
		t[0] = -p;
		d2 = 0.0;
		c = 0.0;
		g = 0.0;
		for (i = 0; i <= n-1; i++)
		{
			q = x[i]-z-p;
			d2 = d2+q * q;
			c = c+y[i] *q;
			g = g+(x[i]-z) *q * q;
		}
		c = c/d2;
		p = g/d2;
		q = d2/d1;
		d1 = d2;
		a[1] = c * t[1];
		a[0] = c * t[0]+a[0];
	}
	for (j = 2; j <= m-1; j++)
	{
		s[j] = t[j-1];
		s[j-1] = -p * t[j-1]+t[j-2];
		if (j >= 3)
			for (k = j-2; k >= 1; k--)
		    {
				s[k] = -p * t[k]+t[k-1]-q * b[k];
		    }
		s[0] = -p * t[0]-q * b[0];
		d2 = 0.0;
		c = 0.0;
		g = 0.0;
		for (i = 0; i <= n-1; i++)
		{
			q = s[j];
			for (k = j-1; k >= 0; k--)
			{
				q = q *(x[i]-z)+s[k];
			}
			d2 = d2+q * q;
			c = c+y[i] *q;
			g = g+(x[i]-z) *q * q;
		}
		c = c/d2;
		p = g/d2;
		q = d2/d1;
		d1 = d2;
		a[j] = c * s[j];
		t[j] = s[j];
		for (k = j-1; k >= 0; k--)
		{
			a[k] = c * s[k]+a[k];
			b[k] = t[k];
			t[k] = s[k];
		}
	}
	dt[0] = 0.0;
	dt[1] = 0.0;
	dt[2] = 0.0;
	for (i = 0; i <= n-1; i++)
	{
		q = a[m-1];
		for (k = m-2; k >= 0; k--)
		{
			q = a[k]+q *(x[i]-z);
		}
		p = q-y[i];
		if (fabs(p) > dt[2])
		{
			dt[2] = fabs(p);
		}
		dt[0] = dt[0]+p * p;
		dt[1] = dt[1]+fabs(p);
	}
	return ;
}