#include <iostream>
#include <math.h>
#include <cmath>
#include <vector>
#include <Windows.h>
#include <stdlib.h>
#include <fstream>
#include <string.h>
using namespace std;

typedef double(*functiontype)(double x);

typedef struct Node
{
    double x, y;
} Node;

typedef double(*method)(double x, Node* Array, int Count, double* DD_massiv);

typedef struct Interval
{
    double InitialNode, EndNode;
} Interval;



void ValueUniformTable(functiontype* f, Node* Array, double Initial, double End,
    int CountNodes) // Точек на 1 больше чем отрезков, массив
// хранит все точки (их пары образуют
// отрезки)
{ // Создание равномерной таблицы значений
    double step = abs(Initial - End) / (CountNodes - 1);
    Array[0].x = Initial;
    Array[0].y = (*f)(Array[0].x);
    for (int i = 1; i < CountNodes; i++)
    {
        Array[i].x = Array[i - 1].x + step;
        Array[i].y = (*f)(Array[i].x);
    }
}

double Myfunc(double x)
{
    return (x * (sqrt(4 - x * x)));
}


double trapezoid_formula(Node* Array, functiontype* f, int CountSegments)
{ //Теор. порядок точности = 2
    int i;
    double area = 0;

    for (i = 0; i < CountSegments; i++)
        area += (Array[i + 1].x - Array[i].x) * ((*f)(Array[i + 1].x) + (*f)(Array[i].x)) / 2;
    return area;
}

double Gauss_formula(Node* Array, functiontype* f, int CountSegments)
{ //Теор. порядок точности = 4

    /*порядок точности составной квадратурной формулы Гаусса по n узлам равен 2n.
  стр 216
  квадр формула Гаусс точно на полиномах степени 2n-1 стр 209
  211 стр формула погр
  Алгебраическая степень точности квадратурной формулы, построенной по n узлам,
  не может превосходить 2n−1.
  Как следствие, если производные подинтегральной функции не рас-
  тут «слишком быстро» с ростом их порядка, то при увеличении числа

  узлов и гладкости интегрируемой функции порядок точности квадра-
  турных формул Гаусса может быть сделан сколь угодно высоким.
  202 стр теорема*/
    int i;
    double x1, x2, area = 0;
    for (i = 0; i < CountSegments; i++)
    {
        x1 = (Array[i + 1].x + Array[i].x) / 2 + (Array[i + 1].x - Array[i].x) * sqrt(3) / 6;
        x2 = (Array[i + 1].x + Array[i].x) / 2 - (Array[i + 1].x - Array[i].x) * sqrt(3) / 6;
        area += (Array[i + 1].x - Array[i].x) * ((*f)(x1) + (*f)(x2)) / 2;
    }

    return area;
}

void PrintNodes(Node* Array, int CountSegments)
{
    int i;
    for (i = 0; i < CountSegments - 1; i++)
        cout << "(" << (Array[i].x) << ":" << (Array[i + 1].x) << ")" << endl;
}

double orig_integral(double Initial, double End)
{
    return (pow((4 - Initial * Initial), 1.5) - pow((4 - End * End), 1.5)) / 3;
}

void get_mas_pract_error(functiontype* f, int CountSegments2, int CountSegments, Node* Array2,
    Node* Array, double*& massiv_errors)
{
    int p_trap = 2, p_Gauss = 4;
    double numerator_trap, numerator_Gauss, denominator_Gauss, denominator_Trap, x, practical_eror_trap,
        practical_eror_trap2, practical_eror_Gauss;
    double practical_eror_Gauss2, C_trap, C_Gauss;
    functiontype Func = &Myfunc;
    double z = CountSegments2, y = CountSegments;
    double  step = 2 / y; //Пришлось так сделать чтобы дабл на
    //дабл делился, а то если на инт то 0
    //получается

    numerator_trap = (trapezoid_formula(Array2, &Func, CountSegments2)
        - trapezoid_formula(Array, &Func, CountSegments))
        * pow(2, p_trap);
    numerator_Gauss = (Gauss_formula(Array2, &Func, CountSegments2)
        - Gauss_formula(Array, &Func, CountSegments))
        * pow(2, p_Gauss);
    denominator_Gauss = pow(step, p_Gauss) * (pow(2, p_Gauss) - 1);
    denominator_Trap = pow(step, p_trap) * (pow(2, p_trap) - 1);
    C_trap = numerator_trap / denominator_Trap;
    C_Gauss = numerator_Gauss / denominator_Gauss;
    //cout << C_trap << "   " << C_Gauss;

    practical_eror_trap = C_trap * pow(step, p_trap);
    //practical_eror_trap2 = C_trap * pow(step2, p_trap);
    practical_eror_Gauss = C_Gauss * pow(step, p_Gauss);
    //practical_eror_Gauss2 = C_Gauss * pow(step2, p_Gauss);

    massiv_errors[0] = practical_eror_trap;
    // massiv_errors[1] = practical_eror_trap2;
    massiv_errors[1] = practical_eror_Gauss;
    // massiv_errors[3] = practical_eror_Gauss2;

    // cout <<  " numerator: " << numerator <<   " denominator: " <<  denominator
    // << endl;
}

