#include <iostream>
#include <math.h>
#include <cmath>
#include <vector>
#include <Windows.h>
#include <stdlib.h>
#include <fstream>
#include <string.h>
#include <time.h>
#include <iomanip>
using namespace std;

const double PI = 3.141592653589793238463;

//Для хранения сетки
typedef struct Node
{
    double x, y;
} Node;

//Задание функции
typedef double(*functiontype)(double x);
double Myfunc(double x)
{
    return x*cos(x);
    //return x * x;
    //return sin(x);
}
functiontype Func = &Myfunc;

//Приблизительное значение 2 производной
double approximate_derivative_2(functiontype* f, double x) {
    double delta = 1e-5;
    return (((*f)(x) - 2 * (*f)(x + delta) + (*f)(x + 2 * delta)) / (delta*delta));
}

//Приблизительное значение 2 производной
double approximate_derivative_1(functiontype* f, double x) {
    double delta = 1e-5;
    return (((*f)(x + delta) - (*f)(x)) / delta);
}

//Возвращает значение сплайна в точке
double Spline(double x, int i, double *Array_steps, double *x_massiv, Node* Array)
{
    double h2, g1, g2, x1, x2, y1, y2;

    h2 = Array_steps[i + 1]; //этот массив с  индексацией от 1, размер - CountSegments, т.е. кол-в точек - 1
    g1 = x_massiv[i];       g2 = x_massiv[i + 1]; //тут обращение в конце к n+1, а в степс к n при том же индексе, всё норм
    x1 = Array[i].x;        x2 = Array[i + 1].x;
    y1 = Array[i].y;        y2 = Array[i + 1].y;

    return (
        (x*x*x)*((g2 - g1) / (6 * h2)) +
        (x*x)*((g1*x2 - g2 * x1) / (2 * h2)) +
        (x)*((g2*x1*x1 - g1 * x2*x2) / (2 * h2) + (g1*h2 - g2 * h2) / 6 + (y2 - y1) / h2) +
        ((g1*x2*x2*x2 - g2 * x1*x1*x1) / (6 * h2) + (g2*h2*x1 - g1 * h2*x2) / 6 + (y1*x2 - y2 * x1) / h2)
        );

}

//Равномерная сетка
void ValueUniformTable(functiontype* f, Node* Array, double Initial, double End, int CountSegments)
{
    int i;
    double h;

    Array[0].x = Initial; Array[0].y = (*f)(Initial);
    h = abs(End - Initial) / CountSegments;

    cout << "(" << Array[0].x << ":" << Array[0].y << ")" << endl;
    for (int i = 1; i <= CountSegments; i++)
    {
        Array[i].x = Array[i - 1].x + h;
        Array[i].y = (*f)(Array[i].x);
        cout << "(" << Array[i].x << ":" << Array[i].y << ")" << endl;
    }

}


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


//Неравномерная сетка
void ValueIrregularTable(functiontype* f, Node* Array, double Initial, double End, int CountSegments, double *Array_steps)
{
    int i;
    double alpha, h;

    Array[0].x = Initial;           Array[CountSegments].x = End;
    Array[0].y = (*f)(Array[0].x);  Array[CountSegments].y = (*f)(Array[CountSegments].x);
    h = (Array[CountSegments].x - Array[0].x) / CountSegments;
    Array_steps[0] = 0;
    alpha = 2 * PI / CountSegments;

    for (i = 1; i < CountSegments; i++) {
        Array_steps[i] = h + (9 * h / 10)*cos(alpha*i);     //hi = h_i-h_i-1, то есть шаг назад
        Array[i].x = Array[i - 1].x + Array_steps[i];
        Array[i].y = (*f)(Array[i].x);
    }
    Array_steps[CountSegments] = Array[CountSegments].x - Array[CountSegments - 1].x;

    //вывод в консоль
    for (i = 0; i <= CountSegments; i++)
        cout << "(" << Array[i].x << ":" << Array[i].y << ")" << endl;

}

/* Принимает на вход коэффициенты СЛАУ в виде массива matrix_coeffs - массив вида [[a1,b1,c1,d1],[a2,b2,c2,d2],...,[an,bn,cn,dn].
   Заполняет Двумерный массив прогоночных коэффициентов coeffs_massiv (0 - Кси, 1 - Эта).
   Находит решение в виде массива x_massiv
 */
