#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 struct Interval
{
    double InitialNode, EndNode;
} Interval;

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

double Myfunk(double x)
{
    return (x*x);
}

double exponenta(double x)
{ // Экспонента
    return exp(x);
}

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);
    //std::cout << Array[0].x << "      " << Array[0].y << endl;
    for (int i = 1; i < Count; i++)
    {
        Array[i].x = Array[i - 1].x + step;
        Array[i].y = (*f)(Array[i].x);
        //std::cout << Array[i].x << "      " << Array[i].y << endl;
    }
}

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

void get_koef_polynom(double* mas, double* res, int n, int k)
{ // Функция на вход берет массив из n-1 точек и
  // возвращает кф полином получившегося при
  // раскрытии скобок !!!(И еще надо домножить на
  // игрик)!!!
    int i;
    double* anx = new double[n + 5];
    double* ann = new double[n + 5];
    res[0] = mas[0];
    res[1] = 1;
    int j;
    for (j = 0; j <= n; j++)
    {
        anx[j] = 0.0;
    }

    for (j = 0; j < n; j++)
    {
        ann[j] = 0.0;
    }
    for (i = 1; i < k; 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];
        }
}

//Возвращает число, возведенное в 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 round(double x) {
    if (x < 0.0000000001)
        x = 0.0;
    return x;
}

// Полином Лагранжа в точке
double PolynomLG_dot(double x, Node *Array, int Count)
{
    double Polynom = 0;
    for (int i = 0; i < Count; i++)
    { // Счетчик по функциям l_i
        double Li = 1;
        for (int j = 0; j < Count; j++)
        { // Счетчик по x_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 Result = Array[0].y;
    for (int i = 1; i < Count; i++)
    {
        double Polynom = 0;
        for (int j = 0; j <= i; j++)
        { // Следующее слагаемое полинома
            double Composition = 1;
            for (int k = 0; k <= i; k++)
            { // Считаем знаменатель разделенной разности
                if (k != j) Composition *= (Array[j].x - Array[k].x);
            }
            Polynom += Array[j].y / Composition; // Считаем разделенную разность
        }
        for (int k = 0; k < i; k++) Polynom *= (x - Array[k].x); // Домножаем разделенную разность на скобки (x-x[0])...(x-x[i-1])
        Result += Polynom; // Полином
    }
    return Result;
}

//Подходит для построенных полиномов, считает их знач в данных нам n узлах, но если их слишком мало? Надо добавить!
void table_in_file(Node* Array, int n, int z, double* polynom_koeff, string file_name, string file_name_pogr) { // Функция берет таблицу иксов и игриков, количество точек в которых считаем знач полинома (мб убрать) и коэфф полинома,
                                                                // По итоге имеем файл с 2 столбацами: х и у, по которому гну пло
    int k, j;
    double y_value, x_value;

    ofstream fout(file_name);
    ofstream pogr(file_name_pogr);


    for (k = 0; k < n; k++) { //n иксов подставляев в полином степени n-1
        x_value = Array[k].x;
        fout << x_value << " ";
        pogr << x_value << " ";
        y_value = 0.0;
        for (j = 0; j < z; j++) { //убрал n+1 потому что вылезал за границы памяти но поч так хз
            //cout << "x: " << x_value << " ";
            //cout << polynom_koeff[j] << " ";
            //cout << "x: ^  " << get_degree(x_value, j) << " ";  //Чтобы не умножать ахуенно маленькое число на икс в большой степени надо знать какой степени у нас многочлен , а то e^-17 умножаем на x^10
            y_value += polynom_koeff[j] * get_degree(x_value, j); //На 0 индексе стоит свободные член, умножаем на икс в 0 стпени = 1
            //cout << y_value << endl << endl << endl;

        }
        fout << y_value << endl;
        //cout << "y_value: " << y_value << endl;
        //cout << "Array[k].y: " << Array[k].y << endl;
        //cout << "Array[k].y - y_value:  " << Array[k].y - y_value << endl;

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

    }

    fout.close();
}

//То же самое но значение считает сразу вточке
void table_in_file2(Node* Array, int n, int z, double* polynom_koeff, string file_name, string file_name_pogr, string method_name) { // Функция берет таблицу иксов и игриков, количество точек в которых считаем знач полинома (мб убрать) и коэфф полинома,
                                                                // По итоге имеем файл с 2 столбацами: х и у, по которому гну пло
    int k, j;
    double y_value, x_value;

    ofstream fout(file_name);
    ofstream pogr(file_name_pogr);




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

        if (method_name == "Lagrange")
            y_value = PolynomLG_dot(x_value, Array, n);
        else
            y_value = PolynomN_dot(x_value, Array, n);

        fout << y_value << endl;
        //cout << "y_value: " << y_value << endl;
        //cout << "Array[k].y: " << Array[k].y << endl;
        //cout << "Array[k].y - y_value:  " << Array[k].y - y_value << endl;

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

    }

    fout.close();
}

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

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

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

    fout.close();
}



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

