Численное интегрирование по Симпсону

//Вызов
auto[sum_s, count_iterations_s] = finding_integral_simpson(f, 1., 1.4, e);
    wcout << L"Интегрирование по формуле Симпсона" << endl;
    wcout << L"Значение интеграла " << fixed << setprecision(12u) << sum_s << endl;
    wcout << L"Достигнуто за " << fixed << setprecision(12u) << count_iterations_s << L" итераций" << endl;
    wcout << endl;
 
//Интегрирование Симпсона
tuple<double, size_t> finding_integral_simpson(const function<double(double)> &f, const double &left, const double &right, const double &e)
{
    double x;
    double h;
    double sum;
    double max_derivative;
    size_t count_rectangle_all;
    size_t rectangle_now;
    size_t n = 1u;
 
    //Поиск M
    x = left;
    h = 0.1;
    count_rectangle_all = (right - left) / h;
    rectangle_now = 0u;
    vector<pair<double, double>> Data_;
    function_point_cloud(f, Data_, -10., 10., 1000);
    max_derivative = abs(get_derivative_n_in_point(Data_, x, 4u));
    for (; rectangle_now <= count_rectangle_all; x += h, ++rectangle_now)
    {
        double now_d = abs(get_derivative_n_in_point(Data_, x, 4u));
        max_derivative = now_d > max_derivative ? now_d : max_derivative;
    }
 
    for (;; ++n)
    {
        //Подготовка к интегрированию
        h = (right - left) / pow(2, n);
        count_rectangle_all = pow(2, n);
        rectangle_now = 0u;
        sum = 0.;
        x = left + h;
 
        //Интегрирование
        bool even = false;
        for (; rectangle_now < count_rectangle_all - 1; x += h, ++rectangle_now)
        {
            sum += (even ? 2. : 4.) * f(x);
            even = !even;
        }
        sum += f(left) + f(x + h);
 
        //Проверка точности
        //double e_fact = max_derivative * pow((right - left), 5) / (180. * pow((count_rectangle_all), 4));
        //double e_fact = pow(h, 4) * max_derivative * (right - left) / 2880.;
        double e_fact = max_derivative * pow((right - left), 5) / (2880. * pow(count_rectangle_all, 4));
        if (e_fact < e)
        {
            break;
        }
    }
    return { sum * h / 3., n };
}
 
//Взятие любой производной
double get_derivative_n_in_point(const vector<pair<double, double>> &points, const double &x, const size_t &n)
{
    vector<pair<double, double>> box_old(points);
    vector<pair<double, double>> box_new;
    for (size_t p = 0u; p < n; ++p)
    {
        box_new = vector<pair<double, double>>(box_old.size() - 2);
        for (long long u = 1u; u < box_old.size() - 1; ++u)
        {
            double y = (box_old[u + 1].second - box_old[u - 1].second) / (box_old[u + 1].first - box_old[u - 1].first);
            double x = box_old[u].first;
            box_new[u - 1] = make_pair(x, y);
        }
        box_old = box_new;
    }
    size_t ind = 0u;
    for (const auto &value : box_new)
    {
        if (x <= value.first)
        {
            break;
        }
        ++ind;
    }
    //wcout << ind << endl;
    double res = abs(box_new[ind].first - x) < abs(box_new[ind - 1].first - x) ? box_new[ind].second : box_new[ind - 1].second;
    return res;
}