//最小二乘法
void pcir(double X[],double Y[],double A[],int N,int M,double &DT1,double &DT2,double &DT3)
{
// X[] Y[] 使需要拟合的点，A[]是拟合曲线的系数，
double S[21],T[21],B[21];
double Z,D1,P,C,D2,G,Q,DT;
int i,j,jk;
for(i=1;i<=M;i++)// DO 5 i=1,M
{
A[i]=0.0;
}
if (M>=N) M=N;
if (M>=20) M=20;
Z=0.0;
/* for(i=1;i<=N;i++)//DO 10 i=1,N
{
Z=Z+X[i]/N;
}//10 Z=Z+X(i)/N*/
B[1]=1.0;
D1=N;
P=0.0;
C=0.0;
for(i=1;i<=N;i++)
{
P=P+(X[i]-Z);
C=C+Y[i];
}
C=C/D1;
P=P/D1;
A[1]=C*B[1];
if(M>1)
{
T[2]=1.0;
T[1]=-P;
D2=0.0;
C=0.0;
G=0.0;
for(i=1;i<=N;i++)
{
Q=X[i]-Z-P;
D2=D2+Q*Q;
C=Y[i]*Q+C;
G=(X[i]-Z)*Q*Q+G;
}
C=C/D2;
P=G/D2;
Q=D2/D1;
D1=D2;
A[2]=C*T[2];
A[1]=C*T[1]+A[1];
}
for(j=3;j<=M;j++)
{
S[j]=T[j-1];
S[j-1]=-P*T[j-1]+T[j-2];
if(j>=4)
for(jk=j-2;jk>=2;jk--)
S[jk]=-P*T[jk]+T[jk-1]-Q*B[jk];
S[1]=-P*T[1]-Q*B[1];
D2=0.0;
C=0.0;
G=0.0;
for(i=1;i<=N;i++)
{
Q=S[j];
for(jk=j-1;jk>=1;jk--)
{
Q=Q*(X[i]-Z)+S[jk];
}
D2=D2+Q*Q;
C=Y[i]*Q+C;
G=(X[i]-Z)*Q*Q+G;
}
C=C/D2;
P=G/D2;
Q=D2/D1;
D1=D2;
A[j]=C*S[j];
T[j]=S[j];
for(jk=j-1;jk>=1;jk--)
{
A[jk]=C*S[jk]+A[jk];
B[jk]=T[jk];
T[jk]=S[jk];
}
}
DT1=0.0;
DT2=0.0;
DT3=0.0;
for(i=1;i<=N;i++)
{
Q=A[M];
for(jk=M-1;jk>=1;jk--)
{
Q=Q*(X[i]-Z)+A[jk];
}
DT=Q-Y[i];
if (fabs(DT)>DT3) DT3=fabs(DT);
DT1=DT1+DT*DT;
DT2=DT2+fabs(DT);
}
//AfxMessageBox("hello");
}
