#include <iostream>
#include <vector>
#include <functional>
#include <fstream>
#include <cmath>

using namespace std;

// CONSTANTS
const int PRINT_PRECISION = 5;
const long double NOT_ZERO_EPS = 1e-16;
const long double SEIDEL_EPS = 1e-10;

struct matrix { // matrix class
    vector < vector <long double> > data;
    matrix(int n = 0, int m = -1, long double default_value = 0) {
        if (m < 0) {
            m = n;
        }
        data.resize(n);
        for (int i = 0; i < n; i++) {
            data[i].resize(m, default_value);
        }
    }
    matrix (vector <vector <long double> > _data) {
        data = _data;
    }
    int width() {
        if (!data.size()) {
            return 0;
        }
        return (int)data[0].size();
    }
    int height() {
        return (int)data.size();
    }
    void fill(function<long double(int, int)> generator) {
        for (int i = 0; i < height(); i++) {
            for (int j = 0; j < width(); j++) {
                data[i][j] = generator(i, j);
            }
        }
    }
    void print(int extra_lines = 0) {
        for (int i = 0; i < height(); i++) {
            for (int j = 0; j < width(); j++) {
                cout.precision(PRINT_PRECISION);
                cout << fixed << data[i][j] << " ";
            }
            cout << endl;
        }
        for (int i = 0; i < extra_lines; i++) {
            cout << endl;
        }
    }
    matrix transpose() {
        matrix ans(width(), height());
        for (int i = 0; i < height(); i++) {
            for (int j = 0; j < width(); j++) {
                ans.data[j][i] = data[i][j];
            }
        }
        return ans;
    }
    matrix multiply(matrix A, matrix B) {
        if (A.width() != B.height()) {
            cerr << "Multiply faild. " << A.height() << "x" << A.width() << " " <<
                                          B.height() << "x" << B.width() << endl;
            return matrix();
        }
        matrix C(A.height(), B.width());
        for (int i = 0; i < A.height(); i++) {
            for (int j = 0; j < B.width(); j++) {
                long double sum = 0;
                for (int k = 0; k < A.width(); k++) {
                    sum += A.data[i][k] * B.data[k][j];
                }
                C.data[i][j] = sum;
            }
        }
        return C;
    }
    void swap_lines(int first, int second) {
        if (first >= 0 && first < height() && second >= 0 && second < height()) {
            swap(data[first], data[second]);
        } else {
            cerr << "(swap lines) Bad lines: " << first << " " << second << endl;
        }
    }
    void swap_columns(int line, int first, int second) {
        if (line >= 0 && line < height() &&
                first >= 0 && second >= 0 &&
                first < width() && second < width()) {
            for (int i = 0; i < height(); i++) {
                swap(data[line][first], data[line][second]);
            }
        }
    }
    void multiply_line(int line, long double x) {
        if (line >= 0 && line < height()) {
            for (int i = 0; i < width(); i++) {
                data[line][i] *= x;
            }
        } else {
            cerr << "(multiply line) Bad line: " << line << endl;
        }
    }
    void add_line(int source, int destination, long double x = 1) {
        if (source >= 0 && source < height() && destination >= 0 && destination < height()) {
            for (int i = 0; i < width(); i++) {
                data[destination][i] += data[source][i] * x;
            }
        } else {
            cerr << "(add line) Bad lines: " << source << " " << destination << endl;
        }
    }
    void triangulate(matrix &twin, bool up = false, bool pivoting = false) {
        vector <pair <int, pair <int, int> > > swaps;
        for (int i = 0; i < min(width(), height()); i++) {
            int x = up ? height() - 1 - i : i;
            for (int j = i; j < height(); j++) {
                int y = up ? height() - 1 - j : j;
                if (fabs(data[y][x]) > NOT_ZERO_EPS) {
                    swap_lines(x, y);
                    if (twin.width()) {
                        twin.swap_lines(x, y);
                    }
                    break;
                }
            }
            if (pivoting) {
                int maxi = i;
                for (int j = i + 1; j < width(); j++) {
                    if (fabs(data[i][j]) > fabs(data[i][maxi])) {
                        maxi = j;
                    }
                }
                if (maxi != i) {
                    swap_columns(i, i, maxi);
                    swaps.push_back(make_pair(i, make_pair(i, maxi)));
                }
            }
            if (fabs(data[x][x]) > NOT_ZERO_EPS) {
                for (int j = i + 1; j < height(); j++) {
                    int y = up ? height() - 1 - j : j;
                    long double delta = -data[y][x] / data[x][x];
                    add_line(x, y, delta);
                    if (twin.width()) {
                        twin.add_line(x, y, delta);
                    }
                }
            }
        }
        if (pivoting) {
            for (int i = (int)swaps.size() - 1; i >= 0; i--) {
                int line = swaps[i].first;
                int x = swaps[i].second.first;
                int y = swaps[i].second.second;
                swap_columns(line, x, y);
            }
        }
    }
    void triangulate(bool up = false, bool pivoting = false) {
        matrix tmp;
        triangulate(tmp, up, pivoting);
    }
    long double determinant() {
        matrix tmp(data);
        tmp.triangulate();
        long double ans = 1;
        for (int i = 0; i < min(tmp.width(), tmp.height()); i++) {
            ans *= tmp.data[i][i];
        }
        return ans;
    }
    void normalize_diag(matrix &twin) {
        for (int i = 0; i < min(width(), height()); i++) {
            if (abs(data[i][i]) > NOT_ZERO_EPS) {
                if (twin.width()) {
                    twin.multiply_line(i, 1. / data[i][i]);
                }
                multiply_line(i, 1. / data[i][i]);
            }
        }
    }
    void normalize_diag() {
        matrix tmp;
        normalize_diag(tmp);
    }
    matrix reverse_matrix() {
        matrix tmp(data);
        if (fabs(tmp.determinant()) < NOT_ZERO_EPS) {
            cerr << "Error: Cannot find reverse matrix. det = 0" << endl;
            return matrix();
        }
        matrix ans(tmp.height(), tmp.width());
        ans.fill([](int i, int j) {
            return 1 ? i == j : 0;
        });
        tmp.triangulate(ans);
        tmp.triangulate(ans, 1);
        tmp.normalize_diag(ans);
        return ans;
    }
    matrix slow_solve_equation(matrix B, bool pivoting = false) { // cool, but slow Gauss
        matrix tmp(data);
        tmp.triangulate(B, 0, pivoting);
        tmp.triangulate(B, 1, pivoting);
        tmp.normalize_diag(B);
        return B;
    }
    matrix solve_equation(matrix B, bool pivoting = false) { // ordinary Gauss
        matrix tmp(data);
        tmp.triangulate(B, 0, pivoting);
        tmp.normalize_diag(B);
        matrix res(B);
        for (int i = height() - 1; i >= 0; i--) {
            for (int j = i + 1; j < width(); j++) {
                B.data[i][0] -= tmp.data[i][j] * res.data[j][0];
            }
            res.data[i][0] = B.data[i][0];
        }
        return res;
    }
    matrix solve_equation_seidel(matrix B, long double w = 1) { // Seidel or iteration algorithm
        matrix T = matrix(data).transpose();
        matrix A = matrix(data);
        A = matrix().multiply(T, A);
        B = matrix().multiply(T, B);
        matrix res(B.height(), 1, 0);
        matrix p;
        bool convenient = false;
        while (!convenient) {
            p = matrix(res);
            for (int i = 0; i < A.width(); i++) {
                long double sum = 0;
                for (int j = 0; j < i; j++) {
                    sum += A.data[i][j] * res.data[j][0];
                }
                for (int j = i; j < A.width(); j++) {
                    sum += A.data[i][j] * p.data[j][0];
                }
                res.data[i][0] = p.data[i][0] + w * (B.data[i][0] - sum) / A.data[i][i];
            }
            long double sum = 0;
            for (int i = 0; i < width(); i++) {
                sum += pow(res.data[i][0] - p.data[i][0], 2);
            }
            convenient = sum < SEIDEL_EPS;
        }
        return res;
    }
};

