#include <iostream>
#include <math.h>
#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;
}

//Возвращает число, возведенное в 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 / (x*x+1));
}

void ValueUniformTable(functiontype* f, Node* Array, double Initial, double End, int Count) // Передаем массив в который все будет записано
{ // Создание равномерной таблицы значений
    double step = abs(Initial - End) / (Count - 1);
    Array[0].x = Initial;
    Array[0].y = (*f)(Array[0].x);
    for (int i = 1; i < Count; i++) {
        Array[i].x = Array[i - 1].x + step;
        Array[i].y = (*f)(Array[i].x);
    }
}

void ValueChebyshevTable(functiontype* f, Node* Array, double Initial, double End, int Count)
{ // Создание таблицы Чебышевских значений
    for (int i = 0; i < Count; i++) {
        Array[i].x = ((End + Initial) / 2)
            + ((End - Initial) / 2) * cos(((2 * i + 1) * PI) / (2 * (Count + 1)));
        Array[i].y = (*f)(Array[i].x);
    }
}

double* DividedDifference(int Count, Node* Array) {
    double* tmp_mas = new double[Count];
    double* dd = new double[Count];
    for (int i = 0; i < Count; i++) {
        tmp_mas[i] = Array[i].y;
        dd[i] = 0.0;
    }
    dd[0] = Array[0].y;
    for (int i = 1; i < Count; i++) {
        for (int j = 0; j < Count - i; j++) {
            tmp_mas[j] = (tmp_mas[j + 1] - tmp_mas[j]) / (Array[i + j].x - Array[j].x);
        }
        dd[i] = tmp_mas[0];
    }
    return dd;
}

double W(double x, int n, Node* Array)
{ // Полином вида: (x - x1) * (x - x2) * ... * (x - xn)
    double w = 1;
    for (int i = 0; i < n; i++) {
        w *= (x - Array[i].x);
    }
    return w;
}

void get_koef_polynom(double* mas, double* res, int CountBrackets) // k - кол во скобок для раскрытия (сколько надо, напр в ньютоне нужно не все), CountBrackets - максимальное кол-во скобок
{
    int i;
    double* anx = new double[CountBrackets + 1];
    double* ann = new double[CountBrackets + 1];
    res[0] = (mas[0]);
    res[1] = 1;
    int j;
    for (j = 0; j <= CountBrackets; j++) {
        anx[j] = 0.0;
        ann[j] = 0.0;
    }


    for (i = 1; i < CountBrackets; i++)
        for (j = 0; j <= i + 1; j++)
        {
            anx[j + 1] = res[j];
            ann[j] = (res[j] * (mas[i]));
            res[j] = (anx[j] + ann[j]);
        }
}

// Полином Лагранжа в точке
double PolynomLG_dot(double x, Node* Array, int Count, double *DD_massiv)
{
    double Polynom = 0;
    for (int i = 0; i < Count; i++) {
        double Li = 1;
        for (int j = 0; j < Count; j++) {
            if (i != j)
                Li *= (x - Array[j].x) / (Array[i].x - Array[j].x);
        }
        Polynom += Array[i].y * Li;
    }
    return Polynom;
}

//Полином Ньютона в точке
double PolynomN_dot(double x, Node* Array, int Count, double *DD_massiv)
{
    double Result = Array[0].y, dd;
    for (int i = 1; i < Count; i++) {
        dd = 0.0;
        dd = DD_massiv[i];
        for (int k = 0; k < i; k++)
            dd = dd * (x - Array[k].x);
        Result += dd; // Домножаем разделенную разность на скобки (x-x[0])...(x-x[i-1])
    }

    return Result;
}

