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

//-----------------------------
//Interpolation Functions
//-----------------------------
//19.04.2004 $ Serkov Alexander
//-----------------------------
#include "ILib.h"
#include <stdio.h>
#include <math.h>
//-----------------------------
#define FROM  0.0
#define TO    4.0
#define STEP  0.3
//-----------------------------
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;

	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;
}
//-----------------------------
                      //-----------------------------
DOUBLE Funct(DOUBLE x)
{
	return (exp(x)*exp(x)*x);
}
//-----------------------------
int main()
{
	POINTS *points;
	POLINOM *p;
	DOUBLE x;

	points=BuildPoints(Funct,FROM,TO,STEP);
	if(points!=NULL)
		printf("POINTS: %lu\n",points->Count);
	p=BuildPolinom(points);
	FreePoints(points);

	if(p==NULL)
	{
		printf("ERROR: Not Enought Memory.");
		return -1;
	}

	for(x=FROM;x<=(TO+STEP*2);x+=(STEP/2))
	{
		printf("X: %e  Yreal: %e  Ypolin: %e\n",
			x,Funct(x),
			CalcPolinom(p,x));
	}
	FreePolinom(p);

	return 0;
}
//----