cal.cpp

vc下实现

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
double **TwoArrayAlloc(int, int);
void TwoArrayFree(double **);
double SMOOTH5(double *x,double *y,double **z,int n,int m,double **a,int p,int q,double *dt1,double *dt2,double *dt3)
{
	int i,j,k,l,kk;
	double apx[20],apy[20],bx[20],by[20],u[20][20],**v,xx,yy;
	double t[20],t1[20],t2[20],d1,d2,g,g1,g2,x1,x2,y1,dd,dt;

	v=TwoArrayAlloc(20,m);
	xx=0;
	for(i=1;i<=n;i++)
		xx=xx+x[i-1]/n;
	yy=0;
	for(i=1;i<=m;i++)
		yy=yy+y[i-1]/m;
	d1=n;
	apx[0]=0;
	for(i=1;i<=n;i++)
		apx[0]=apx[0]+x[i-1]-xx;
	apx[0]=apx[0]/d1;
	for(j=1;j<=m;j++)
	{
		v[0][j-1]=0;
		for(i=1;i<=n;i++)
			v[0][j-1]=v[0][j-1]+z[i-1][j-1];
		v[0][j-1]=v[0][j-1]/d1;
	}
	if(p>1)
	{
		d2=0;
		apx[1]=0;
		for(i=1;i<=n;i++)
		{
			g=x[i-1]-xx-apx[0];
			d2=d2+g*g;
			apx[1]=apx[1]+(x[i-1]-xx)*g*g;
		}
		apx[1]=apx[1]/d2;
		bx[1]=d2/d1;
		for(j=1;j<=m;j++)
		{
			v[1][j-1]=0;
			for(i=1;i<=n;i++)
			{
				g=x[i-1]-xx-apx[0];
				v[1][j-1]=v[1][j-1]+z[i-1][j-1]*g;
			}
			v[1][j-1]=v[1][j-1]/d2;
		}
		d1=d2;
	}
	for(k=3;k<=p;k++)
	{
		d2=0;
		apx[k-1]=0;
		for(j=1;j<=m;j++)
			v[k-1][j-1]=0;
		for(i=1;i<=n;i++)
		{
			g1=1;
			g2=x[i-1]-xx-apx[0];
			for(j=3;j<=k;j++)
			{
				g=(x[i-1]-xx-apx[j-2])*g2-bx[j-2]*g1;
				g1=g2;
				g2=g;
			}
			d2=d2+g*g;
			apx[k-1]=apx[k-1]+(x[i-1]-xx)*g*g;
			for(j=1;j<=m;j++)
				v[k-1][j-1]=v[k-1][j-1]+z[i-1][j-1]*g;
		}
		for(j=1;j<=m;j++)
			v[k-1][j-1]=v[k-1][j-1]+z[i-1][j-1]*g;
		apx[k-1]=apx[k-1]/d2;
		bx[k-1]=d2/d1;
		d1=d2;
	}
	d1=m;
	apy[0]=0;
	for(i=1;i<=m;i++)
		apy[0]=apy[0]+y[i-1]-yy;
	apy[0]=apy[0]/d1;
	for(j=1;j<=p;j++)
	{
		u[j-1][0]=0;
		for(i=1;i<=m;i++)
			u[j-1][0]=u[j-1][0]+v[j-1][i-1];
		u[j-1][0]=u[j-1][0]/d1;
	}
	if(q>1)
	{
		d2=0;
		apy[1]=0;
		for(i=1;i<=m;i++)
		{
			g=y[i-1]-yy-apy[0];
			d2=d2+g*g;
			apy[1]=apy[1]+(y[i-1]-yy)*g*g;
		}
		apy[1]=apy[1]/d2;
		by[1]=d2/d1;
		for(j=1;j<=p;j++)
		{
			u[j-1][1]=0;
			for(i=1;i<=m;i++)
			{
				g=y[i-1]-yy-apy[0];
				u[j-1][1]=u[j-1][1]+v[j-1][i-1]*g;
			}
			u[j-1][1]=u[j-1][1]/d2;
		}
		d1=d2;
	}
	for(k=3;k<=q;k++)
	{
		d2=0;
		apy[k-1]=0;
		for(j=1;j<=p;j++)
			u[j-1][k-1]=0;
		for(i=1;i<=m;i++)
		{
			g1=1;
			g2=y[i-1]-yy-apy[0];
			for(j=3;j<=k;j++)
			{
				g=(y[i-1]-yy-apy[j-2])*g2-by[j-2]*g1;
				g1=g2;
				g2=g;
			}
			d2=d2+g*g;
			apy[k-1]=apy[k-1]+(y[i-1]-yy)*g*g;
			for(j=1;j<=p;j++)
				u[j-1][k-1]=u[j-1][k-1]+v[j-1][i-1]*g;
		}
		for(j=1;j<=p;j++)
			u[j-1][k-1]=u[j-1][k-1]/d2;
		apy[k-1]=apy[k-1]/d2;
		by[k-1]=d2/d1;
		d1=d2;
	}
	v[0][0]=1;
	v[1][0]=-apy[0];
	v[1][1]=1;
	for(i=1;i<=p;i++)
	for(j=1;j<=q;j++)
		a[i-1][j-1]=0;
	for(i=3;i<=q;i++)
	{
		v[i-1][i-1]=v[i-2][i-2];
		v[i-1][i-2]=-apy[i-2]*v[i-2][i-2]+v[i-2][i-3];
		if(i>=4)
			for(k=i-2;k>=2;k--)
				v[i-1][k-1]=-apy[i-2]*v[i-2][k-1]+v[i-2][k-2]-by[i-2]*v[i-3][k-1];
		v[i-1][0]=-apy[i-2]*v[i-2][0]-by[i-2]*v[i-3][0];
	}
	for(i=1;i<=p;i++)
	{
		if(i==1)
		{
			t[0]=1;
			t1[0]=1;
		}
		else if(i==2)
		{
			t[0]=-apx[0];
			t[1]=1;
			t2[0]=t[0];
			t2[1]=t[1];
		}
		else 
		{
			t[i-1]=t2[i-2];
			t[i-2]=-apx[i-2]*t2[i-2]+t2[i-3];
			if(i>=4)
				for(k=i-2;k>=2;k--)
					t[k-1]=-apx[i-2]*t2[k-1]+t2[k-2]-bx[i-2]*t1[k-1];
			t[0]=-apx[i-2]*t2[0]-bx[i-2]*t1[0];
			t2[i-1]=t[i-1];
			for(k=i-1;k>=1;k--)
			{
				t1[k-1]=t2[k-1];
				t2[k-1]=t[k-1];
			}
		}
		for(j=1;j<=q;j++)
		for(k=i;k>=1;k--)
		for(l=j;l>=1;l--)
			a[k-1][l-1]=a[k-1][l-1]+u[i-1][j-1]*t[k-1]*v[j-1][l-1];
	}
	*dt1=0;
	*dt2=0;
	*dt3=0;
	for(i=1;i<=n;i++)
	{
		x1=x[i-1]-xx;
		for(j=1;j<=m;j++)
		{
			y1=y[j-1]-yy;
			x2=1;
			dd=0;
			for(k=1;k<=p;k++)
			{
				g=a[k-1][q-1];
				for(kk=q-1;kk>=1;kk--)
					g=g*y1+a[k-1][kk-1];
				g=g*x2;
				dd=dd+g;
				x2=x2*x1;
			}
			dt=dd-z[i-1][j-1];
			if(fabs(dt)>(*dt3))
				*dt3=fabs(dt);
			*dt1=(*dt1)+dt*dt;
			*dt2=(*dt2)+fabs(dt);
		}
	}
	TwoArrayFree(v);
	return(1);
}
void main()
{
	int i,j,n,m,p,q;
	double *x,*y,**z,**a,dt1,dt2,dt3;

	n=11;
	m=21;
	p=6;
	q=5;

	x=(double *)calloc(n,sizeof(double));
	if(x==NULL)
		exit(0);
	y=(double *)calloc(m,sizeof(double));
	if(y==NULL)
		exit(0);
	z=TwoArrayAlloc(n,m);
	a=TwoArrayAlloc(p,q);

/*	for(i=1;i<=n;i++)
		x[i-1]=0.2*(i-1);
	for(i=1;i<=m;i++)
		y[i-1]=0.1*(i-1);
	for(i=1;i<=n;i++)
	for(j=1;j<=m;j++)
		z[i-1][j-1]=exp(x[i-1]*x[i-1]-y[j-1]*y[j-1]);
*/
	SMOOTH5(x,y,z,n,m,a,p,q,&dt1,&dt2,&dt3);
	for(i=1;i<=p;i++)
	{
		for(j=1;j<=q;j++)
			printf("%.10f\t",a[i-1][j-1]);
		printf("\n");
	}
	printf("\n%.10e\t%.10e\t%.10e\n",dt1,dt2,dt3);
	
	free(x);
	free(y);
	TwoArrayFree(z);
	TwoArrayFree(a);
}