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

{
    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;
    }
}
double Myfunc(double x)
{
    return (x * (sqrt(4 - x * x)));
}
functiontype Func = &Myfunc;
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;
}
double orig_integral(double Initial, double End)
{
    return (pow((4 - Initial * Initial), 1.5) - pow((4 - End * End), 1.5)) / 3;
}
void get_order_accuracy(int CountSegments, functiontype* f, double Initial, double End, double* massiv_exp_accuracy) {
    /*ofstream Trap_File("D:/exp_accuracy_T.txt");
    ofstream Gauss_File("D:/exp_accuracy_T.txt");*/

    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;


    Node* Array2Nodes = new Node[2 * CountSegments + 1];
    Node* ArrayNodes = new Node[CountSegments + 1];
    ValueUniformTable(f, Array2Nodes, Initial, End, 2 * CountSegments + 1);
    ValueUniformTable(f, ArrayNodes, Initial, End, CountSegments + 1);
    Gauss_square_exp = Gauss_formula(ArrayNodes, f, CountSegments);
    Gauss_square_exp2 = Gauss_formula(Array2Nodes, f, 2 * CountSegments);
    trapezoid_square_exp = trapezoid_formula(ArrayNodes, f, CountSegments);
    trapezoid_square_exp2 = trapezoid_formula(Array2Nodes, f, 2 * CountSegments);
    orig_square_exp = orig_integral(Initial, End);

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

    /*Trap_File << CountSegments << " " << exp_p_T;
    Gauss_File << CountSegments << " " << exp_p_G;*/

    massiv_exp_accuracy[0] = exp_p_T;
    massiv_exp_accuracy[1] = exp_p_G;

}
void graphic_exp(functiontype* f, double Initial, double End, int СountDots) {
    ofstream Trap_File("D:/exp_accuracy_T.txt");
    ofstream Gauss_File("D:/exp_accuracy_G.txt");
    double* massiv_exp_accuracy = new double[2];
    for (int n = 1; n < СountDots; n++) {
        massiv_exp_accuracy[0] = 0;
        massiv_exp_accuracy[1] = 0;
        get_order_accuracy(n, f, Initial, End, massiv_exp_accuracy);
        Trap_File << n << " " << massiv_exp_accuracy[0] << endl;
        Gauss_File << n << " " << massiv_exp_accuracy[1] << endl;
    }

}
void PrintNodes(Node* Array, int CountSegments)
{
    int i;
    for (i = 0; i < CountSegments - 1; i++)
        cout << "(" << (Array[i].x) << ":" << (Array[i + 1].x) << ")" << endl;
}
void get_mas_pract_error(functiontype* f, int CountSegments, double*& massiv_data, double Initial, double End)
{
    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;

    Node* Array4 = new Node[4 * CountSegments + 1];
    Node* Array2 = new Node[2 * CountSegments + 1];
    Node* Array = new Node[CountSegments + 1];
    ValueUniformTable(f, Array4, Initial, End, 4 * CountSegments + 1);
    ValueUniformTable(f, Array2, Initial, End, 2 * CountSegments + 1);
    ValueUniformTable(f, Array, Initial, End, CountSegments + 1);

    numerator_trap = (trapezoid_formula(Array2, f, CountSegments) - trapezoid_formula(Array, f, CountSegments)) * pow(2, p_trap);
    numerator_Gauss = (Gauss_formula(Array2, f, 2 * CountSegments) - Gauss_formula(Array, f, CountSegments)) * pow(2, p_Gauss);

    our_numerator_trap = abs(trapezoid_formula(Array2, f, 2 * CountSegments) - trapezoid_formula(Array, f, CountSegments));
    our_numerator_Gauss = abs(Gauss_formula(Array2, f, 2 * CountSegments) - Gauss_formula(Array, f, CountSegments));
    our_denominator_Gauss = abs(Gauss_formula(Array2, f, 2 * CountSegments) - Gauss_formula(Array4, f, 4 * CountSegments));
    our_denominator_trap = abs(trapezoid_formula(Array2, f, 2 * CountSegments) - trapezoid_formula(Array4, f, 4 * CountSegments));

    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");
    functiontype Func = &Myfunc;
    Interval Interval; Interval.InitialNode = 0; Interval.EndNode = 2;
    int CountSegments, СountDots = 500;
    cout << "Введите число интервалов разбиения: " << endl; cin >> CountSegments; cout << endl;
    double CountNodes = CountSegments + 1;

    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];
    get_mas_pract_error(&Func, CountSegments, massiv_data, Interval.InitialNode, Interval.EndNode);

    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;
    /* -----------------------Конец Рунге-----------------------------------*/


    /* -----------------Экспериментальный порядок точности   -----------------------------*/
    double* massiv_exp_accuracy = new double[2];
    get_order_accuracy(CountSegments, &Func, Interval.InitialNode, Interval.EndNode, massiv_exp_accuracy);
    cout << "Эксп. Порядок точности для Гаусса и Трапеций соответственно:" << endl;
    cout << massiv_exp_accuracy[1] << " " << massiv_exp_accuracy[0] << endl << endl;

    //График
    graphic_exp(&Func, Interval.InitialNode, Interval.EndNode, СountDots);
    /* -------------------------   Конец --------------------------------------*/

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