4MO Lab2

#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;

using namespace std;
#define eps 1e-3
int n;
vector<double> X_;


inline double EuclideanNorm(vector<double> v){
    double sum = 0;
    rep(i , sz(v))sum += v[i] * v[i];
    return sqrt(sum);
}

inline vector<double> Substract(vector<double> & a , vector<double> & b){
    vector<double> res(n , 0.);
    rep(i , n)res[i] = a[i] - b[i];
    return res;
}

inline double Residual(vector<vector<double> > & A , vector<double> & X , vector<double> & B){
    vector<double> res(n , 0.);
    rep(i , n){
        double sum = 0;
        rep(j , n){
            sum += A[i][j] * X[j];
        }
        res[i] = sum - B[i];
    }
    return EuclideanNorm(res);
}

inline double AnsDiff(vector<double> & X){
    return EuclideanNorm(Substract(X , X_));
}

inline void DiagonallyDomination(vector<vector<double> > & A){
    rep(i , n){
        int sgn = A[i][i] < 0 ? -1 : 1;
        double sum = -abs(A[i][i]);
        rep(j , n){
            sum += abs(A[i][j]);
        }
        if(sum > abs(A[i][i])){
            A[i][i] = sgn * (sum + (rand() % 5));
        }
    }
}

inline void solve1(vector<vector<double> > A , vector<double> B){//Gauss
    vector<double> X(n);
    rep(j , n){
        repf(i , j + 1 , n){
        /*  int mx = j;
            rep(k , n)if(abs(A[k][j]) > abs(A[mx][j]))mx = k;
            swap(A[i] , A[mx]);
            swap(B[i] , B[mx]);*/
            double k = A[i][j] / A[j][j];
            rep(r , n){
                A[i][r] -= k * A[j][r];
            }
            B[i] -= k * B[j];
        }
    }

    repn(i , n){
        double res = B[i];
        repf(j , i + 1 , n){
            res -= A[i][j] * X[j];
        }
        res /= A[i][i];
        X[i] = res;
        if(X[i] == -0.)X[i] = 0;
    }
    rep(i , n)cout << X[i] << " ";
    cout << nl;
}

inline void solve2(vector<vector<double> > A , vector<double> B){//Jacobi
    vector<double> X = B;
    vector<double> T = B;
    int cnt = 0;
    do{
        ++cnt;
        X = T;
        rep(i , n){
            double sum = 0;
            rep(j , n){
                if(i == j)continue;
                sum += A[i][j] * X[j];
            }
            T[i] = (B[i] - sum) / A[i][i];
        }
        cout << "Step " << cnt << ", residual is " << Residual(A , T , B) << ", diff " << AnsDiff(T) << nl;
    }
    while(EuclideanNorm(Substract(X , T)) > eps);
    rep(i , n)cout << T[i] << " ";
    cout << nl;
}

inline void solve3(vector<vector<double> > A , vector<double> B){//Seidel
    vector<double> X = B;
    vector<double> T = B;
    int cnt = 0;
    do{
        ++cnt;
        X = T;
        rep(i , n){
            double sum = 0;
            rep(j , n){
                if(i == j)continue;
                sum += A[i][j] * T[j];
            }
            T[i] = (B[i] - sum) / A[i][i];
        }
        cout << "Step " << cnt << ", residual is " << Residual(A , T , B) << ", diff " << AnsDiff(T) << nl;
    }while(EuclideanNorm(Substract(X , T)) > eps);
    rep(i , n)cout << T[i] << " ";
    cout << nl;
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);
    cout.precision(5);
    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));
    vector<vector<double> > A(n , vector<double>(n));
    rep(i , n){
        rep(j , n){
            A[i][j] = ((rand() % 500) / 100.) - 2.5;
         // A[i][j] = 1. / (i + j + 1);
        }
    }
    DiagonallyDomination(A);
    rep(i , n){
        rep(j , n){
            cout << setw(6) << A[i][j] << " ";
        }
        cout << nl;
    }
    vector<double> b(n);
    rep(i , n){
        X_.pb(0);
        cin >> X_[i];
    }
    rep(i , n){
        double cur = 0;
        rep(j , n){
            cur += X_[j] * A[i][j];
        }
        b[i] = cur;
    }
    cout << "B: ";
    rep(i , n)cout << b[i] << " ";
    cout << nl;
    solve1(A , b);
    solve2(A , b);
    solve3(A , b);
    return 0;
}