void tridiagonal_matrix_algorithm(double** matrix_coeffs, double* Array_steps, Node* Array, int CountSegments, double* x_massiv, double B) {
    int i;
    double denominator;
    double K, E, K_prev = 0, E_prev = 0; //Прогоночные коэффициенты
    double y1, y2, h2;
    int matrix_size = CountSegments - 1;

    double** coeffs_massiv;
    coeffs_massiv = new double*[matrix_size + 1];
    for (i = 0; i <= matrix_size ; i++)
        coeffs_massiv[i] = new double[2];

    coeffs_massiv[0][0] = 0; coeffs_massiv[0][1] = 0; //K1,E1 = 0
//  matrix_coeffs[0][0] = 0; matrix_coeffs[matrix_size - 1][2] = 0; //а1=сn=0 - убрать если работает без этого, написано в др функции

    //Прямой ход
    for (i = 1; i <= matrix_size; i++) { //Ищем K_i+1 и E_i+1 на iом шаге
        denominator = (matrix_coeffs[i - 1][0] * K_prev + matrix_coeffs[i - 1][1]);
        K = -(matrix_coeffs[i - 1][2]) / denominator;
        E = (matrix_coeffs[i - 1][3] - matrix_coeffs[i - 1][0] * E_prev) / denominator;
        K_prev = K;
        E_prev = E;
        coeffs_massiv[i][0] = K; coeffs_massiv[i][1] = E; //лежат со сдвигом из за индексации
    }

    //Обратный ход
    //x_massiv[matrix_size-1] = coeffs_massiv[matrix_size-1][1]; //это gn-1. [g0,...gn] - n+1 значение, перепишу в соот-вии с этим
    x_massiv[matrix_size+1] = 0; //Зануляю для кореткной работы алгоримта, считая что это xn+1 в терминах прогонки

    for (i = matrix_size; i >= 1; i--) {
        x_massiv[i] = coeffs_massiv[i][0] * x_massiv[i + 1] + coeffs_massiv[i][1];
    }

    //Определяю gn
    h2 = Array_steps[CountSegments]; //вроде так, h0 условно существует и h1 лежит на 1 инд, h_n-1 na matrix size
    y1 = Array[CountSegments - 1].y;
    y2 = Array[CountSegments].y;
    x_massiv[CountSegments] = 6 * (y2 - y1 - B * h2) / (h2*h2) - 2 * x_massiv[CountSegments - 1]; //gn - записываю чему действительно является в моей задачи (посл ур)

    cout << endl << "Размер матрицы = " << matrix_size << ", Гамм" << matrix_size+2 << endl;
    for (i = 0; i < matrix_size+2; i++) //+1 потому что matrix_size = CS-1 тут мб??
        cout << "g_ " << i << ": " << x_massiv[i] << endl;

    cout << endl << endl;
    for (i = 0; i < CountSegments; i++) //+1 потому что matrix_size = CS-1 тут мб??
        cout << "e,k_i " << i << ": " << coeffs_massiv[i][0] << " " << coeffs_massiv[i][1] << endl;


}

//Составляем СЛАУ с учётом нач данных и закидываем это всё в один массив
void get_matrix_coeffs(double Initial, double End, int CountSegments, functiontype* f, double** matrix_coeffs, Node* Array, double *Array_steps, double B) {
    int i;
    double h1, h2, h_prev, d2, g2;
    double x1, x2, y1, y2, y3, denominator;
    int matrix_size = CountSegments - 1;

    for (i = 1; i <= matrix_size; i++) { //в системе i=1,..n-1, g0,gn - известны
        h1 = Array_steps[i]; //i
        h2 = Array_steps[i + 1];; //i+1
        h_prev = h2;
        matrix_coeffs[i - 1][0] = h1;
        matrix_coeffs[i - 1][1] = 2 * (h1 + h2);
        matrix_coeffs[i - 1][2] = h2;
        matrix_coeffs[i - 1][3] = 6 * ((Array[i + 1].y - Array[i].y) / h2 - (Array[i].y - Array[i - 1].y) / h1);
    }

    /*//Переопределние последнего уравнения для моего случая:
    h1 = Array_steps[matrix_size];
    y1 = Array[matrix_size - 1].y;
    y2 = Array[matrix_size].y;

    matrix_coeffs[matrix_size - 1][0] = h1;
    matrix_coeffs[matrix_size - 1][1] = 2 * h1;
    matrix_coeffs[matrix_size - 1][2] = 0; //cn=0
    matrix_coeffs[matrix_size - 1][3] = 6 * (B - (y2 - y1) / h1);
    matrix_coeffs[0][0] = 0; //а1=0*/

    matrix_coeffs[0][0] = 0; //а1=0
    matrix_coeffs[matrix_size-1][2] = 0; //cn=0

    cout  << "(" << matrix_coeffs[0][1] << ")" << "  " << "(" << matrix_coeffs[0][2] << ")" << "  " <<  "(" << "0" << ")" << endl;
    for (i = 1; i < matrix_size-1; i++) {
        cout << "(" << matrix_coeffs[i][0] << ")" << "  " << "(" << matrix_coeffs[i][1] << ")" << "  " << "(" << matrix_coeffs[i][2] << ")" << "  " << endl;
    }
    cout << "(" << "0" << ")" << "  " << "(" << matrix_coeffs[matrix_size-1][0] << ")" << "  " << "(" << matrix_coeffs[matrix_size-1][1] << ")" << "  " << endl;
}

