Интерполяция функции по Ньютону

#include <cstdlib>
#include <iostream>
#include <math.h>
#include <conio.h>
#include <locale>
 
using namespace std;
///////////////////////////////////////////////////////////////////////////////
// Некоторые глобальные переменные
///////////////////////////////////////////////////////////////////////////////
const int maxN = 50;
int N;
double x[maxN]; // Точки на [A,B]
double yd[maxN],yt[maxN],Temp[maxN],fuctor[maxN]; // Коэфициенты при многочленах для Интерполяции
double mnog[maxN]; // Многочлен по теореме Виетта
double inter[maxN]; // Интерполяционный многочлен Ньютона
///////////////////////////////////////////////////////////////////////////////
// Структура для теоремы Виетта
struct kf
{
    int por;
    double zna4;
};
kf koef1[maxN*maxN];
kf koef2[maxN*maxN];
// 1/n!
double fuct_1(int q)
{
    double temp=1;
    for(int j=2;j<=q;j++)temp/=j;
    return (temp);
}
///////////////////////////////////////////////////////////////////////////////
// Свободная функция
double fFree(double x)
{
    return ( exp( 0.5*x ) );
}
///////////////////////////////////////////////////////////////////////////////
// Построение многочленов по теореме Виетта
///////////////////////////////////////////////////////////////////////////////
// Находим конечные разности
void dy()
{
    double yt[maxN], yk[maxN];
    for(int i = 0; i<=N; i++)yt[i] = Temp[i];
    for(int k=0;k<=N;k++)
    {
        yd[k] = yt[0];
        for(int i=0;i<N-k;i++)
        {
            yk[i] = yt[i+1] - yt[i];
            yt[i] = yk[i];
        }
    }
}
// Находим коэфициенты при многочленах
void Ci()
{
    dy();
    for(int i=0;i<=N;i++)yd[i]*=pow(h,-i)*fuctor[i];
}
// Находит сумму попарных произведений
double sum(int n,int q)
{
    double s=0;
    for(int ii=0;ii<=n;ii++)s+=koef1[ii].zna4;
    if( q%2 == 0 ) s = -s;
    return ( s );
}
// Коэфициенты многочлена по теореме Виетта
void TVietta(int w)
{
    for(int i=0; i <= N; i++)mnog[i] = 0;
    for(int i=0; i <= N-w; i++)
    {
        koef1[i].zna4 = x[i];
        koef1[i].por = i;
    }
    mnog[N-w+1] = 1;    
    for(int i=0; i<=N-w; i++)mnog[N-w]+=x[i];mnog[N-w]=-mnog[N-w];
    int l = N-w;
    int soka=1;
 
while( soka <= N-w )
{
    int m = l;
    l = -1;
    for(int i = 0; i <= m; i++)
    {
        for(int j = koef1[i].por+1; j <= N-w; j++)
        {
            l++;
            koef2[l].zna4 = koef1[i].zna4 * x[j];
            koef2[l].por = j;
        }
    }
    for(int i = 0; i <= l; i++)
    {
        koef1[i].zna4 = koef2[i].zna4;
        koef1[i].por = koef2[i].por;
    }
    mnog[N-w-soka] = sum(l,soka);
    soka++;
}
}
// Интерполяция по Ньютону, собственно строим сам многочлен
void Interpol_P(int r)
{
    Ci();
    for(int e = 0; e <= N; e++)inter[e] = 0;
    for(int e = 0; e <= N; e++)
    {
        TVietta( e+1 );
        for(int i = 0; i <= N-e; i++)inter[i] += yd[N-e]*mnog[i];
    }
    if(r != 0)for(int e = 0; e <= N; e++)printf("X^%1.1i:  %4.25f\n",e,inter[e]);
}
///////////////////////////////////////////////////////////////////////////////
// Интерполяция свободной функции
///////////////////////////////////////////////////////////////////////////////
void Interpol_fFree()
{
   for(int i=0;i<=N;i++)
   {
           Temp[i]=fFree(x[i]);
           fuctor[i]=fuct_1(i);
   }
   Interpol_P(1);
}
int main()
{
    N = 10;
    h = (2.0-1.0)/N;
    for(int i = 0; i <= N; i++)x[i] = 1 + i * h;
    Interpol_fFree();
    return (0);
}