Используя интерполяционную формулу Ньютона, вычислить значение функции при заданных значениях аргумента с точностью до 0.001

// Выводим массив на консоль
// _array - массив для вывода, 
// n -размер массива
void PrintArray(double* _array, int n)
{
    for (int i = 0; i < n; i++)
    {
        printf("%lf\t", _array[i]);
    }
}
//Считываем массив с консоли
//n размер массива
double* ReadArray(int n)
{
    double* _array = new double[n];
    for (int i = 0; i < n; i++)
    {
        scanf_s("%lf", &_array[i]);
    }
    return _array;
}
int main()
{
    double *x, *y, *a, *b, *d;
    int n, numb;
    double q, p, h, f;
    printf("Input n :\r\n");
    scanf_s("%i", &n);
    printf("Input x :\r\n");
    x = ReadArray(n);
    PrintArray(x, n);
    printf("Input y :\r\n");
    y = ReadArray(n);
    printf("\r\nX = :\r\n");
    PrintArray(x, n);
    printf("\r\nY = :\r\n");
    PrintArray(y, n);
    a = new double[n], d = new double[n], b = new double[n];
    for (int i = 0; i < n - 1; i++)
        a[i] = y[i + 1] - y[i];
    for (int i = 0; i < n - 1; i++)
        b[i] = a[i + 1] - a[i];
    for (int i = 0; i < n - 1; i++)
        d[i] = b[i + 1] - b[i];
    for (int i = 0; i < n - 1; i++)
        printf("X = %lf \t Y = %lf \t A = %lf\t B = %lf \t D = %lf \r\n", x[i], y[i], a[i], b[i], d[i]);
    h = x[2] - x[1];
    printf("Input f:\r\n");
    scanf_s("%lf", &f);
    printf("Input number х:");
    scanf_s("%i", &numb);
    q = (f - x[numb]) / h;
    p = y[numb] + (q*a[numb]) + (q*(q - 1) / 2)*b[numb] + (q*(q - 1)*(q - 2) / 6)*d[numb];
    printf("P = %lf\r\n", p);
    system("Pause");
    delete[] a, d, b, x, y;
    return 0;
}

// Выводим массив на консоль
// _array - массив для вывода, 
// n -размер массива
void PrintArray(double* _array, int n)
{
    for (int i = 0; i < n; i++)
    {
        printf("%lf\t", _array[i]);
    }
}
//Считываем массив с консоли
//n размер массива
double* ReadArray(int n)
{
    double* _array = new double[n];
    for (int i = 0; i < n; i++)
    {
        scanf_s("%lf", &_array[i]);
    }
    return _array;
}
int main()
{
    /*srand(time(NULL));*/
    double *x, *y, *a, *b, *d;
    int n, numb;
    double q, p, h, f;
    printf("Input n :\r\n");
    scanf_s("%i", &n);
    printf("Input x :\r\n");
    x = ReadArray(n);
    PrintArray(x, n);
    printf("Input y :\r\n");
    y = ReadArray(n);
    printf("\r\nX = :\r\n");
    PrintArray(x, n);
    printf("\r\nY = :\r\n");
    PrintArray(y, n);
    a = new double[n], d = new double[n], b = new double[n];
    for (int i = 0; i < n - 1; i++)
        a[i] = y[i + 1] - y[i];
    for (int i = 0; i < n - 2; i++)
        b[i] = a[i + 1] - a[i];
    for (int i = 0; i < n - 3; i++)
        d[i] = b[i + 1] - b[i];
    for (int i = 0; i < n - 3; i++)
        printf("X = %lf \t Y = %lf \t A = %lf\t B = %lf \t D = %lf \r\n", x[i], y[i], a[i], b[i], d[i]);
    h = x[2] - x[1];
    printf("Input f:\r\n");
    scanf_s("%lf", &f);
    printf("Input number х:");
    scanf_s("%i", &numb);
    q = (f - x[numb]) / h;
    p = y[numb] + (q*a[numb]) + (q*(q - 1) / 2)*b[numb] + (q*(q - 1)*(q - 2) / 6)*d[numb];
    printf("P = %lf\r\n", p);
    system("Pause");
    delete[] a, d, b, x, y;
    return 0;
}