#include <iostream>
using namespace std;
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <fstream>
#include <math.h>
#include <iomanip>
#define PI  3.14159265358979323846

using namespace std;

void fill_matrix(float **matrix, int n, float t, float h, float *vector_u, float *vector_v, float *vector_b) {
    for (int i = 0; i < n; i ++) {
        for (int j = 0; j < n; j ++) {
             matrix[i][j] = 0;
        }
    }
    for (int i = 0; i < n; i ++) {
        //столбец b
        if (i < n/2) {
            vector_b[i] = vector_v[i]/t;
        } else {
            vector_b[i] = vector_u[i-(n/2)+1]/t;
        }
        for (int j = 0; j < n; j ++) {
            //диагональ
            if (i == j) {
                matrix[i][j] = 1/t;
            }
            if (i > 0 && i < (n/2)-1) {
                matrix[i][i-1] = -vector_u[i]/h;
                matrix[i][i+1] = vector_u[i]/h;
            }
            if (i == (n/2)-1) {
                matrix[i][i-1] = -vector_u[i]/h;
            }
            if (i == n/2) {
                matrix[i][i+1] = vector_v[i-(n/2)+1]/h;
            }
            if (i > (n/2) && i < n-1) {
                matrix[i][i-1] = -vector_v[i-(n/2)+1]/h;
                matrix[i][i+1] = vector_v[i-(n/2)+1]/h;
            }
            //поддиагональная диагональ :D
            if (i >= n/2 && i < n-2) {
                matrix[i][i-(n/2)] = -(vector_u[i-(n/2)+1]+vector_v[i-(n/2)+1])/h;
                matrix[i][i-(n/2)+2] = (vector_u[i-(n/2)+1]+vector_v[i-(n/2)+1])/h;
            }
            if (i == n-2) {
                matrix[i][i-(n/2)] = -(vector_u[i-(n/2)+1]+vector_v[i-(n/2)+1])/h;
            }
            if (i == n-1) {
                matrix[i][i-(n/2)-1] = vector_u[n/2]/h;
                matrix[i][i-(n/2)] = -4*vector_u[n/2]/h;
            }
            //наддиагональная диагональ :D
            if (i == 0) {
                matrix[i][i+(n/2)] = 4*vector_v[i]/h;
                matrix[i][i+(n/2)+1] = -vector_v[i]/h;
            }
            if (i == 1) {
                matrix[i][i+(n/2)] = (vector_u[i]+vector_v[i])/h;
            }
            if (i > 1 && i < (n/2)) {
                matrix[i][i+(n/2)-2] = -(vector_u[i]+vector_v[i])/h;
                matrix[i][i+(n/2)] = (vector_u[i]+vector_v[i])/h;
            }
        }
    }
    return;
}

std::ofstream ofs ("output.txt", std::ofstream::out);

void print_matrix(float **matrix, int n) {
    for (int i = 0; i < n; i ++) {
        for (int j = 0; j < n; j ++) {
            ofs << matrix[i][j] << setw(13) << ' ' ;
        }
        ofs << endl;
    }
    ofs << endl;
    return;
}

void print_vector(float *p_vector, int n) {
    for (int i = 0; i < n; i ++) {
        ofs << p_vector[i] << ' ' ;
    }
    ofs << endl;
}

float** matrix_return(float cos, float sin, int n, int c, int u) {
    int i, j;
    float** mat = (float**)malloc(n * sizeof(float*));
    for (i = 0; i < n; i ++) {
        mat[i] = (float*)malloc(n * sizeof(float));
    }
    for (i = 0; i < n; i ++) {
        for (j = 0; j < n; j ++) {
            mat[i][j] = 0;
        }
        mat[i][i] = 1;
    }
    mat[c][c] = cos;
    mat[c][u+1] = -sin;
    mat[u+1][c] = sin;
    mat[u+1][u+1] = cos;
    return (mat);
}

float** multiplication(float **mat1, float **mat2, int n) {
    int i, j, k;
    float** mat3 = (float**)malloc(n * sizeof(float*));
    for (i = 0; i < n; i ++) {
        mat3[i] = (float*)malloc(n * sizeof(float));
    }
    for (i = 0; i < n; i ++) {
        for (j = 0; j < n; j ++) {
            mat3[i][j] = 0;
        }
    }
    for (i = 0; i < n; i ++) {
        for (j = 0; j < n; j ++) {
            for (k = 0; k < n; k ++) {
                mat3[i][j] += mat1[i][k] * mat2[k][j];
            }
        }
    }
    return (mat3);
}

float* multiplication1(float **mat1, float *m_vector, int n) {
    int i, j, k;
    float* mat3 = (float*)malloc(n*sizeof(float));
    for (i = 0; i < n; i ++) {
        mat3[i] = 0;
    }
    for (i = 0; i < n; i ++) {
        for (k = 0; k < n; k ++) {
            mat3[i] += mat1[i][k] * m_vector[k];
        }
    }
    return(mat3);
}

