Построить таблицу приближенных значений, используя составную квадраторную формулу центральных прямоугольников

#include <iostream>
#include <math.h>
#include <iomanip>
#include <string>
#define double long double
#define sigma 0.577215
 
using namespace std;
 
// Печать матриц виде таблиц
int PrintTable2(double **mas, int size, string title = "", int width1 = 30, int width2 = 30)
{
    cout << title << endl;
    for(int i = 0; i < (width1 + width2) + 9; i++)
        cout << "=";
    cout << endl;
 
    for (int i = 0; i < size; i ++)
    {
        cout.unsetf (ios::fixed);
        cout << setprecision(2) << setw(2) << i << " ||" << setprecision(10) << setw(width1) << mas[i][0] ;
        cout.setf(ios::fixed);
        cout << " |" << setw(width2) << mas[i][1] << " |"  << endl;
    }
 
    for(int i = 0; i < (width1 + width2) + 9; i++)
        cout << "=";
    cout << endl;
    return true;
}
int PrintTable4(double **mas, int size, string title = "", int width1 = 30, int width2 = 30, int width3 = 30, int width4 = 30)
{
    cout << title << endl;
    for(int i = 0; i < (width1 + width2 + width3 + width4) + 13; i++)
        cout << "=";
    cout << endl;
 
    for (int i = 0; i < size; i ++)
    {
        cout.unsetf (ios::fixed);
        cout << setprecision(2) << setw(2) << i << " ||" << setprecision(10) << setw(width1) << mas[i][0] ;
        cout.setf(ios::fixed);
        cout << " |" << setw(width2) << mas[i][1] << " |" << setw(width3) << mas[i][2] << " |" << setw(width4) << setprecision(10)  << mas[i][3] << " |" << endl;
    }
 
    for(int i = 0; i < (width1 + width2 + width3 + width4) + 13; i++)
        cout << "=";
    cout << endl;
    return true;
}
double Factorial(int n)
{
    return !n ? 1 : n * Factorial(n - 1);
}
double Ci(double x,double eps = 0.000001)
{
    double Ci = 0;
    double Cipred = 1 + eps;
    int minus = -1;
    int n = 1;
    double buf;
    while((fabs(Ci - Cipred) > eps))
    {
        Cipred = Ci;
        buf = minus * pow(x, (2 * n))/2;
        buf /= n;
        buf /= Factorial(2*n);
        Ci += buf;
        minus = minus * (-1);
        n++;
    }
    return Ci + log(x) + sigma;
}
 
double CiIntegral(double t)
{
    return (cos(t) - 1)/t;
}
 
double Lagrange(double **mas, int n, double x)
{
    double sum = 0;
    double proizved = 1;
    for(int i = 0; i < n; i++)
    {
        proizved = 1;
        for(int j = 0; j < n; j++)
        {
            if(i == j) continue;
 
            proizved *= (x - mas[j][0]) / (mas[i][0] - mas[j][0]);
        }
        sum += mas[i][1] * proizved;
    }
    return sum;
}
 
double Max(double **mas, int size)
{
    double max = mas[0][1];
    for(int i = 1; i < size; i++)
    {
        if(mas[i][1] > max) max = mas[i][1];
    }
    return max;
}
 
double Min(double **mas, int size)
{
    double min = mas[0][1];
    for(int i = 1; i < size; i++)
    {
        if(mas[i][1] < min) min = mas[i][1];
    }
    return min;
}
//составляем составную квадраторную формулу центральных прямоугольников
double Integral(double a,double b, int fragments = 100)
{
    double step = (b - a) / fragments;
    double sum = 0;
    for(int i = 1; i < fragments; i++)
    {
        sum += step * CiIntegral( ((a + (i - 1) * step) + (a + (i) * step)) / 2);
    }
    return sum;
}
// Нахождение производной
double Derivative(double x, double h = 0.001)
{
    return (-Ci(x + 2*h) + 8*Ci(x + h) - 8*Ci(x - h) + Ci(x - 2*h))/(12 * h);
}
// Нахождение обратной функции
double Find_root(double b, double y, double eps)
{
    double a;
    double gz, gzu;
    for(;;)
    {
        gz = Ci(b) - y;
        gzu = Ci(b - gz/Derivative(b)) - y;
        a = b - gz/Derivative(b)  +  gz * gzu  /  Derivative(b) * (gzu - gz);
 
        if (fabs(b - a) <= eps) break;
        b = a;
    }
    return b;
}
 
