/*
//Эксп ПТ
plot "H:/exp_accuracy_T.txt" with lines title "Эксп ПТ Трапец", "H:/exp_accuracy_G.txt" with lines title "Эксп ПТ Гаусс"

//Точность Рунге
plot "H:/Runge_accuracy_T.txt" with lines title "Точность Рунге Трапец","H:/Runge_accuracy_G.txt" with lines title "Точность Рунге Гаус"

//Погреш Рунге
plot "H:/Runge_error_T.txt" with lines title "Погр Рунге Трап", "H:/Runge_error_G.txt" with lines title "Погр Рунге Гаусс"

plot "H:/abs_error_T.txt" with lines title "Абс погр Трапец", "H:/abs_error_G.txt" with lines title "Абс погр Гаусс"
*/

#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)
{
    return x * sin(x);
}
functiontype Func = &Myfunc;
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) / 6) * ((*f)(Array[i + 1].x) + (4 * (*f)(((Array[i].x) + (Array[i + 1].x)) / 2)) + (*f)(Array[i].x));
    }
    return area;
}
double Gauss_formula(Node* Array, functiontype* f, int CountSegments)
{ //Теор. порядок точности = 4


    int i;
    double x1, x2, x3, area = 0;

    double t1, t2, t3 = 0;
    t3 = -sqrt(0.6); t1 = sqrt(0.6); t2 = 0;


    double a1 = 5/9;
    double a2 = 8/9;
    double a3 = 5/9;

    for (i = 0; i < CountSegments; i++)
    {
        x1 = ((Array[i + 1].x + Array[i].x) / 2) + ((Array[i + 1].x - Array[i].x) * t1) / 2;
        x2 = ((Array[i + 1].x + Array[i].x) / 2);
        x3 = ((Array[i + 1].x + Array[i].x) / 2) + ((Array[i + 1].x - Array[i].x) * t3) / 2;
        //area += ((Array[i + 1].x - Array[i].x) / 2) * ((a1 * (*f)(x1) + a2 * (*f)(x2) + a3 * (*f)(x3)));
        area += ((Array[i + 1].x - Array[i].x) / 2) * ((a1 * (*f)(t1) + a2 * (*f)(0) + a3 * (*f)(t3)));
    }

    return area;
}
double orig_integral(double Initial, double End)
{
    return (sin(End) - End * (cos(End))) - (sin(Initial)) - Initial * (cos(Initial));
}

//Заполняет переданный массив эксп порядоком точности для трапец и Гаусса соот-но
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 = 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(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((trapezoid_formula(Array2, f, 2 * CountSegments) - trapezoid_formula(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[0] = R_trap;
    massiv_data[1] = 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(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));

    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[0] = our_accuracy_trap;
    order_massiv[1] = 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 graph_Runge_error(functiontype* f, double Initial, double End, int СountDots) {
    ofstream Trap_error_file("H:/Runge_error_T.txt");
    ofstream Gauss_error_file("H:/Runge_error_G.txt");

    double* Runge_eror_massiv = new double[2];
    for (int n = 1; n < СountDots; n++) {
        Runge_eror_massiv[0] = 0;
        Runge_eror_massiv[1] = 0;
        Runge_Err(f, n, Runge_eror_massiv, Initial, End);
        Trap_error_file << n << " " << Runge_eror_massiv[0] << endl;
        Gauss_error_file << n << " " << Runge_eror_massiv[1] << endl;
    }
}

//Строит график абсолютной погрешности
void graph_abs_error(Node* Array, functiontype* f, double Initial, double End, int СountDots) {
    ofstream Trap_error_file("H:/abs_error_T.txt");
    ofstream Gauss_error_file("H:/abs_error_G.txt");

    for (int n = 1; n < СountDots; n++) {
        Trap_error_file << n << " " << abs(orig_integral(Initial, End) - trapezoid_formula(Array, f, n)) << endl;
        Gauss_error_file << n << " " << abs(orig_integral(Initial, End) - Gauss_formula(Array, f, n)) << endl;
    }
    Trap_error_file.close();
    Gauss_error_file.close();
}


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 = 0;
    Interval.EndNode = 4;
    //Interval.EndNode = 1.5;
    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);

    /*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 << "Абс. погрешность вычисления по формуле Гаусса по 3 узлам:" << 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];
    double* order_massiv = new double[4];
    Runge_Err(&Func, CountSegments, massiv_data, Interval.InitialNode, Interval.EndNode);
    Runge_order_accuracy(CountSegments, &Func, Interval.InitialNode, Interval.EndNode, order_massiv);

    cout << "Оценка погрешности по Рунге для формулы Симпсона:" << endl;
    cout << massiv_data[0] << endl << endl;
    cout << "Оценка погрешности по Рунге для формулы Гаусса по 3 узлам:" << endl;
    cout << massiv_data[1] << endl << endl;

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

    /* -----------------Порядок точности на основе правила Рунге -----------------------------*/
    cout << "Оценка точности на основе правила Рунге для формулы Симпсона:" << endl;
    cout << order_massiv[0] << endl << endl;
    cout << "Оценка точности на основе правила Рунге для для формулы Гаусса по 3 узлам:" << endl;
    cout << order_massiv[1] << endl << endl;

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

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

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

    //График абсолютной погрешности
    graph_abs_error(ArrayUniformNodes, &Func, Interval.InitialNode, Interval.EndNode, СountDots);


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