float** fill_matrix_for_reverse(float **buf, float **single, float **res, int n, int c) {
    int i, j, k;
    float t, sum;
    float** buf1 = (float**)malloc(n * sizeof(float*));
    for (i = 0; i < n; i ++) {
        buf1[i] = (float*)malloc(n * sizeof(float));
    }
    for (i = 0; i < n; i ++) {
        for (j = 0; j < n; j++) {
            buf1[i][j] = buf[i][j];
        }
    }
    for (j = n - 1; j >= 0; j --) {
        sum = 0;
        for (k = 0; k < n; k ++) {
            if (k != j) {
                sum += buf1[j][k];
            }
        }
        t = (single[j][c] - sum )/ buf1[j][j];
        res[j][c] = t;
        for (i = j; i >= 1; i --) {
            buf1[i-1][j] *= t;
        }
    }
    return res;
}

float** reverse_matrix(float **buf, int n) {
    int i, j;
    float** buf_new = (float**)malloc(n * sizeof(float*));
    for (i = 0; i < n; i ++) {
        buf_new[i] = (float*)malloc(n * sizeof(float));
    }
    for (i = 0; i < n; i ++) {
        for (j = 0; j < n; j++) {
            buf_new[i][j] = buf[i][j];
        }
    }
    float** single = (float**)malloc(n * sizeof(float*));
    for (i = 0; i < n; i ++) {
        single[i] = (float*)malloc(n * sizeof(float));
    }
    for (i = 0; i < n; i ++) {
        for (j = 0; j < n; j ++) {
            single[i][j] = 0;
        }
        single[i][i] = 1;
    }
    float** matrix;
    int k, t, p;
    p = 0;
    for(j = 0; j < n; j ++) {
        t = j + 1;
        k = 0;
        for(i = j; i < n - 1; i ++) {
            float x = buf[j][j];
            float y = buf[t][j];
            float cos = (x / sqrt(x*x + y*y));
            float sin = - (y / sqrt(x*x + y*y));
            matrix = matrix_return(cos, sin, n, j, k+p);
            buf = multiplication(matrix, buf, n);
            single = multiplication(matrix, single, n);
            t ++;
            k ++;
        }
        p ++;
    }
    float** result = (float**)malloc(n * sizeof(float*));
    for (i = 0; i < n; i ++) {
        result[i] = (float*)malloc(n * sizeof(float));
    }
    for(j = 0; j < n; j ++) {
        result = fill_matrix_for_reverse(buf, single, result, n, j);
    }
    return (result);
}

int main()
{
    int i, j, N, K, k;
    float r;
    printf("%s\n", "Input N:");
    scanf("%d", &N);
    printf("%s\n", "Input K:");
    scanf("%d", &K);
    float t = 1.0/N;
    float h = 1.0/K;
    cout << "t: " << t << endl;
    cout << "h: " << h << endl;
    ofs << "N: " << N << endl;
    ofs << "K: " << K << endl;
    ofs << "t: " << t << endl;
    ofs << "h: " << h << endl;

    float *vector_u;
    vector_u = (float*)malloc(K*sizeof(float));
    float *vector_v;
    vector_v = (float*)malloc(K*sizeof(float));
    float *vector_b;
    vector_b = (float*)malloc(2*K*sizeof(float));

    for (i = 1; i <= K; i ++) {
        vector_u[i] = (float)i/K;
    }
    vector_u[0] = 0;

    for (i = 0; i < K; i ++) {
        vector_v[i] = (float)(K-i)/K;
    }
    vector_v[K] = 0;

    //print_vector(vector_u, K+1);
    //print_vector(vector_v, K+1);
    float **matrix;
    matrix = (float**)malloc(2*K*sizeof(float*));
	for(int i = 0; i < 2*K; i ++) {
		matrix[i] = (float*)malloc(2*K*sizeof(float));
	}
	ofs << "Integral: " << endl;
    for (k = 0; k <= N; k ++) {
        float integral = 0;
        for (i = 1; i <= K-1; i ++) {
            integral += (vector_u[i]*vector_u[i]+vector_v[i]*vector_v[i]);
        }
        integral *= h;
        integral += (h/2)*(vector_u[K]*vector_u[K]+vector_v[0]*vector_v[0]);
        //ofs << "Iteration: " << k << "; Integral: " << integral << endl;
        ofs << integral << ", ";
        fill_matrix(matrix, 2*K, t, h, vector_u, vector_v, vector_b);

        //print_vector(vector_b, 2*K);
        //print_matrix(matrix, 2*K);

        float** rev_m = reverse_matrix(matrix, 2*K);
        //print_matrix(rev_m, 2*K);
        float* res_v = multiplication1(rev_m, vector_b, 2*K);
        //print_vector(res_v, 2*K);
        for (i = 1; i <= K; i ++) {
            vector_u[i] = res_v[i+K-1];
        }
        vector_u[0] = 0;
        for (i = 0; i < K; i ++) {
            vector_v[i] = res_v[i];
        }
        vector_v[K] = 0;
    }
    ofs << endl;
    ofs << "vector_u:" << endl;
    print_vector(vector_u, K+1);
    ofs << "vector_v:" << endl;
    print_vector(vector_v, K+1);

    return 0;
}