#include <vector>
#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;

// Multiplication A * x
vector<double> Multiplication(vector<double> low, vector<double> main, vector<double> up, vector<double> x) {
    vector<double> result(x.size());
    for (int i = 0; i < x.size(); i++) {
        result[i] = main[i] * x[i];
        if (i != 0) {
            result[i] += low[i] * x[i - 1];
        }
        if (i != x.size() - 1) {
            result[i] += up[i] * x[i + 1];
        }
    }
    return result;
}

// Norm(lhs - rhs)
double Norm(vector<double> lhs, vector<double> rhs) {
    double ttl = 0;
    for (int i = 0; i < rhs.size(); i++) {
        ttl += (lhs[i] - rhs[i]) * (lhs[i] - rhs[i]);
    }
    return sqrt(ttl);
}

// Matvec for tay: tay * A * x
vector<double> Matvec(vector<double> low, vector<double> main, vector<double> up, vector<double> x, double tay) {
    vector<double> result = Multiplication(low, main, up, x);
    for (int i = 0; i < x.size(); i++) {
        result[i] *= tay;
    }
    return result;
}

vector<double> Matvec(vector<double> low, vector<double> main, vector<double> up, vector<double> x) {
    vector<double> d_1(main.size()), null(main.size());
    for (int i = 0; i < main.size(); i++) {
        d_1[i] = 1 / main[i];
    }
    return Multiplication(null, d_1, null, Multiplication(low, main, up, x));
}

// SI_Solver for MTI
void SI_Solver(vector<double> low, vector<double> main, vector<double> up, vector<double> f, vector<double> res, double tay) {
    // Norm (Ax - f)
    vector<double> norm_vec;
    double norm0 = Norm(Multiplication(low, main, up, res), f);
    ofstream mti("MTI.txt"), mti_answ("MTI_answ.txt");
    while (1) {
        vector<double> v = Matvec(low, main, up, res, tay);
        double norm = Norm(v, f);
        if (norm / norm0 < 1e-4) {
            break;
        }
        for (int i = 0; i < res.size(); i++) {
            res[i] += f[i] - v[i];
        }
        mti << norm / norm0 << "\n";
    }
    for (int i = 0; i < res.size(); i++) {
        mti_answ << res[i] << "\n";
    }
}

// SI_Solver for Jacobi
void SI_Solver(vector<double> low, vector<double> main, vector<double> up, vector<double> f, vector<double> res) {
    double norm0 = Norm(Multiplication(low, main, up, res), f);
    ofstream jac("jacobi.txt"), jac_answ("jacobi_answ.txt");
    while (1) {
        vector<double> v = Matvec(low, main, up, res);
        double norm = Norm(v, f);
        jac << norm / norm0 << "\n";
        if (norm / norm0 < 1e-4) {
            break;
        }
        for (int i = 0; i < res.size(); i++) {
            res[i] += f[i] - v[i];
        }
    }
    for (int i = 0; i < res.size(); i++) {
        jac_answ << res[i] << "\n";
    }

}

int main()
{
    double tay = 0.1651512666112257;
    vector<double> up(1000, -1);
    vector<double> low(1000, -1);
    vector<double> res(1000);
    vector<double> main(1000, 2.01);
    vector<double> f(1000), ft(1000), d_1f(1000);
    main[0] = 12;
    for (int i = 494; i < 505; i++) {
        f[i] = 1;
        ft[i] = tay;
        d_1f[i] = 1 / main[i];
    }

    // MTI
    SI_Solver(low, main, up, ft, res, tay);

    //Jacobi
    SI_Solver(low, main, up, d_1f, res);
}