tridiagonal matrix algorithm

#include<bits/stdc++.h>

#define INF 1000010000
#define nl '\n'
#define pb push_back
#define ppb pop_back
#define mp make_pair
#define fi first
#define se second
#define pii pair<int,int>
#define pdd pair<double,double>
#define all(c) (c).begin(), (c).end()
#define SORT(c) sort(all(c))
#define sz(c) (c).size()
#define rep(i,n) for( int i = 0; i < n; ++i )
#define repi(i,n) for( int i = 1 ; i <= n; ++i )
#define repn(i,n) for( int i = n - 1 ; i >= 0 ; --i )
#define repf(j,i,n) for( int j = i ; j < n ; ++j )
#define die(s) {std::cout << s << nl;}
#define dier(s) {std::cout << s; return 0;}
#define vi vector<int>
typedef long long ll;
#define MATRIX vector<vector<double> >
#define MAXN 5050
using namespace std;
int n;
vi dj = {-1 , 0 , 1};

inline vector<double> Tridiagonal(MATRIX A , vector<double> b){
    vector<double> ans(n , 0);
    vector<double> alpha(n , 0) , beta(n , 0);
    alpha[1] = -A[0][1] / A[0][0];
    beta[1] = b[0] / A[0][0];
    repi(i , n - 2){
        alpha[i + 1] = -A[i][i + 1] / (A[i][i - 1] * alpha[i] + A[i][i]);
        beta[i + 1] = (b[i] - (A[i][i - 1] * beta[i])) / (A[i][i - 1] * alpha[i] + A[i][i]);
    }
    ans[n - 1] = (b[n - 1] - (A[n - 1][n - 2] * beta[n - 1])) / (A[n - 1][n - 1] + (A[n - 1][n - 2] * alpha[n - 1]));
    repn(i , n - 1){
        ans[i] = (alpha[i + 1] * ans[i + 1]) + beta[i + 1];
    }
    return ans;
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);
    cout.precision(3);
    cout.setf(ios::fixed);
    cin >> n;
    srand(std::chrono::duration_cast<std::chrono::nanoseconds>((std::chrono::high_resolution_clock::now().time_since_epoch())).count() + time(0));
    MATRIX A(n , vector<double>(n , 0.));
    rep(i , n){
        rep(j , 3){
            if(i + dj[j] < 0 || i + dj[j] >= n)continue;
            //cin >> A[i][i + dj[j]];
            double tmp = (double)(rand() % 100) / (double)((rand() % 100) + 1);
            if(rand() % 2)tmp = -tmp;
            A[i][i + dj[j]] = tmp;
        }
    }
    vector<double> X(n , 0);
    rep(i , n){
    //  cin >> X[i];
        double tmp = (double)(rand() % 100) / (double)((rand() % 100) + 1);
        if(rand() % 2)tmp = -tmp;
        X[i] = tmp;
    }
    vector<double> B(n , 0);
    rep(i , n){
        rep(j , 3){
            if(i + dj[j] < 0 || i + dj[j] >= n)continue;
            B[i] += A[i][i + dj[j]] * X[i + dj[j]];
        }
    }
    rep(i , n){
        rep(j , n){
            if(j && A[i][j] >= 0)cout << "+ ";
            else cout << "  ";
            cout << setw(6) << A[i][j] << "x" << j + 1 << " ";
        }
        cout << " = " << B[i] << nl;
    }
    vector<double> ans = Tridiagonal(A , B);
    cout << nl;
    rep(i , n)cout << setw(4) << "x" << i + 1 << " = " << ans[i] << " ";
    return 0;
}