void PolynomLG(Node* Array, int Count, double*& mas_coeff_polynom) //Заполняет массив коэффициентами интерполяционного полинома
{
    int i, j, p, c;
    double* tmp_massiv_coeff = new double[Count]; // Хранит коэффиценты Li
    double* brackets_massiv = new double[Count - 1]; //Для раскрытия (x-x0)...(x-xn), умноженных на коэффициент

    double denominator;

    for (j = 0; j < Count; j++) {
        mas_coeff_polynom[j] = 0.0;
    }

    for (i = 0; i < Count; i++) {
        c = 0; //Счетчик для числителя, т.к. иначе при счете по P будет дырка в массиве и get_koef_polynom не сможет посчитать (т.е. с делает сдвиг на  1 ячейку)
        for (j = 0; j < Count - 1; j++) {
            brackets_massiv[j] = 0.0;
            tmp_massiv_coeff[j] = 0.0;
        }
        tmp_massiv_coeff[Count - 1] = 0.0;
        denominator = Array[i].y;
        for (p = 0; p < Count; p++) {
            if (p != i) {
                denominator *= 1.0 / (Array[i].x - Array[p].x);
                brackets_massiv[c] = 0.0 - Array[p].x;
                c++;
            }
        }

        get_koef_polynom(brackets_massiv, tmp_massiv_coeff, Count - 1); //Раскрываем все скобки из возможных (по Лагранжу так) , -1 т.к. без i-ой
        for (j = 0; j < Count; j++) {
            mas_coeff_polynom[j] += (tmp_massiv_coeff[j] * denominator);
            mas_coeff_polynom[j] = (mas_coeff_polynom[j]);
        }
    }
    for (j = 0; j < Count; j++) {
        cout << mas_coeff_polynom[j] << " ";
    }
}

void PolynomN(Node* Array, int Count, double*& mas_coeff_polynom, double *DD_massiv)
{
    int j;
    double dd;
    double* tmp_massiv_coeff = new double[Count]; // массив А0...An конечных коэфф для одного шага цикла сколько точек, столько и коэфф
    double* brackets_massiv = new double[Count - 1]; //массив X0,,,Хn

    int i;

    for (int j = 0; j < Count - 1; j++) {
        mas_coeff_polynom[j] = 0.0;
        brackets_massiv[j] = 0.0;
    }

    mas_coeff_polynom[Count - 1] = 0.0;

    for (int j = 0; j < Count - 1; j++) { //Кладем все иксы
        brackets_massiv[j] = 0.0 - Array[j].x;
    }

    mas_coeff_polynom[0] = Array[0].y; // 0 индекс хранит свободный член итого многочлена

    for (i = 1; i < Count; i++) {
        for (j = 0; j < Count; j++) {
            tmp_massiv_coeff[j] = 0.0;
        }

        dd = DD_massiv[i]; //Ищем итую РР

        get_koef_polynom(brackets_massiv, tmp_massiv_coeff, i); //Возвращает кф при перемножении (х-Х0)...(х-Хi+1)

        for (j = 0; j < Count; j++) {
            tmp_massiv_coeff[j] = (tmp_massiv_coeff[j]);
            mas_coeff_polynom[j] += (tmp_massiv_coeff[j] * dd);
        }
    }
    for (j = 0; j < Count; j++) {
        cout << (mas_coeff_polynom[j]) << " ";
    }
}

void table_in_file(Node* MyArray, Node* Array, int MyCount, int Count, string file_name, string file_name_error, string method_name, double *DD_massiv)
{ //Строит таблицу для инт-го полинома и табл погрешности между ориг и ним
    int k;
    double y_value, x_value;

    ofstream fout(file_name); //Лежит таблица от интерполяции (по заданным точкам) по выбранному методу для графика (точки сами выбрали)
    ofstream pogr(file_name_error);

    method Func = PolynomN_dot;
    if (method_name == "Lagrange")
        method Func = PolynomLG_dot;

    for (k = 0; k < MyCount; k++) {
        x_value = MyArray[k].x;
        fout << x_value << " ";
        pogr << x_value << " ";
        y_value = Func(x_value, Array, Count, DD_massiv);

        fout << y_value << endl;

        pogr << abs(y_value - MyArray[k].y) << endl;
    }

    fout.close();
}

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

void ApproximateError(Node* Array, Node* MyArray, int MyCount, int Count) //Значение MaxDerivative получаем из питона и вручную подст
{ // Приближенная погрешность
    double MaxDerivative;
    ifstream N("D:/test1.txt");
    N >> MaxDerivative;
    N.close();


    //ifstream Max("D:/Max.txt");
    //Max >> MaxDerivative;
    //Max.close();

    ofstream AppError("D:/AppError.txt");

    double x_value, y_value;

    for (int k = 0; k < MyCount; k++) {
        x_value = MyArray[k].x;
        AppError << x_value << "   ";

        y_value = abs((MaxDerivative * W(x_value, Count, Array)) / Factorial(Count + 1));

        AppError << y_value << endl;
    }
    AppError.close();
}

void display_PolyDiff(double* mas_coeff_polynomLag, double* mas_coeff_polynomNew, int CountNodes)
{
    for (int i = 0; i < CountNodes; i++)
        cout << mas_coeff_polynomLag[i] - mas_coeff_polynomNew[i] << " ";
}

