Численное интегрирование - C (СИ)

#include <iostream>
#include <cmath>
#include <windows.h>
 
//Пределы интегрирования
const double a = 0;
const double b = 1;
//Точность интегрирования
const double e = 0.0005;
//Множитель аргумента подинтегральной функции
const double k = 10;
 
//Знак-множитель аргумента подинтегральной функции
enum sign
{
    plus = 1,
    minus = -1
};
 
//Подинтегральная функция
long double f(long double x, sign sn)
{
    return exp(sn*k*x);
}
 
//Вычисление приближенного значения интеграла на отрезке [x, x+step]
long double q(long double x, long double step, sign sn)
{
    return (step/2.0*(f(x, sn)+f((x+step), sn)));
}
 
//Вычисление приближенного значения интеграла на отрезке [a, b]
void CalcInt(sign sn)
{
    long double h=0.1;      //Первоначальный шаг
    long double h1=h/2.0;   //Вдвое меньший шаг
    long double Q=0;        //Искомое значение интеграла
    long double x=a;        //Текущий узел сетки
    long double Q1, Q2;     //Значения интеграла, вычисленного на текущем
                            //отрезке с шагом h и h1, соответственно
    long double LExp, RExp; //Левая и правая части условия проверки
                            //точности решения и оптимальности шага
    size_t cnt=0;           //Число операций суммирования
    while(true)
    {
        while((x+h)<b)      //Пока не достигли правого конца интервала
        {
            //Вычисляем приближенные значения интеграла на текущем отрезке
            Q1=q(x, h, sn);                 //С шагом h
            Q2=q(x, h1, sn)+q(x+h1, h1, sn);//С шагом h1
            LExp=fabs(Q1-Q2)/3.0;
            RExp=h*e/(b-a);
            //Проверяем, можно ли увеличить шаг
            if(10*LExp<=RExp)
            {
                //Увеличиваем шаг вдвое
                x+=h;
                h1=h;
                h*=2;
                Q+=Q2;
                ++cnt;
            }
            //Сохраняется ли необходимая точность
            else if(LExp<=RExp)
            {
                //Уменьшаем шаг вдвое
                x+=h;
                h=h1;
                h1/=2.0;
                Q+=Q2;
                ++cnt;
            }
            else
            {
                //Уменьшаем шаг вдвое
                h=h1;
                h1/=2.0;
            }
        }
        //Вычисляем значение интеграла на последнем отрезке
        h=b-x;
        h1=h/2.0;
        Q1=q(x, h, sn);
        Q2=q(x, h1, sn)+q(x+h1, h1, sn);
        //Проверяем, удастся ли сохранить точность на последнем отрезке
        if(fabs(Q1-Q2)/3.0<=h*e/(b-a))
        {
            //Если да, то суммируем значение последнего отрезка
            Q+=Q2;
            ++cnt;
            //...и выходим из цикла
            break;
        }
        else
        {
            //Иначе -  уменьшаем шаг вдвое
            h=h1;
            h1/=2.0;
        }
    }
    //Вывод полученных результатов
    std::cout << "Приближенное значение интеграла функции exp("
        << (sn*k) << "x) на отрезке [" << a << "; " << b << "]"
        << "\n\tравно " << Q << ". Число операций суммирования равно "
        << cnt << std::endl << std::endl;
}
 
int main()
{
    SetConsoleCP(1251);
    SetConsoleOutputCP(1251);
    CalcInt(plus);
    CalcInt(minus);
    system("pause");
    return EXIT_SUCCESS;
}