#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-1; 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-1; 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 CountSegments4, Node* Array4, int CountSegments2, int CountSegments, Node* Array2, Node* Array, double*& massiv_data)
{
    functiontype Func = &Myfunc;
    int p_trap = 2, p_Gauss = 4;
    double numerator_trap, numerator_Gauss, denominator_Gauss, denominator_Trap, R_trap, R_Gauss;
    double our_numerator_trap, our_numerator_Gauss, our_denominator_Gauss, our_denominator_trap, our_accuracy_trap, our_accuracy_Gauss;
    double y = CountSegments, step = 2 / y;

    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);

    our_numerator_trap = abs(trapezoid_formula(Array2, &Func, CountSegments2) - trapezoid_formula(Array, &Func, CountSegments));
    our_numerator_Gauss = abs(Gauss_formula(Array2, &Func, CountSegments2) - Gauss_formula(Array, &Func, CountSegments));
    our_denominator_Gauss = abs(Gauss_formula(Array2, &Func, CountSegments2) - Gauss_formula(Array4, &Func, CountSegments4));
    our_denominator_trap = abs(trapezoid_formula(Array2, &Func, CountSegments2) - trapezoid_formula(Array4, &Func, CountSegments4));

    denominator_Gauss = (pow(2, p_Gauss) - 1);
    denominator_Trap = (pow(2, p_trap) - 1);
    R_trap = numerator_trap / denominator_Trap;
    R_Gauss = numerator_Gauss / denominator_Gauss;

    our_accuracy_trap = log(our_numerator_trap / our_denominator_trap) / log(2);
    our_accuracy_Gauss = log(our_numerator_Gauss / our_denominator_Gauss) / log(2);

    massiv_data[0] = R_trap;
    massiv_data[1] = R_Gauss;
    massiv_data[2] = our_accuracy_trap;
    massiv_data[3] = our_accuracy_Gauss;

}

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 << "Гаусс:" << 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_data = new double[4];
    Node* Array4Nodes = new Node[4*CountSegments+1];
    Node* Array2Nodes = new Node[2*CountSegments+1];
    Node* ArrayNodes = new Node[CountSegments+1];
    ValueUniformTable(&Func, Array4Nodes, Interval.InitialNode, Interval.EndNode, 4*CountSegments+1);
    ValueUniformTable(&Func, Array2Nodes, Interval.InitialNode, Interval.EndNode, 2*CountSegments+1);
    ValueUniformTable(&Func, ArrayNodes, Interval.InitialNode, Interval.EndNode, CountSegments+1);
    get_mas_pract_error(&Func, 4*CountSegments, Array4Nodes,2*CountSegments, CountSegments, Array2Nodes, ArrayNodes, massiv_data);
    //PrintNodes(Array2Nodes, 2 * CountSegments + 1);

    cout << "Погрешность по Рунге для " << CountSegments << " частей по формуле трапеций:" << endl;
    cout << massiv_data[0] << endl << endl;
    cout << "Погрешность по Рунге для " << CountSegments << " частей по формуле Гаусса:" << endl;
    cout << massiv_data[1] << endl << endl;
    cout << "Точность по Антону " << CountSegments << " Трапеция" << endl;
    cout << massiv_data[2] << endl << endl;
    cout << "Точность по Антону " << CountSegments << " Гаусс" << endl;
    cout << massiv_data[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;
}