#include <iostream>
#include <math.h>
#include <cmath>
#include <vector>
#include <Windows.h>
#include <stdlib.h>
#include <fstream>
#include <string.h>
using namespace std;
const double PI = 3.1415926535897932384626433832795;
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;

int Factorial(int n)
{ // Факториал
    int x = 1;
    for (int i = 1; i <= n; i++) {
        x *= i;
    }
    return x;
}

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

    }
}

//Возвращает число, возведенное в cтепень i
double get_degree(double x, int degree)
{
    int i;
    double y = x;
    if (degree == 0)
        return 1;

    for (i = 0; i < degree - 1; i++)
        y *= x;
    return y;
}

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

void orig_table_in_file(Node* Array, int Count)
{ // Функция берет таблицу иксов и игриков, количество точек в которых считаем знач полинома,
    // По итогу имеем файл с 2 столбацами: х и у
    int k;
    double y_value, x_value;

    ofstream fout("D:/Original_graphic.txt");

    for (k = 0; k < Count; k++) { //n иксов подставляев в полином степени n-1
        x_value = Array[k].x;
        fout << x_value << " ";
        y_value = Array[k].y;
        fout << y_value << endl;
    }

    fout.close();
}

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 Gause_formula(Node* Array, functiontype* f, int CountSegments) { //Теор. порядок точности =

/*порядок точности составной квадратурной формулы Гаусса по 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 ;
}

double get_C(functiontype* f, int CountSegments2, int CountSegments, Node* Array2, Node* Array, string method_name) {
    int p;
    double numerator, denominator;
    functiontype Func = &Myfunc;
    double step2 = 2 / (CountSegments2);

    if (method_name == "trap") {
        p = 2;
        numerator = (trapezoid_formula(Array2, &Func, CountSegments2) - trapezoid_formula(Array, &Func, CountSegments)) * pow(2, p);
    }
    else {
        p = 4;
        numerator = (Gause_formula(Array2, &Func, CountSegments2) - Gause_formula(Array, &Func, CountSegments)) * pow(2, p);
    }

    denominator = pow(step2, p) * (pow(2, p) - 1);

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

    return numerator / denominator;

}


int main()
{
    setlocale(LC_ALL, "RUS");
    Interval Interval;
    int CountSegments, MyNodes = 5000;
    double C;
    functiontype Func = &Myfunc;
    cout << "Enter the n: " << endl;
    cin >> CountSegments;
    double CountNodes = CountSegments + 1;
    Interval.InitialNode = 0;
    Interval.EndNode = 2;

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


    //Реализация правила Рунге
    Node* Array2Nodes = new Node[int(CountNodes)];
    Node* ArrayNodes = new Node[int(CountNodes/2)];
    ValueUniformTable(&Func, Array2Nodes, Interval.InitialNode, Interval.EndNode, CountNodes);
    ValueUniformTable(&Func, ArrayNodes, Interval.InitialNode, Interval.EndNode, int(CountNodes/2));
    C = get_C(&Func, CountSegments, int(CountNodes / 2), Array2Nodes, ArrayNodes, "trap");
    cout << "C = " << C << endl;

    /*cout << endl;
    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 Gause_square = Gause_formula(ArrayUniformNodes, &Func, CountSegments);
    cout << Gause_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 - Gause_square) << endl << endl;

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





    system("pause");
    return 0;
}