/*
plot "D:\graphic.txt" with lines title "График U(x)"
plot "D:\graphic.txt" with lines title "График U(x)", "D:\graphic.txt" with points title "" pointtype 3 pointsize 1 lc 200 200 200
*/
#include <iostream>
#include <math.h>
#include <cmath>
#include <stdlib.h>
#include <fstream>
using namespace std;
const double PI = 3.141592653589793238463;

//Задание функции
typedef double(*functiontype1)(double x);
typedef double(*functiontype2)(double x, double s);
double f(double x)
{
    return 28*x * x * x + x - 56.33333333;
}
double K(double x, double s)
{
    return x * x*x -s;
}
functiontype1 Func1 = &f;
functiontype2 Func2 = &K;

//Задает набор узлов, равномерно удаленных друг от друга
void ValueUniformTable(double* Array, double Initial, double End, int CountDots)
{
    double h = abs(End - Initial) / (CountDots - 1);
    Array[0] = Initial;
    for (int i = 1; i < CountDots; i++)
        Array[i] = Array[i - 1] + h;
}

void Gaussian_elimination(double** coeffs_massiv, double* y_Massiv, double* U_massiv, int n) {

    double* tmp_mas;
    tmp_mas = new double[n];

    int i, j, k, index;
    double tmp, max, r;

    for (j = 0; j <= n - 2; j++) {

        max = abs(coeffs_massiv[j][j]);
        index = j;

        for (i = index + 1; i <= n - 1; i++) {
            if (abs(coeffs_massiv[i][index]) > max) {
                max = abs(coeffs_massiv[i][index]);
                index = i;
            }
        }

        if (index != j) {
            for (i = 0; i <= n - 1; i++) {
                tmp_mas[i] = coeffs_massiv[j][i];
                coeffs_massiv[j][i] = coeffs_massiv[index][i];
                coeffs_massiv[index][i] = tmp_mas[i];
            }
            tmp = y_Massiv[j];
            y_Massiv[j] = y_Massiv[index];
            y_Massiv[index] = tmp;
        }

        for (i = j + 1; i <= n - 1; i++) {
            r = -coeffs_massiv[i][j] / coeffs_massiv[j][j];

            for (k = j; k <= n - 1; k++) {
                coeffs_massiv[i][k] = coeffs_massiv[i][k] + r * coeffs_massiv[j][k];
            }
            y_Massiv[i] = y_Massiv[i] + r * y_Massiv[j];
        }
    }


    U_massiv[n - 1] = y_Massiv[n - 1] / coeffs_massiv[n - 1][n - 1];
    for (i = n - 2; i >= 0; i--) {
        tmp = 0.0;
        for (int j = i + 1; j <= n - 1; j++)
            tmp += coeffs_massiv[i][j] * U_massiv[j];
        U_massiv[i] = (y_Massiv[i] - tmp) / coeffs_massiv[i][i];
    }
}


//Запись в файл (x_i, U_i)
void table_in_file(double* x_Array, int CountDots, double* U_massiv)
{
    ofstream fout("D:/graphic.txt");
    for (int i = 0; i < CountDots; i++)
        fout << x_Array[i] << " " << U_massiv[i] << endl;

    fout.close();
}

//Массив, хранящий коэффиценты системы
void get_coeffs_massiv(functiontype2* K, functiontype1* f, int CountDots, double** coeffs_massiv, double* x_Array, double* y_Array) {
    int i, j;
    for (i = 0; i < CountDots; i++) {
        for (j = 0; j < CountDots; j++) {
                coeffs_massiv[i][j] = (x_Array[j] - x_Array[j - 1]) * (*K)(x_Array[i], x_Array[j]);
                if  (j == 0) coeffs_massiv[i][j] = 0;
                if (i == j) coeffs_massiv[i][j] += 1;

        }
        y_Array[i] = (*f)(x_Array[i]);
    }
}

int main()
{
    setlocale(LC_ALL, "RUS");
    int i, j;

    /*---------------------------------Входные данные-------------------------------------*/
    double Initial = 0, End = 4 ;
    int CountDots = 50;


    /*---------------------------------Задаём необходимые массивы-------------------------*/
    double** coeffs_massiv; //массив хранящий коэффиценты СЛАУ
    coeffs_massiv = new double* [CountDots];
    for (i = 0; i < CountDots; i++)
        coeffs_massiv[i] = new double[CountDots];

    double* U_massiv = new double[CountDots]; //Хранит U_i
    double* y_Array = new double[CountDots]; //Правая часть СЛАУ
    double* x_Array = new double[CountDots]; //Узлы


    /*---------------------------------Заполняем и решаем---------------------------------*/
    ValueUniformTable(x_Array, Initial, End, CountDots); //Строим сетку
    get_coeffs_massiv(&Func2, &Func1, CountDots, coeffs_massiv, x_Array, y_Array); //Строим СЛАУ
    Gaussian_elimination(coeffs_massiv, y_Array, U_massiv, CountDots);
    table_in_file(x_Array, CountDots, U_massiv);


    /*---------------------------------Вывод в консоль------------------------------------*/

    /*cout << endl << "Полученная СЛАУ: " << endl;
    for (i = 0; i < CountDots; i++) {
        for (int j = 0; j < CountDots; j++) {
            cout << coeffs_massiv[i][j] << " ";
        }
        cout << y_Array[i] << endl;
        cout << endl;
    }*/

    cout << endl << "Узлы: " << endl;
    for (j = 0; j < CountDots; j++) {
        cout << x_Array[j] << " ";
    }

    cout << endl << endl << "U_i: " << endl;
    for (j = 0; j < CountDots; j++) {
        cout << U_massiv[j] << " ";
    }

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