void AbsoluteError(functiontype *Function, double Initial, double End, Node *Array, int Count)
{ // Абсолютная погрешность

    for (double x = Initial; x <= End; x += 0.01f)
    {
        abs((*Function)(x) - PolynomLG_dot(x, Array, Count));
    }
}

void ApproximateError(functiontype *Function, double Initial, double End, Node *Array, int Count)
{ // Приближенная погрешность

    for (double x = Initial; x <= End; x += 0.01f)
    {
        abs(exp(End) * W(x, Count, Array) / Factorial(Count + 1));
        abs(exp(End) * W(x, Count, Array) / Factorial(Count + 1));
    }
}

void PolynomLG(Node* Array, int n, double* &konmas)
{
    int i, j, l, p, z, c;
    double* massiv2 = new double[n]; // массив А0...An конечных коэфф для одного
                                     // шага цикла сколько точек, столько и коэфф
    double* massiv = new double[n - 1]; //массив X0,,без Хi,Хn домноженных на коэфф.  ПОЭТОМУ Н-1
    //double* konmas = new double[n]; // будет коннечный массив коэфф
    double k;
    for (j = 0; j < n; j++)
    {
        konmas[j] = 0.0;
    }
    for (i = 0; i < n; i++)
    {
        c = 0;
        for (j = 0; j < n - 1; j++)
        {
            massiv[j] = 0.0;
            massiv2[j] = 0.0;
        }
        massiv2[n - 1] = 0.0;
        k = Array[i].y;
        for (p = 0; p < n; p++)
        {
            if (p != i)
            {
                k *= 1.0 / (Array[i].x - Array[p].x);
                massiv[c] = 0.0 - Array[p].x;
                c++;
            }
        }
        cout << endl;
        get_koef_polynom(massiv, massiv2, n - 1, n - 1);
        for (l = 0; l < n; l++)
        {
            konmas[l] += massiv2[l] * k;
        }
        /*cout << "konmas  ";
        for (z = 0; z < n; z++)
        {
            cout << konmas[z] << ' ';
        }*/
    }

    //table_in_file(Array, n, konmas); //
    cout << endl;
}
double DividedDifference(int i, Node* Array)
{ // Разделенная разность
    double DD = 0;
    for (int j = 0; j <= i; j++)
    {
        double tmp = 1;
        for (int k = 0; k <= i; k++)
        {
            if (k != j)
                tmp *= Array[j].x - Array[k].x;
        }
        DD += Array[j].y / tmp;
    }

    return DD;
}

void PolynomN(Node* Array, int n, double* &konmas)
{
    int j;
    double dd;
    double* massiv2 = new double[n]; // массив А0...An конечных коэфф для одного
                                     // шага цикла сколько точек, столько и коэфф
    double* massiv = new double[n - 1]; //массив X0,,,Хn                 +1 для зануления
    //double* konmas = new double[n]; // будет коннечный массив коэфф
    int i;

    for (int j = 0; j < n - 1; j++)
    {
        konmas[j] = 0.0;
        massiv[j] = 0.0;
    }

    konmas[n - 1] = 0.0;

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

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

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

        dd = DividedDifference(i, Array); //Ищем итую РР
    //  cout << "DD " << dd << endl << endl;
        get_koef_polynom(massiv, massiv2, n - 1,
            i); //Возвращает кф при перемножении (х-Х0)...(х-Хi+1),
        for (j = 0; j < n; j++)
        {
            konmas[j] += massiv2[j] * dd;

        }

        cout << endl;

        cout << endl;
    }

    /*for (j = 0; j < n; j++)
    {
        konmas[j] = round(konmas[j]);
        cout << konmas[j] << "   ";
    }*/
    cout << endl;


}

