написанные программы / интерполяция функции полиномом лагранжа по схеме горнера / var0 / Interpolation / Interpolation / ILib

//-----------------------------
//Interpolation Functions
//-----------------------------
//19.04.2004 $ Serkov Alexander
//-----------------------------
#include "ILib.h"
//-----------------------------
POINTS *AllocPoints(UINT32 n)
{
	POINTS *points;

	points=malloc(sizeof(POINTS));
	if(points==NULL)
		return NULL;
	points->Points=malloc(sizeof(POINT)*n);
	if(points->Points==NULL)
	{
		free(points);
		return NULL;
	}
	points->Count=n;
	points->L=NULL;

	return points;
}
//-----------------------------
POINTS *BuildPoints(MATHFUNC f,DOUBLE from,DOUBLE to,DOUBLE step)
{
	POINTS *points;
	UINT32 count,i;

	if((f==NULL)||(to<=from)||(step<=0.0))
		return NULL;
	//TODO: Не очень точная формула.
	//		Может выходить за границы
	//		диапазона, хотя не очень это
	//		и важно. (Наверное...)
	count=(UINT32)((to-from)/step+1);
	
	points=AllocPoints(count);
	if(points==NULL)
		return NULL;

	for(i=0;i<count;i++)
	{
		points->Points[i].x=from;
		points->Points[i].y=f(from);
		from+=step;
	}

	return points;
}
//-----------------------------
void FreePoints(POINTS *points)
{
	if(points!=NULL)
	{
		if(points->Points!=NULL)
			free(points->Points);
		if(points->L!=NULL)
			free(points->L);
		free(points);
	}
}
//-----------------------------
void *BuildL(POINTS *points)
{
	UINT32 i,j;

	if((points==NULL)||(points->Points==NULL)||(points->Count<2))
		return NULL;
	points->L=malloc(sizeof(DOUBLE)*(points->Count));
	if(points->L==NULL)
		return NULL;
	for(i=0;i<points->Count;i++)
	{
		points->L[i]=1.0;
		for(j=0;j<points->Count;j++)
		{
			if(i==j) continue;
			points->L[i]*=(points->Points[i].x
						-points->Points[j].x);
		}
		points->L[i]=points->Points[i].y/points->L[i];
	}

	return points->L;
}
//-----------------------------
DOUBLE __AllSumsX(POINTS *points,
				UINT32 from,
				UINT32 count,
				UINT32 except)
{
	UINT32 i;
	DOUBLE ret;

	if(from==except)
		from++;
	if(count==1)
		return points->Points[from].x;
	if(count<=0)
		return 0.0;

	ret=0.0;
	//TODO: не до points->Count,
	//		а до чего то меньшего.
	//		Подумать. (В целях оптимизации,
	//		потому как в __AllSumsX
	//		и так всё учитывается.)
	for(i=from+1;i<points->Count;i++)
		if(i!=except)
			ret+=__AllSumsX(points,i,count-1,except);
	return (ret*points->Points[from].x);
}
//-----------------------------
DOUBLE AllSumsX(POINTS *points,UINT32 count,UINT32 except)
{
	DOUBLE res=0.0;
	UINT32 i;
	
	for(i=0;i<points->Count;i++)
		if(i!=except)
			res+=__AllSumsX(points,i,count,except);
	return res;
}
//-----------------------------
DOUBLE CalcLx(POINTS *points,UINT32 Li,UINT32 st)
{
	if((points==NULL)||
		(points->Points==NULL)||
		(st>=points->Count)||
		(Li>=points->Count))
			return 0.0;
	if(st==(points->Count-1))
		return 1.0;

	return AllSumsX(points,points->Count-st-1,Li);
}
//-----------------------------
POLINOM *BuildPolinom(POINTS *points)
{
	POLINOM *p;
	UINT32 i,j;
	DOUBLE sign;

	if(BuildL(points)==NULL)
		return NULL;
	p=malloc(sizeof(POLINOM));
	if(p==NULL)
		return NULL;
	p->p=malloc(sizeof(DOUBLE)*points->Count);
	if(p==NULL)
	{
		free(p);
		return NULL;
	}
	p->Count=points->Count;
	sign=((points->Count%2)==1)?1.0:-1.0;
	for(i=0;i<points->Count;i++)
	{
		p->p[i]=0.0;
		for(j=0;j<points->Count;j++)
			p->p[i]+=sign*points->L[j]*CalcLx(points,j,i);
		sign=-sign;
	}
	return p;
}
//-----------------------------
void FreePolinom(POLINOM *p)
{
	if(p!=NULL)
	{
		if(p->p!=NULL)
			free(p->p);
		free(p);
	}
}
//-----------------------------
DOUBLE CalcPolinom(POLINOM *p,DOUBLE x)
{
	DOUBLE res=0.0;
	UINT32 i;

	if((p==NULL)||(p->p==NULL))
		return 0.0;

	//По Горнеру...
	for(i=0;i<p->Count;i++)
	{
		res*=x;
		res+=p->p[p->Count-i-1];
	}
	return res;
}
//-----------------------------