/*
Строит график функции на всём интервале с заданным количеством точек Count_dots
*/
void orig_table_in_file(functiontype* f, int Count_dots, double Initial, double End) {

    int k, i;
    double x_value, step;

    step = (End - Initial) / (Count_dots - 1);
    x_value = Initial - step;
    ofstream fout("D:/original_graphic.txt");

    for (i = 0; i < Count_dots; i++) {
        x_value += step;
        fout << x_value << " ";
        fout << (*f)(x_value) << endl;
    }
    fout.close();
}

//Используя функцию Spline записываем в файл таблицу значений спалйна на всём интервале в Count_dots точках
void spline_table_in_file(Node* Array, int Count_dots, int Count_Segments, double Initial, double End, double *Array_steps, double *x_massiv) { // Функция берет таблицу иксов и игриков, количество точек в которых считаем знач полинома (мб убрать) и коэфф полинома,

    int i = 0;
    double step = (End - Initial) / (Count_dots - 1), y_value, x_value;

    ofstream fout("D:/spline_graphic.txt");
    x_value = Initial;
    for (i = 0; i < Count_Segments; i++) {
        while (x_value < Array[i + 1].x) {
            fout << x_value << " ";
            y_value = Spline(x_value, i, Array_steps, x_massiv, Array);
            fout << y_value << endl;
            x_value += step;
        }
    }
    fout << End << " ";
    fout << Spline(End, Count_Segments - 1, Array_steps, x_massiv, Array) << endl;

    fout.close();
}

int main()
{
    setlocale(LC_ALL, "RUS");
    functiontype Func = &Myfunc;
    double Initial = -5, End = 5, B, A;
    int CountSegments, i, Countdots = 5000;
    cout << "Введите N: ";
    //cin >> CountSegments;
    CountSegments = 10;
    cout << "Точек: " << CountSegments + 1 << endl << endl;
    cout << "Введите A,B: " << endl << endl;
    //  cin >> A >> B;
    A = approximate_derivative_2(&Func, Initial); B = approximate_derivative_1(&Func, End);

    double *Array_steps = new double[CountSegments+1];
    double *Array_steps_cheb = new double[CountSegments + 1];
    //Построение сетки
    //Равномерная сетка
    cout << "Равномерная Сетка: " << endl;
    Node* ArrayUniformNodes = new Node[CountSegments + 1];
    ValueUniformTable(&Func, ArrayUniformNodes, Initial, End, CountSegments);

    //Чебышёв
    Node* ArrayChebyshevNodes = new Node[CountSegments + 1]; // Массив с Чебышевскими узлам
    ValueChebyshevTable(&Func, ArrayChebyshevNodes, Initial, End, CountSegments+1, Array_steps_cheb); // Заполнение таблицы Чебышевских значений

    //Неравномерная сетка
    cout << endl << "Неравномерная Сетка: " << endl;
    Node* ArrayIrregularNodes = new Node[CountSegments + 1];
    ValueIrregularTable(&Func, ArrayIrregularNodes, Initial, End, CountSegments, Array_steps);
    cout << endl;

    //Заполнение массива коэффициентами матрицы
    double** matrix_coeffs;
    matrix_coeffs = new double*[CountSegments - 1];
    for (i = 0; i < CountSegments - 1; i++)
        matrix_coeffs[i] = new double[3];


    for (i = 0; i <= CountSegments; i++) {
        cout << "h_" << i << ": " << Array_steps_cheb[i] << endl;
    }

    cout << endl << "Полученная матрица: " << endl;
    get_matrix_coeffs(Initial, End, CountSegments, &Func, matrix_coeffs, ArrayChebyshevNodes, Array_steps_cheb, B);//  ДЛЯ ЧЕБЫШЕВА ПРОВЕРИТЬ ВСЕ МАССИВЫ ЧТОБЫ НОРМ ПО НЕМУ ЗАДАВАЛИСЬ (ЭРЭЙ СТЕПС)

    cout << "Правая часть: " << endl << endl;
    for (i = 0; i < CountSegments-1; i++) {
        cout << "d_" << i << ": " << matrix_coeffs[i][3] << endl;
    }


    //Получение ответа в x_massiv
    double* x_massiv = new double[CountSegments + 1];
    x_massiv[0] = A; //g0
    tridiagonal_matrix_algorithm(matrix_coeffs, Array_steps_cheb, ArrayChebyshevNodes, CountSegments, x_massiv, B); //ПОч CountSegments - 1 было? счиьаю ошибкой

    //Вывод интерполируемой функции
    orig_table_in_file(&Func, Countdots, Initial, End);

    //Вывод сплайна на экран
    spline_table_in_file(ArrayChebyshevNodes, Countdots, CountSegments, Initial, End, Array_steps_cheb, x_massiv);

    //Приближение производных на кроях
    cout << endl;
    cout << "1 производная на левом конце: " << approximate_derivative_1(&Func, Initial) << endl;
    cout << "2 производная на левом конце: " << approximate_derivative_2(&Func, Initial) << endl;
    cout << "1 производная на правом конце: " << approximate_derivative_1(&Func, End) << endl;
    cout << "2 производная на правом конце: " << approximate_derivative_2(&Func, End) << endl;


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