int main(int argc, char *argv[])
{
    string data_file;
    string result_file;
    string action_type;
    int N;
    long double x = 0, M = 0;
    if (argc == 5) {
        N = strtoll(argv[1], NULL, 10);
        action_type = string(argv[2]);
        if (action_type == "file") {
            data_file = string(argv[3]);
            result_file = string(argv[4]);
        } else if (action_type == "formula") {
            x = atof(argv[3]);
            M = atof(argv[4]);
        } else {
            cout << "Wrong usage: N Type(file/formula) [[data_file_name result_file_name], [x, M]]" << endl;
            return 0;
        }
    } else {
        cout << "Wrong usage: N Type(file/formula) [[data_file_name result_file_name], [x, M]]" << endl;
        return 0;
    }
    cout << "Start" << endl;
    cout << "N: " << N << endl;
    matrix A(N, N), B(N, 1);
    if (action_type == "file") {
        cout << "Data file: " << data_file << endl;
        cout << "Result file: " << result_file << endl;
        auto data_fill = [=](int, int) {
            static ifstream fin(data_file);
            long double ans;
            fin >> ans;
            return ans;
        };
        auto res_fill = [=](int, int) {
            static ifstream fin(result_file);
            long double ans;
            fin >> ans;
            return ans;
        };
        A.fill(data_fill);
        B.fill(res_fill);
    } else {
        cout << "M: " << M << endl;
        cout << "x: " << x << endl;
        auto A_gen = [=](int i, int j) {
            long double q = 1.001 - 2 * M * pow(10, -3);
            if (i == j) {
                return pow(q - 1, i + j + 2);
            } else {
                return pow(q, i + j + 2) + 0.1 * (j - i);
            }
        };
        auto B_gen = [=](int i, int){
            return  N * exp(x / (i + 1)) * cos(x);
        };
        A.fill(A_gen);
        B.fill(B_gen);
    }
    cout << "Data:" << endl;
    A.print(1);
    cout << "Answers:" << endl;
    B.print(1);
    cout << "Determinant: " << A.determinant() << endl;
    if (abs(A.determinant()) < NOT_ZERO_EPS) {
        cerr << "Error: det = 0" << endl;
        return 0;
    }
    cout << "Reverse matrix:" << endl;
    A.reverse_matrix().print(1);
    cout << "Gauss solution:" << endl;
    A.solve_equation(B, 1).print(1);
    cout << "Seidel solution:" << endl;
    A.solve_equation_seidel(B, 1.99).print(1);
    return 0;
}