/*
plot "D:\U(x).txt" with lines title "График U(x)", "D:\dots.txt" with points title "" pointtype 3 pointsize 1 lc 200 200 200

plot "D:\U(x).txt" with lines title "График U(x)", "D:\dots_4.30.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;
/*
K(x,s) = x^2+s^2    , U(x) = x^3    ;  [0,5], f(x) = x^3 + 624/4*x^2 + 15625/6
K(x,s) = x^3-s      , U(x) = x+5    ; [0,4] ; f(x) = 28*x^3+x-169/3:
K(x,s) = sin(x+s)   , U(x)=x^2      ; [0,5] ; f(x) = x^2+10sin(x+5)-2cos(x)-23cos(x+5):
K(x,s) = x + s      , U(x) = 1 - |x|; [-1,1]; f(x) = 1 + x - |x|:
*/
//Задание функции
typedef double(*functiontype1)(double x);
typedef double(*functiontype2)(double x, double s);
double K(double x, double s)
{
    //return x * x + s * s;
    //return x*x*x - s;
    //return sin(x+s);
    return x + s;
}
double U(double x)
{
    //return x * x * x;
    //return x + 5;
    //return x*x;
    return 1 - abs(x);
}
double f(double x)
{
    //return x * x * x + 625  * x * x /4+ 15625 / 6;
    //return 28*x*x*x + x - 169/3;
    //return x*x + 10*sin(x+5) - 2*cos(x) - 23*cos(x+5);
    return 1 + x - abs(x);
}
functiontype1 Func1 = &f;
functiontype2 Func2 = &K;
functiontype1 Func3 = &U;

//Задает набор узлов, равномерно удаленных друг от друга
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];
    }
}

void table_in_file(functiontype1* U, double* x_Array, int CountDots, double* U_massiv, double* U_massiv2)
{
    int i;
    ofstream fout1("D:/dots.txt");
    ofstream fout3("D:/dots_4.30.txt");
    ofstream fout2("D:/U(x).txt");

    //Запись точек (x_i,U_i)
    for (i = 0; i < CountDots; i++) {
        fout1 << x_Array[i] << " " << U_massiv[i] << endl;
        fout3 << x_Array[i] << " " << U_massiv2[i] << endl;
    }

    //Запись настоящей функции U(x)
    double h = abs(x_Array[CountDots-1] - x_Array[0]) / (300 - 1);
        for (i = 0; i < 300; i++)
            fout2 << x_Array[0] + h* i << " " <<(*U)(x_Array[0] + h * i) << endl;

    fout1.close();
    fout2.close();
    fout3.close();
}


//Массив, хранящий коэффиценты системы
void get_coeffs_massiv(functiontype2* K, functiontype1* f, int CountDots, double** coeffs_massiv, double** coeffs_massiv2,
    double* x_Array, double* y_Array, double* y_Array2, double  h) {

    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;

            //4.30
            coeffs_massiv2[i][j] = h * (*K)(x_Array[i], x_Array[j]);
            if (j == 0 || j == CountDots - 1)
                coeffs_massiv2[i][j] *= 0.41666666;

            if (j == 2 || j == CountDots - 2)
                coeffs_massiv2[i][j] *= 1.08333333;

            if (i == j)
                coeffs_massiv2[i][j] += 1;

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

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

    /*---------------------------------Входные данные-------------------------------------*/
    //double Initial = 0, End = 5;
    //double Initial = 0, End = 4 ;
    //double Initial = 0, End = 5 ;
    double Initial = -1, End = 1;

    int CountDots = 5;
    double h = abs(End - Initial) / (CountDots - 1);

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

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


    /*---------------------------------Заполняем и решаем---------------------------------*/
    ValueUniformTable(x_Array, Initial, End, CountDots); //Строим сетку
    get_coeffs_massiv(&Func2, &Func1, CountDots, coeffs_massiv, coeffs_massiv2, x_Array, y_Array, y_Array2, h); //Строим СЛАУ
    cout << endl << "Полученная СЛАУ: " << endl;
    for (i = 0; i < CountDots; i++) {
        for (int j = 0; j < CountDots; j++) {
            cout << coeffs_massiv2[i][j] << " ";
        }
        cout << y_Array[i] << endl;
        cout << endl;
    }
    Gaussian_elimination(coeffs_massiv, y_Array, U_massiv, CountDots); //Обыч система
    Gaussian_elimination(coeffs_massiv2, y_Array2, U_massiv2, CountDots); //4.30
    table_in_file(&Func3, x_Array, CountDots, U_massiv,U_massiv2);


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

    cout << endl << "Полученная СЛАУ: " << endl;
    for (i = 0; i < CountDots; i++) {
        for (int j = 0; j < CountDots; j++) {
            cout << coeffs_massiv2[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 << endl << "U_i2: " << endl;
    for (j = 0; j < CountDots; j++) {
        cout << U_massiv2[j] << " ";
    }

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