#include <iostream>
#include <math.h>
#include <cmath>
#include <vector>
#include <Windows.h>
#include <stdlib.h>
#include <fstream>
#include <string.h>
#include <time.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)
{

    if (sin(x) < 0) return -sin(x);
    else return sin(x);

}
functiontype Func = &Myfunc;
double trapeciya(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
    //clock_t time;
    //time = clock();


    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;
    }
    //time = clock() - time;
    //cout << (double)time << "CLOCKS_PER_SEC"; //время выполнения "каких-то действий"
    return area;

}
double orig_integral(double Initial, double End)
{

    return (2 - cos(Initial) - cos(End));
    //return(cos(Initial) - cos(End));  //для обоих краев больших 0
    //return (cos(End) - cos(Initial)); // для обоих краев меньших 0
}

//Заполняет переданный массив эксп порядоком точности для трапец и Гаусса соот-но
void exp_order_accuracy(int CountSegments, functiontype* f, double Initial, double End, double* massiv_exp_accuracy) {

    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 = trapeciya(ArrayNodes, f, CountSegments);
    trapezoid_square_exp2 = trapeciya(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(abs(R_G / R_G2)) / log(2);
    exp_p_T = log(abs(R_T / R_T2)) / log(2);

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

}

//Заполняет переданный массив погрешностью полученной по правилу Рунге для трапец и Гаусса соот-но
void Runge_Err(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;

    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 = abs((trapeciya(Array2, f, 2 * CountSegments) - trapeciya(Array, f, CountSegments)))* pow(2, p_trap);
    numerator_Gauss = abs((Gauss_formula(Array2, f, 2 * CountSegments) - Gauss_formula(Array, f, CountSegments)))* pow(2, p_Gauss);

    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;

    massiv_data[2] = R_trap;
    massiv_data[3] = R_Gauss;
}

//Заполняет переданный массив порядком точности на основе правила Рунге для трапец и Гаусса соот-но
void Runge_order_accuracy(int CountSegments, functiontype* f, double Initial, double End, double* order_massiv) {

    double our_numerator_trap, our_numerator_Gauss, our_denominator_Gauss, our_denominator_trap, our_accuracy_trap, our_accuracy_Gauss;

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

    our_numerator_trap = abs(trapeciya(Array2, f, 2 * CountSegments) - trapeciya(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(trapeciya(Array2, f, 2 * CountSegments) - trapeciya(Array4, f, 4 * CountSegments));

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

    order_massiv[4] = our_accuracy_trap;
    order_massiv[5] = our_accuracy_Gauss;

}

//Строит график эксп точности
void graphic_exp_acc(functiontype* f, double Initial, double End, int СountDots) {
    ofstream Trap_File("H:/exp_accuracy_T.txt");
    ofstream Gauss_File("H:/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;
        exp_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 graphic_Runge_acc(functiontype* f, double Initial, double End, int СountDots) {
    ofstream Trap_File("H:/Runge_accuracy_T.txt");
    ofstream Gauss_File("H:/Runge_accuracy_G.txt");
    double* massiv_Runge_accuracy = new double[2];
    for (int n = 1; n < СountDots; n++) {
        massiv_Runge_accuracy[0] = 0;
        massiv_Runge_accuracy[1] = 0;
        Runge_order_accuracy(n, f, Initial, End, massiv_Runge_accuracy);
        Trap_File << n << " " << massiv_Runge_accuracy[0] << endl;
        Gauss_File << n << " " << massiv_Runge_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;
}


int main()
{
    setlocale(LC_ALL, "RUS");
    functiontype Func = &Myfunc;
    Interval Interval; Interval.InitialNode = -1; Interval.EndNode = 2;
    int CountSegments, СountDots = 1000;
    cout << "Введите число интервалов разбиения: " << endl; cin >> CountSegments; cout << endl;
    int CountNodes = CountSegments + 1;
    Node* ArrayUniformNodes = new Node[CountNodes];
    ValueUniformTable(&Func, ArrayUniformNodes, Interval.InitialNode, Interval.EndNode, CountNodes);
    double Value_trapeciya = trapeciya(ArrayUniformNodes, &Func, CountSegments);
    double Value_Gauss_formula = Gauss_formula(ArrayUniformNodes, &Func, CountSegments);
    cout << "полученные значения интеграла по формуле трапеции и Гаусса с  2 точками:" << endl;
    cout << Value_trapeciya << endl;
    cout << Value_Gauss_formula << endl << endl;

    cout << "действительное значение интеграла:" << endl;
    double Value_orig_integral = orig_integral(Interval.InitialNode, Interval.EndNode);
    cout << Value_orig_integral << endl << endl<<endl;

    cout << "сравнение значения интеграла со значением вычисленным по правилу трапеции и Гаусса соответственно:" << endl;
    cout << abs(Value_orig_integral - Value_trapeciya) << endl;
    cout << abs(Value_orig_integral - Value_Gauss_formula) << endl << endl;

    double* Massiv_of_Data = new double[6];
    for (int i = 0; i < 6; i++)
        Massiv_of_Data[i] = 0;

    Runge_order_accuracy(CountSegments, &Func, Interval.InitialNode, Interval.EndNode, Massiv_of_Data);
    exp_order_accuracy(CountSegments, &Func, Interval.InitialNode, Interval.EndNode, Massiv_of_Data);
    Runge_Err(&Func, CountSegments, Massiv_of_Data, Interval.InitialNode, Interval.EndNode);

    cout << "Экспериментальный порядок точности метода трапеции и Гаусса:" << endl;
    cout << Massiv_of_Data[0] << "   " << Massiv_of_Data[1] << endl << endl;

    cout << "Практическая оценка погрешности на основе правила Рунге для метода трапеции и Гаусса:" << endl;
    cout << Massiv_of_Data[2] << "    " << Massiv_of_Data[3] << endl << endl;

    cout << "Оценка порядка точности на основе правила Рунге для метода трапеции и Гаусса:" << endl;
    cout << Massiv_of_Data[4] << "    " << Massiv_of_Data[5] << endl << endl;

    graphic_exp_acc(&Func, Interval.InitialNode, Interval.EndNode, СountDots);
    graphic_Runge_acc(&Func, Interval.InitialNode, Interval.EndNode, СountDots);

    system("pause");
    return 0;
}