#include <iostream>
#include <vector>
#include <cmath>
#include <string>

// Función para imprimir un vector
void printVector(const std::vector<double>& vec) {
    for (double v : vec) {
        std::cout << v << " ";
    }
    std::cout << std::endl;
}

// Función para calcular la norma infinita de la diferencia entre dos vectores
double infNorm(const std::vector<double>& a, const std::vector<double>& b) {
    double max = 0.0;
    for (size_t i = 0; i < a.size(); ++i) {
        max = std::max(max, std::abs(a[i] - b[i]));
    }
    return max;
}

// Método de Jacobi
std::vector<double> jacobi(const std::vector<std::vector<double>>& A, const std::vector<double>& b, const std::vector<double>& x0, double tol, int max_iter) {
    size_t n = b.size();
    std::vector<double> x = x0;
    std::vector<double> x_old(n);

    for (int k = 0; k < max_iter; ++k) {
        x_old = x;
        
        for (size_t i = 0; i < n; ++i) {
            double sum = 0.0;
            for (size_t j = 0; j < n; ++j) {
                if (j != i) {
                    sum += A[i][j] * x_old[j];
                }
            }
            x[i] = (b[i] - sum) / A[i][i];
        }
        
        if (infNorm(x, x_old) < tol) {
            std::cout << "Jacobi method converged in " << k + 1 << " iterations." << std::endl;
            return x;
        }
    }
    
    std::cout << "Jacobi method did not converge within the maximum number of iterations." << std::endl;
    return x;
}

// Método de Gauss-Seidel
std::vector<double> gaussSeidel(const std::vector<std::vector<double>>& A, const std::vector<double>& b, const std::vector<double>& x0, double tol, int max_iter) {
    size_t n = b.size();
    std::vector<double> x = x0;

    for (int k = 0; k < max_iter; ++k) {
        std::vector<double> x_old = x;

        for (size_t i = 0; i < n; ++i) {
            double sum = 0.0;
            for (size_t j = 0; j < n; ++j) {
                if (j != i) {
                    sum += A[i][j] * x[j];
                }
            }
            x[i] = (b[i] - sum) / A[i][i];
        }

        if (infNorm(x, x_old) < tol) {
            std::cout << "Gauss-Seidel method converged in " << k + 1 << " iterations." << std::endl;
            return x;
        }
    }

    std::cout << "Gauss-Seidel method did not converge within the maximum number of iterations." << std::endl;
    return x;
}

// Función principal para resolver el sistema de ecuaciones lineales
std::vector<double> solveLinearSystem(const std::vector<std::vector<double>>& A, const std::vector<double>& b, const std::vector<double>& x0, double tol, int max_iter, const std::string& method) {
    if (method == "jacobi") {
        return jacobi(A, b, x0, tol, max_iter);
    } else if (method == "gauss-seidel") {
        return gaussSeidel(A, b, x0, tol, max_iter);
    } else {
        throw std::invalid_argument("Unknown method. Use 'jacobi' or 'gauss-seidel'.");
    }
}

int main() {
    // Definición de la matriz A y el vector b
    std::vector<std::vector<double>> A = {{4, -1}, {-1, 3}};
    std::vector<double> b = {3, 3};
    std::vector<double> x0 = {0, 0};  // Valor inicial
    double tol = 1e-6;  // Tolerancia
    int max_iter = 100;  // Máximo número de iteraciones

    // Resolver usando el método de Jacobi
    std::vector<double> x_jacobi = solveLinearSystem(A, b, x0, tol, max_iter, "jacobi");
    std::cout << "Solution using Jacobi method: ";
    printVector(x_jacobi);

    // Resolver usando el método de Gauss-Seidel
    std::vector<double> x_gauss_seidel = solveLinearSystem(A, b, x0, tol, max_iter, "gauss-seidel");
    std::cout << "Solution using Gauss-Seidel method: ";
    printVector(x_gauss_seidel);

    return 0;
}