Интерполяционная формула Бесселя

// Программа для реализации формулы Бесселя
#include<iostream>
#include <cstdlib>
#include <iomanip> 
using namespace std;
 
 
// вычисление u, упомянутый в формуле
float ucal(float u, int n)
{
    if (n == 0)
        return 1;
 
    float temp = u;
    for (int i = 1; i <= n / 2; i++)
        temp = temp * (u - i);
 
    for (int i = 1; i < n / 2; i++)
        temp = temp * (u + i);
 
    return temp;
}
 
// расчет факториала заданного числа n
int fact(int n)
{
    int f = 1;
    for (int i = 2; i <= n; i++)
        f *= i;
 
    return f;
}
 
int main()
{
    // количество заданных значений
    int const n = 6;
    float x[n];
    const int PI = 3.14;
    int a, b;
    cout << "vvedite a:" << endl;
    cin >> a;
    cout << "vvedite b:" << endl;
    cin >> b;
    cout << "\n" << endl;
    const int m = 6;
    double yy[m];
    cout << "\n znacheniya yy=" << endl;
    for (int i = 0; i < m; i++) {
        //int x;
        x[i] = (b - a)*i*PI;
        yy[i] = x[i] * sin(x[i]);
        cout << yy[i] << endl;
    }
    float y[m][1];
    for (int i = 0; i < m; i++) {
        y[i][1] = 0;
        y[i][0] = yy[i];
    }
 
 
    // Расчет центральной разностной таблицы 
    for (int i = 1; i < n; i++)
        for (int j = 0; j < n - i; j++)
            y[j][i] = y[j + 1][i - 1] - y[j][i - 1];
 
    // Отображение центральной разностной таблицы
    cout << "\n" << endl;
    cout << "\n Centralnaya raznostnaya tablica:" << endl;
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n - i; j++)
            cout << setw(4) << y[i][j] << "\t";
        cout << endl;
    }
 
    // значение для интерполяции при
    float value = x[3];
 
    // Инициализация и сумма
    float sum = (y[2][0] + y[3][0]) / 2;
 
    // k является источником, то есть f(0)
    int k;
    if (n % 2) // происхождение для нечетного 
        k = n / 2;
    else
        k = n / 2 - 1;
 
    float u = (value - x[k]) / (x[1] - x[0]);
 
    // Решение по формуле Бесселя
    for (int i = 1; i < n; i++) {
        if (i % 2)
            sum = sum + ((u - 0.5) *
                ucal(u, i - 1) * y[k][i]) / fact(i);
        else
            sum = sum + (ucal(u, i) *
                (y[k][i] + y[--k][i]) / (fact(i) * 2));
    }
 
    cout << "Value at " << value << " is " << sum << endl;
         cout << value*sin(value) << endl;
    system("pause");
    return 0;
}