Интерполяционная Ф-ла Ньютона

#include <stdlib.h>
#include <conio.h>
#include <stdio.h>
#include <iostream>
#include <locale.h>
#include <math.h>
using namespace std;
int factorial (int f);
 
int   main ()
{
        setlocale(LC_CTYPE,"Russian");
       
                        float X[]={ 300,400,500,600 };
                        float Y[]={52.88,65.61,78.07,99.24};
                       
                        int n=3;
 
       
                        float testX;
                        cout<<"Вводим проверяемый Х\n";
                        cin>>testX;
                        if(testX>X[3]||testX<X[0])
                        {cout<<"Недопустимое значение\n";return 0;}
 
 
                float ** razn= new float *[n+n+1]; /*создание динамического массива для таблицы конечных разностей
                                                                                  [n+n+1]столько у нас строк в таблице*/
        int k=0;
        for (int i=0;i<=n;++i)
        {
                razn[i] = new float [n+1+k];
                --k;// выделяем память для каждой строки определенное кол столбцов [n+n+1+k]
                //к- уменьшает число столбцов в строчке на 1, т.к.
                //таблица разностей имеет диагональный вид
        }
        for ( int i=0;i<=n;++i)// заполняем нулевой столбец массива значениями функции в узлах
        {
                razn [i][0]=Y[i];
                cout<<razn[i][0]<<" "<<endl;
        }
        cout<<endl<<"----------------"<<endl;
        int r=n;
 
        for (int j=1;j<=n+n;++j)
        {
                for(int i=0;i<r;++i)
                {
                        int s=-i;
                        razn[i][j]=razn[-(s-1)][j-1]-razn[i][j-1];//получаем таблицу разностей
                        cout<<"r[ "<<i<<"][ "<<j<<"]"<<razn[i][j]<<" "<<endl;
                }
                cout<<endl<<"----------------"<<endl;
                r--;//это чтобы уменьшать число столбцов, в строке, приводить к диагональному виду таблицу разности
               if(r==0)break;
        }
 
float sum,sum1,sum3,sum2;
float st=X[1]-X[0];
 
float *q = new float;
float tmp=1;
 
sum=Y[n];int p=1,u=n;
for(int i=0;i<=n;i++)
 
        {
                        q[i]=(testX-X[u]);
                        tmp*=q[i];
                        sum1=razn[u-1][i+1]*tmp;
                        sum2=pow(st,p)*factorial(p);
                        sum3=sum1/sum2;
                        sum=sum+sum3;
                        cout<<"sum= "<<sum<<"\n";
                        p++;u--;
                        if(u-1<0)break;//чтобы не выйти за границы таблицы конечн.разностей
        }      
 
        cout<<"Результата вычислений Р("<<testX<<")= "<<sum<<endl;
        return 0;
}
       
int factorial (int f)
{
int product = 1;
while (f> 0) {product = f *product;f--;}
return product;
}