Численное интегрирование и его оптимизация

#pragma once///////////////////////////////////////////////////////////////////////////
////////////////////////////INTEGRATE//////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
 
double Simpson(double f_sh(double), double a, double b, int n = 64);//формула Симпсона (интегрирование)
double Simpson_auto(double f_sh(double), double a, double b, double eps = 1.e-9);//формула Симпсона (интегрирование с автошагом)
double bul(double f_sh(double), double a, double b, int n = 64);//формула Буля (интегрирование)
double bul_auto(double f_sh(double), double a, double b, double eps = 1.e-9);//формула Буля (интегрирование с автошагом)
 
double FAbs(double x);//модуль числа

#include "diff.h"
 
 
double FAbs(double x)//модуль
{
    if(x > 0.)
        return x;
    return -x;
}
 
double Simpson(double f_sh(double), double a, double b, int n)//h*(f0+f2n+2*(f2k)+4*(f2k-1))/3
{//интеграл от функции
    int N = 2 * n;
    double integ = 0, h = (b - a) / N, x0 = a;
    integ = f_sh(b) + f_sh(a);
    for (int i = 1; i < N; i++)
    {
        x0 += h;
        integ += 2 * f_sh(x0) * ((i + 1) % 2) + 4 * f_sh(x0) * (i % 2);
    }
    return integ * h / 3;
}
 
 
double Simpson_auto(double f_sh(double), double a, double b, double eps)//h*(f0+f2n+2*(f2k)+4*(f2k-1))/3
{
    double res1 = Simpson(f_sh, a, b, 1);
    double res2 = Simpson(f_sh, a, b, 2);
    if (FAbs(res1 - res2) > eps)
        return Simpson_auto(f_sh, a, (a + b) / 2, eps) + Simpson_auto(f_sh, (a + b) / 2, b, eps);
    else
        return res2;
}
 
 
double bul(double f_sh(double), double a, double b, int n)//формула Буля (интегрирование)
{
    int N = 4 * n;
    double integ, h = (b - a) / N, x0 = a;
    integ = 7 * (f_sh(b) + f_sh(a));
    for (int i = 1; i < N; i++)
    {
        x0 += h;
        integ += 32 * f_sh(x0) * (i % 2) + 12 * f_sh(x0) * ((i - 2) % 4 ? 0 : 1) + 14 * f_sh(x0) * ((i - 4) % 4 ? 0 : 1);
    }
    return integ * 2 * h / 45;
}
 
 
double bul_auto(double f_sh(double), double a, double b, double eps)
{
    double res1 = bul(f_sh, a, b, 1);
    double res2 = bul(f_sh, a, b, 2);
    if (FAbs(res1 - res2) > eps)
        return bul_auto(f_sh, a, (a + b) / 2, eps) + bul_auto(f_sh, (a + b) / 2, b, eps);
    else
        return res2;
}

double FAbs(double x)//модуль
{
    if(x > 0)
        return x;
    return -x;
}
 
double bul(double f_sh(double), double a, double b, int n)//формула Буля (интегрирование)
{
    int N = 4 * n;
    double integ, h = (b - a) / N, x0 = a;
    integ = 7 * (f_sh(b) + f_sh(a));
    for (int i = 1; i < N; i++)
    {
        x0 += h;
        integ += 32 * f_sh(x0) * (i % 2) + 12 * f_sh(x0) * ((i - 2) % 4 ? 0 : 1) + 14 * f_sh(x0) * ((i - 4) % 4 ? 0 : 1);
    }
    return integ * 2 * h / 45;
}
 
 
double bul_auto(double f_sh(double), double a, double b, double eps)
{
    double h2 = (b - a) / 8, h1 = h2 + h2;//для меньшего кол-ва операций деления
    double f2 = f_sh(a + h1), f4 = f_sh(a + h1 + h1), f6 = f_sh(b - h1);//для меньшего количества вызовов функции
    double res1 = 7 * (f_sh(b) + f_sh(a)) + 32 * (f2 + f6) + 12 * f4;//почти что интеграл с шагом h1
    double dr = 32 * (f_sh(a + h2) + f_sh(a + h1 + h2) + f_sh(b - h1 - h2) + f_sh(b - h2)) - 20 * (f2 + f6) + 2 * f4;// разница почти что интегралов с шагом h1 и h2
    double res2 = (res1 + dr) * h1 / 45; res1 *= 2 * h1 / 45;//вот теперь это интегралы с шагом h2 и h1
    if (FAbs(res2 - res1) > eps)
        return bul_auto(f_sh, a, (a + 2 * h1), eps) + bul_auto(f_sh, (a + 2 * h1), b, eps);
    else
        return res2;
}

#include <iostream>
#include <cmath>
 
double trapezoidal(double (*f)(double), double a, double b, int n) {
    double h = (b-a)/n;
    double sum = 0.5*(f(a) + f(b));
    for (int i = 1; i < n; i++) {
        double x = a + i*h;
        sum += f(x);
    }
    return h*sum;
}
 
 
int main()
{
    auto f = [](double x){return std::pow(std::sin(x),2.)-std::log(x*x);};
    std::cout.precision(12);
    std::cout << trapezoidal(f, 0.1, 10., 1300000);
}

#include <iostream>
#include <ctime>
#include "MyMath/diff.h"
 
using namespace std;
 
double f(double x)
{
    return sin(x) * sin(x) - log(x * x);
}
 
double trapezoidal(double (*f)(double), double a, double b, int n) {
    double h = (b - a) / n;
    double sum = 0.5 * (f(a) + f(b));
    for (int i = 1; i < n; i++) {
        double x = a + i * h;
        sum += f(x);
    }
    return h * sum;
}
 
int main()
{
    time_t start, finish;
    double a = 0.1, b = 10., eps = 1.e-10;
    int n = 10000;
 
    start = clock();
    for (int i = 0; i < n; i++)
        bul_auto(f, a, b, eps);
    finish = clock();
    cout << "\n" << finish - start << "\n\n" << bul_auto(f, a, b, eps) << "\n\n";
    printf("EPS < %e", eps);
 
    start = clock();
    for (int i = 0; i < n; i++)
        trapezoidal(f, a, b, 4096);
    finish = clock();
    cout << "\n" << finish - start << "\n\n" << bul_auto(f, a, b, eps) << "\n\n";
 
    //printf("%e", FAbs(trapezoidal(f, a, b, 4096) - bul_auto(f, a, b, eps)));
 
    printf("EPS = %e", FAbs(trapezoidal(f, a, b, 4096) - trapezoidal(f, a, b, 4096*2)));
 
    return cin.get();
}