Здравствуйте, уважаемые эксперты!
Помогите, пожалуйста, разобраться, почему неправильно считает программа?
Делаю программу, которая бы считала первую производную с помощью полинома Ньютона.
На вход считывает из файла (файл во вложении): порядок полинома, интервал сетки, исходную сетку (x); значения дифференцируемой ф-ции (y); количество шагов на новой сетке, новую сетку узлов для расчета значения производной.
На выход: новая сетка (x1); значения производной функции на новой сетке.
Поскольку сетка равномерная использую формулу:

Вот программа:
[code h=300]#include <conio.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include<iostream>

// функция
double f(double x)
{
double y;
y = pow(x, 1.0 / 3.0);
return y;
}

// функция расчета конечных разностей
double *ComputationDeltaY(double *y,int N)
{
double *deltaY,*tempDeltaY;
deltaY = new double[N+1];
tempDeltaY = new double[N+1];
deltaY[0]=y[0];
for(int i=0;i<N+1;i++)
{
tempDeltaY[i]=y[i]; // массив Y теперь в массиве tempDeltaY
}
for(int i=0;i<N;i++)
{
for(int j=0;j<N-i;j++)
{
tempDeltaY[j]=tempDeltaY[j+1]-tempDeltaY[j]; // разница между Y
}
deltaY[i]=tempDeltaY[0];
}
delete[] tempDeltaY;
return deltaY;
}

// функция расчета факториала
double Factorial(int a)
{
double Fact=1;
for(int i=1;i<=a;i++)
{
Fact *= i;
}
return Fact;
}

// функция вычисления первой производной
double *PP(double *x1,double *y,int N,int M,double a,double b)
{
double *y1,h,*delta,sum,q;
y1 = new double[M+1];
h = (b - a)/N;
delta = ComputationDeltaY(y,N);
for(int i=0;i<M+1;i++)
{
sum=0;
q = (x1[i]-a)/h;
for(int j=1;j<N+1;j++)
{
double MultiSum=0;
for(int k=0;k<=j-1;k++)
{
double multi=1;
for(int l=0;l<=j-1;l++)
{
if(k!=l)
{
multi *= q - l;
}
}
MultiSum += multi/h;
}
sum += (double) (1.0/Factorial(j))*delta[j-1]*MultiSum;
}
y1[i]=sum;
}
delete[] delta;
return y1;
}

int main(){
setlocale(LC_CTYPE,"Russian");
int M, N, k, i; // вспомогательные переменные
// переменные для массивов узлов, значений функций и шага
double a, b, *x, *y, *x1, *y1;

FILE *f1;
if((f1=fopen("input.txt","rt"))==NULL)
{
printf("Файл не открывается!\n");
exit(1);
}

x = new double[N+1];
y = new double[N+1];
fscanf(f1,"%d",&N); // порядок полинома
fscanf(f1,"%lf",&a); // начало интервала
fscanf(f1,"%lf",&b); // конец интервала
for (int i=0;i<N+1;i++)
{
fscanf(f1,"%lf",&x[i]); // x
}
for(int i=0;i<N+1;i++)
{
fscanf(f1,"%lf",&y[i]); // y
}
fscanf(f1,"%d",&M); // количество шагов по новой сетке
x1 = new double[M+1];
for (int i=0;i<M+1;i++)
{
fscanf(f1,"%lf",&x1[i]); // новый x
}

y1 = PP(x1,y,N,M,a,b);
for (int i=0;i<M+1;i++)
{
printf("y'(%d)=%.4f\n",i,y1[i]);
}
fclose(f1);

delete[] y1;
delete[] y;
delete[] x1;
delete[] x;


return 0;
}

[/code]

Она считает, но неправильно...