lab3

using namespace std;

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>

template <typename T, int N>
struct vec {
    T data[N];
    T &operator[](int i) {
        if (i >= N) throw "Out of range";
        return data[i];
    }
    const T &operator[](int i) const {
        if (i >= N) throw "Out of range";
        return data[i];
    }

    vec<T, N> operator*(const float &value) {
        vec<T, N> result = { 0 };
        for (int i = 0; i < N; i++) {
            result[i] = data[i] * value;
        }
        return result;
    }

    T dot(const vec<T, N> &value) {
        T result = { 0 };
        for (int i = 0; i < N; i++) {
            result += data[i] * value[i];
        }
        return result;
    }

    vec<T, N> &operator *=(const float &value) {
        return *this = operator*(value);
    }

    vec<T, N> &operator +=(const vec<T, N> &value) {
        return *this = operator+(value);
    }

    vec<T, N> operator+(const vec<T, N> &value) {
        vec<T, N> result;
        for (int i = 0; i < N; i++) {
            result[i] = data[i] + value[i];
        }
        return result;
    }

    vec<T, N> operator-(const vec<T, N> &value) {
        vec<T, N> result;
        for (int i = 0; i < N; i++) {
            result[i] = data[i] - value[i];
        }
        return result;
    }
};

typedef vec<double, 4> vec4;
typedef vec<vec4, 4> mat4x4;

typedef vec<double, 5> vec5;
typedef vec<vec5, 4> mat4x5;

void print(const mat4x5 &matrix) {
    for (int i = 0; i < 4; i++) {
        for (int j = 0; j < 5; j++) {
            printf("%-10g", matrix[i][j]);
        }
        printf("\n");
    }
}

void complete_pivoting(mat4x5 matrix) {
    printf("complete pivoting\n");

    mat4x5 temp = { 0 };
    for (int k = 0, p, q; k < 3; k++) {
        double a = 0;

        for (int i = 0; i < 4; i++) {
            for (int j = 0; j < 4; j++) {
                float m = fabs(matrix[i][j]);
                if (m > a) {
                    a = m;
                    p = i;
                    q = j;
                }
            }
        }

        temp[k] = matrix[p];

        for (int i = 0; i < 4; i++) {
            if (i == p) continue;

            matrix[i] += matrix[p] * (-matrix[i][q] / a);
        }

        matrix[p] = { 0 };
        for (int i = 0; i < 4; i++) {
            matrix[i][q] = 0;
        }
    }

    for (int i = 0; i < 4; i++) {
        if (matrix[i][4] != 0) {
            temp[3] = matrix[i];
            break;
        }
    }

    vec4 x = { 0 };
    for (int i = 3; i >= 0; i--) {
        double R = temp[i][4];
        for (int j = 0; j < 4; j++) {
            if ((temp[i][j] != 0) && (x[j] != 0)) {
                R -= temp[i][j] * x[j];
            }
        }
        for (int j = 0; j < 4; j++) {
            if ((temp[i][j] != 0) && (x[j] == 0)) {
                x[j] = R / temp[i][j];
            }
        }
    }

    for (int i = 0; i < 4; i++) {
        printf("x%i = %e\n", i, x[i]);
    }
}

void gauss_seidel(const mat4x5 &matrix) {
    printf("gauss-seidel\n");
    print(matrix);

    vec4 beta;
    for (int i = 0; i < 4; i++) {
        beta[i] = matrix[i][4] / matrix[i][i];
    }

    mat4x4 alfa;
    for (int i = 0; i < 4; i++) {
        for (int j = 0; j < 4; j++) {
            alfa[i][j] = (i != j) ? (-matrix[i][j] / matrix[i][i]) : 0;
        }
    }

    double q = 0;
    for (int i = 0; i < 4; i++) {
        double s = 0;
        for (int j = 0; j < 4; j++) {
            s += fabs(alfa[i][j]);
        }
        if (s > q) q = s;
    }

    double t = fabs((1 - q) * 1e-8 / q);
    printf("|      x1      |      x2      |      x3      |      x4      |\n");

    vec4 X = { 0 };
    while (true) {
        vec4 x0 = X;
        for (int i = 0; i < 4; i++) {
            X[i] = beta[i] + alfa[i].dot(X);
        }

        printf("| %12.5e | %12.5e | %12.5e | %12.5e |\n", X[0], X[1], X[2], X[3]);

        double err = 0;
        for (int i = 0; i < 4; i++) {
            double delta = fabs(x0[i] - X[i]);
            if (delta > err) {
                err = delta;
            }
        }
        if (err <= t) break;
    }
}

int main() {
    mat4x5 a = {
        10,  5,  6, 16, 72,
        19, 38, 18,  0, 76,
        16, 19, 53, 17, 98,
        10, 14, 14,  4, 48
    };

    mat4x5 b = { a[3], a[1], a[2], a[0] };

    complete_pivoting(a);
    printf("\n");
    gauss_seidel(b);
    return 0;
}