double DFunc(functiontype* func, double x,
    int n) //возварщает проивзодную нго порядка как константу типа дабл
{ // Производная функции
    double h = 0.00001;
    if (n == 1)
    {
        return ((*func)(x + h) - (*func)(x)) / h;
    }
    else
    {
        return (DFunc(func, x + h, n - 1) - DFunc(func, x, n - 1)) / h;
    }
}

double test(functiontype* func, double x, int n, long double h)
{ // Производная функции
    h = 0.0001;
    while (n > 1)
    {
        return (test(func, x + h, n - 1, h / pow(10, -2))
            - test(func, x - h, n - 1, h / pow(10, -2)))
            / 2 * h;
    }
    if (n == 1)
    {
        return ((*func)(x + h) - (*func)(x - h)) / 2 * h;
    }
}



int main()
{

    setlocale(LC_ALL,"RUS");

    //ПОСТРОЕНИЕ ДАННЫХ
    Interval Interval;
    int CountNodes, MyNodes = 100;
    functiontype Func = &Myfunk;

    /*cout << "Enter the interval: " << endl;
    cin >> Interval.InitialNode >> Interval.EndNode;
    cout << "Enter the number of nodes: " << endl;
    cin >> CountNodes;*/

    //Ввод вручную
    Interval.InitialNode = -2;
    Interval.EndNode = 2;
    CountNodes = 25;

    double* konmas = new double[CountNodes];

    Node *ArrayUniformNodes = new Node[CountNodes];// Массив с равномерными // узлами
    Node *ArrayUniformNodes2 = new Node[MyNodes];
    //Node* ArrayChebyshevNodes = new Node[CountNodes]; // Массив с Чебышевскими узлами
    Node *ArrayformNodes2 = new Node[100]; //Массив где будет много точек для построение ориг графика

    ValueUniformTable(&Func, ArrayUniformNodes, Interval.InitialNode, Interval.EndNode, CountNodes); // Заполнение таблицы равномерных значений
    ValueUniformTable(&Func, ArrayUniformNodes2, Interval.InitialNode, Interval.EndNode, MyNodes);
    //ValueChebyshevTable(&Func, ArrayChebyshevNodes, Interval.InitialNode, Interval.EndNode, CountNodes); // Заполнение таблицы Чебышевских значений


    // ЗАПУСКИ ФУНКЦИЙ
    //cout << "Лагранж в точке:" << PolynomLG_dot(5.00, ArrayUniformNodes, CountNodes);


    cout << endl;
    //PolynomLG(ArrayUniformNodes, CountNodes, konmas);
    cout << endl;
    //PolynomN(ArrayUniformNodes, CountNodes,konmas);
    table_in_file2(ArrayUniformNodes2, MyNodes, CountNodes, konmas, "D:/Lagrange.txt", "D:/Pogr_Lagrange.txt", "Lagrange");
    table_in_file2(ArrayUniformNodes2, MyNodes, CountNodes, konmas, "D:/Newton.txt", "D:/Pogr_Newton.txt", "Newton");
    table_in_file(ArrayUniformNodes2, MyNodes, CountNodes, konmas, "D:/withkonmas.txt", "D:/Pogr_withkonmas.txt");
    orig_table_in_file(ArrayUniformNodes2, MyNodes);


    // Абсолютная погрешность
    //AbsoluteError(&Func, Interval.InitialNode, Interval.EndNode, ArrayUniformNodes, CountNodes);
    // Приближенная погрешность
    //ApproximateError(&Func, Interval.InitialNode, Interval.EndNode, ArrayUniformNodes, CountNodes);

    //system("pause");
    return 0;
}