int main()
{
    setlocale(LC_ALL, "RUS");
    Interval Interval;
    int CountSegments, MyNodes = 5000;
    double C_Trap, C_Gauss;
    double orig_square_exp, R_G, R_G2, R_T, R_T2, exp_p_G, exp_p_T;
    double Gauss_square_exp, Gauss_square_exp2, trapezoid_square_exp, trapezoid_square_exp2;
    functiontype Func = &Myfunc;
    cout << "Введите число интервалов разбиения: " << endl;
    cin >> CountSegments;
    cout << endl;
    double CountNodes = CountSegments + 1;
    Interval.InitialNode = 0;
    Interval.EndNode = 2;

    Node* ArrayUniformNodes = new Node[CountNodes];
    ValueUniformTable(&Func, ArrayUniformNodes, Interval.InitialNode, Interval.EndNode, CountNodes);

    cout << "Разбиение интервала на равные отрезки интегрирования:" << endl;
    PrintNodes(ArrayUniformNodes, CountNodes);
    cout << endl;

    cout << "Значение интеграла через формулу трапеций:" << endl;
    double trapezoid_square = trapezoid_formula(ArrayUniformNodes, &Func, CountSegments);
    cout << trapezoid_square << endl << endl;

    cout << "Значение интеграла через формулу Гаусса по 2 узлам:" << endl;
    double Gauss_square = Gauss_formula(ArrayUniformNodes, &Func, CountSegments);
    cout << Gauss_square << endl << endl;

    cout << "Значение интеграла:" << endl;
    double orig_square = orig_integral(Interval.InitialNode, Interval.EndNode);
    cout << orig_square << endl << endl;

    cout << "Разность между значениям интеграла и значением по формуле Гаусса:" << endl;
    cout << abs(orig_square - Gauss_square) << endl << endl;

    cout << "Разность между значениям интеграла и значением по формуле Трапеций:" << endl;
    cout << abs(orig_square - trapezoid_square) << endl << endl;

    /* -----------------Реализация правила Рунге-----------------------------*/
    double* massiv_errors = new double[2];
    Node* Array2Nodes = new Node[2*CountSegments+1];
    Node* ArrayNodes = new Node[CountNodes];
    ValueUniformTable(&Func, Array2Nodes, Interval.InitialNode, Interval.EndNode, 2 * CountSegments + 1);
    ValueUniformTable(&Func, ArrayNodes, Interval.InitialNode, Interval.EndNode, CountNodes);
    get_mas_pract_error(&Func, 2 * CountSegments, CountSegments, Array2Nodes, ArrayNodes, massiv_errors);
    cout << "Практ. оценка погр. для разбиения на " << CountSegments << " частей по формуле трапеций:" << endl;
    cout << massiv_errors[0] << endl << endl;
    cout << "Практ. оценка погр. для разбиения на " << CountSegments << " частей по формуле Гаусса:" << endl;
    cout << massiv_errors[1] << endl << endl;

    // cout << "C_Trap2 = " << massiv_errors[1] << endl << endl;
    // cout << "C_Gauss2 = " << massiv_errors[3] << endl << endl;
    /* -----------------------Конец Рунге-----------------------------------*/

    /* -----------------Экспериментальный порядок точности   -----------------------------*/

    Gauss_square_exp = Gauss_formula(ArrayNodes, &Func, CountSegments);
    Gauss_square_exp2 = Gauss_formula(Array2Nodes, &Func, 2 * CountSegments);
    trapezoid_square_exp = trapezoid_formula(ArrayNodes, &Func, CountSegments);
    trapezoid_square_exp2 = trapezoid_formula(Array2Nodes, &Func, 2 * CountSegments);

    orig_square_exp = orig_integral(Interval.InitialNode, Interval.EndNode);

    R_G = abs(orig_square_exp - Gauss_square_exp);
    R_G2 = abs(orig_square_exp - Gauss_square_exp2);
    R_T = abs(orig_square_exp - trapezoid_square_exp);
    R_T2 = abs(orig_square_exp - trapezoid_square_exp2);

    exp_p_G = log(R_G / R_G2) / log(2) ;
    exp_p_T = log(R_T / R_T2) / log(2) ;

    cout << "Эксперим. порядок точности для формул Гаусса и Трапеций соответственно:" << endl;
    cout << exp_p_G << " " << exp_p_T << endl << endl;
    /* -------------------------   Конец --------------------------------------*/

    cout << endl;
    system("pause");
    return 0;
}