int main()
{
    setlocale(LC_ALL, "Russian");
    double a = 0.4;
    double b = 4.0;
    double step = 0.1999999999999999;
    double step2 = 0.1;
    double eps = 0.000001;
    double integration_eps = 0.00001;
    double accuracy = 0; //Макс погрешность интерполирования
 
    // Задаём массив для хранения значений интегрального косинуса
    double **cosinus_table = new double * [static_cast< int >((b-a) / step)];
    int cosinus_table_size = 0;
 
    // Задаём массив для хранения начений погрешностей
    //значений вычесленных с помощью полинома Лагранжа
    double **table_accuracy = new double * [static_cast< int >((b-a) / step2)];
    int table_accuracy_size = 0;
 
    // Заполняем таблицу значениями интегрального косинуса
    int n = 0;
    for(double i = a; i <= b; i += step)
    {
        cosinus_table[n] = new double[2];
        cosinus_table[n][0] = i;
        cosinus_table[n][1] = Ci(i, eps);
        n++;
        cosinus_table_size++;
    }
 
    PrintTable2(cosinus_table, cosinus_table_size, "Значения интегрального косинуса:", 4, 13);
    //////////////////КОНЕЦ 1 ЗАДАНИЯ
 
 
    n = 0;
    for(double i = 0.4; i <= 4; i += 0.1)
    {
        table_accuracy[n] = new double[4];
        table_accuracy[n][0] = i;
        table_accuracy[n][1] = Ci(table_accuracy[n][0], eps);
        table_accuracy[n][2] = Lagrange(cosinus_table, cosinus_table_size, i);
        table_accuracy[n][3] = fabs(table_accuracy[n][2] - table_accuracy[n][1]);
        if(table_accuracy[n][3] > accuracy)
            accuracy = table_accuracy[n][3];
        n++;
        table_accuracy_size++;
    }
    PrintTable4(table_accuracy, table_accuracy_size, "\nСравнивание значений функции со значением полинома\nЛагранжа, с удвоенным шагом (0,1):" , 4, 14, 14, 14);
    cout << "Максимальная погрешность интерполирования = " << setprecision(10) << fixed << accuracy << "\n\n";
    //////////////////КОНЕЦ 2 ЗАДАНИЯ
 
    double **table_integral = new double * [static_cast< int >((b-a) / step)];
    int table_integral_size = 0;
    n = 0;
    for(double i = a; i <= b; i += step)
    {
        table_integral[n] = new double[2];
        table_integral[n][0] = i;
 
        double  buf = 0;
        double  buf_pred = integration_eps + 1;
        int fragments = 50;
        for(;;)
        {
            if(fabs(buf_pred - buf) < integration_eps)
            {
                table_integral[n][1] = buf + log(i) + sigma;
                break;
            }
            buf_pred = buf;
            buf = Integral(0, i, fragments);
            fragments *= 2;
        }
        table_integral_size++;
        n++;
    }
    PrintTable2(table_integral, table_integral_size, "Значения вычесленные с помощью вычисления интеграла:", 4, 13);
    //////////////////КОНЕЦ 3 ЗАДАНИЯ
    // Задаём массив для хранения обратной таблицы
    double **table_reverse = new double * [static_cast< int >((b-a) / step)];
    int table_reverse_size = 0;
    n = 0;
    double reverse_step = fabs((Max(cosinus_table, cosinus_table_size) - Min(cosinus_table, cosinus_table_size)) / (cosinus_table_size - 1) );
    for(double i = -0.178682515480129; i <= Max(cosinus_table, cosinus_table_size) ; i += reverse_step)
    {
        table_reverse[n] = new double[2];
        table_reverse[n][0] = i;
        table_reverse[n][1] = Find_root(3, i, 0.0001);
        table_reverse_size++;
        n++;
    }
    PrintTable2(table_reverse, table_reverse_size, "\n\nОбратная таблица к Сi(x):", 14, 14);
    /////////////////////////////// КОНЕЦ 4 задания
return 0;
}