void Display_table_points(Node* Array, int count)
{
    for (int i = 0; i < count; i++) {
        cout << Array[i].x << "  " << Array[i].y << endl;
    }
}

int main()
{

    setlocale(LC_ALL, "RUS");
    //ДАННЫЕ
    Interval Interval;
    int CountNodes, MyNodes = 100;
    functiontype Func = &Myfunc;

    cout << "Enter the number of nodes: " << endl;
    cin >> CountNodes;
    CountNodes += 1;

    //Передача CountNodes В Питон
    ofstream N("D:/Degree_Derivative.txt");
    N << CountNodes + 1;
    N.close();
    system("pause");

    Interval.InitialNode = -2;
    Interval.EndNode = 2;

    double* mas_coeff_polynomNew = new double[CountNodes];
    double* mas_coeff_polynomLag = new double[CountNodes];



    //По заданным точкам(Nodes)
    Node* ArrayUniformNodes = new Node[CountNodes];                                                     // Массив с равномерными узлами
    ValueUniformTable(&Func, ArrayUniformNodes, Interval.InitialNode, Interval.EndNode, CountNodes);    // Заполнение таблицы равномерных значений

    //Построение графика(MyNodes)
    Node* Graph_ArrayUniformNodes = new Node[MyNodes];                                                  // Массив с равномерными узлами для построение графика (MyNodes)
    ValueUniformTable(&Func, Graph_ArrayUniformNodes, Interval.InitialNode, Interval.EndNode, MyNodes); //Заполнение таблицы для построение графика

    //Получаем массив РР для передачи функциям
    double* DD_massiv = new double[CountNodes];
    DD_massiv = DividedDifference(CountNodes, ArrayUniformNodes);

    //Чебышев
    //Node* ArrayChebyshevNodes = new Node[CountNodes];                                                    // Массив с Чебышевскими узлами
    //ValueChebyshevTable(&Func, ArrayChebyshevNodes, Interval.InitialNode, Interval.EndNode, CountNodes);   // Заполнение таблицы Чебышевских значений
    //table_in_file(Graph_ArrayUniformNodes, ArrayChebyshevNodes, MyNodes, CountNodes, "D:/Lagrange_Cheb.txt", "D:/Pogr_Lagrange_Cheb.txt", "Lagrange", DD_massiv);
    //table_in_file(Graph_ArrayUniformNodes, ArrayChebyshevNodes, MyNodes, CountNodes, "D:/Newton_Cheb.txt", "D:/Pogr_Newton_Cheb.txt", "Newton", DD_massiv);

    cout << endl;
    cout << "Таблица значений по равномерной сетке" << endl;
    Display_table_points(ArrayUniformNodes, CountNodes);
    cout << endl;

    //Заполнение массивов mas_coeff_polynomLag, mas_coeff_polynomNew коэффициентами при соот-х степенях полиномов
    cout << "Полином Лагранжа:" << endl;
    PolynomLG(ArrayUniformNodes, CountNodes, mas_coeff_polynomLag);
    cout << endl << endl << "Полином Ньютона:" << endl;
    PolynomN(ArrayUniformNodes, CountNodes, mas_coeff_polynomNew, DD_massiv);
    cout << endl << endl;

    //Вывод разности между полиномами Лагранжа и Ньютона
    cout << "Разность между полиномами Лагранжа и Ньютона:" << endl;
    display_PolyDiff(mas_coeff_polynomLag, mas_coeff_polynomNew, CountNodes);
    cout << endl;

    //Запись точек в файл для построения
        //Равномерная сетка
    table_in_file(Graph_ArrayUniformNodes, ArrayUniformNodes, MyNodes, CountNodes, "D:/Lagrange.txt", "D:/Pogr_Lagrange.txt", "Lagrange", DD_massiv);
    table_in_file(Graph_ArrayUniformNodes, ArrayUniformNodes, MyNodes, CountNodes, "D:/Newton.txt", "D:/Pogr_Newton.txt", "Newton", DD_massiv);
    orig_table_in_file(Graph_ArrayUniformNodes, MyNodes); //Подставляет 100 точек в ориг функцию

    // Приближенная погрешность
    ApproximateError(ArrayUniformNodes, Graph_ArrayUniformNodes, MyNodes, CountNodes);

    system("pause");
    return 0;
}