#include <iostream>
#include <vector>
#include <cmath>
#include <stdexcept>
#include <iomanip>
#include <sstream>
using namespace std;
vector<vector<double>> gauss_jordan(vector<vector<double>> A, vector<double> B) {
size_t n = [Link]();
vector<vector<double>> M(n, vector<double>(n + 1));
for (size_t i = 0; i < n; i++) {
for (size_t j = 0; j < n; j++) {
M[i][j] = A[i][j];
}
M[i][n] = B[i];
}
for (size_t i = 0; i < n; i++) {
size_t max_row = i;
for (size_t j = i + 1; j < n; j++) {
if (abs(M[j][i]) > abs(M[max_row][i])) {
max_row = j;
}
}
if (max_row != i) {
swap(M[i], M[max_row]);
}
if (abs(M[i][i]) < 1e-10) {
throw runtime_error("Matrix is singular or nearly singular");
}
double pivot = M[i][i];
for (size_t j = i; j < n + 1; j++) {
M[i][j] /= pivot;
}
for (size_t j = 0; j < n; j++) {
if (i != j) {
double factor = M[j][i];
for (size_t k = i; k < n + 1; k++) {
M[j][k] -= factor * M[i][k];
}
}
}
}
vector<double> solution(n);
for (size_t i = 0; i < n; i++) {
solution[i] = M[i][n];
}
return vector<vector<double>>{solution};
}
vector<vector<double>> matrix_multiply(vector<vector<double>> A,
vector<vector<double>> B) {
size_t n = [Link]();
size_t m = B[0].size();
size_t p = [Link]();
vector<vector<double>> result(n, vector<double>(m, 0));
for (size_t i = 0; i < n; i++) {
for (size_t j = 0; j < m; j++) {
for (size_t k = 0; k < p; k++) {
result[i][j] += A[i][k] * B[k][j];
}
}
}
return result;
}
vector<vector<double>> matrix_inverse(vector<vector<double>> A) {
size_t n = [Link]();
vector<vector<double>> M(n, vector<double>(2 * n, 0));
for (size_t i = 0; i < n; i++) {
for (size_t j = 0; j < n; j++) {
M[i][j] = A[i][j];
}
M[i][n + i] = 1;
}
for (size_t i = 0; i < n; i++) {
size_t max_row = i;
for (size_t j = i + 1; j < n; j++) {
if (abs(M[j][i]) > abs(M[max_row][i])) {
max_row = j;
}
}
if (max_row != i) {
swap(M[i], M[max_row]);
}
if (abs(M[i][i]) < 1e-10) {
throw runtime_error("Matrix is singular, cannot compute inverse");
}
double pivot = M[i][i];
for (size_t j = 0; j < 2 * n; j++) {
M[i][j] /= pivot;
}
for (size_t j = 0; j < n; j++) {
if (i != j) {
double factor = M[j][i];
for (size_t k = 0; k < 2 * n; k++) {
M[j][k] -= factor * M[i][k];
}
}
}
}
vector<vector<double>> inverse(n, vector<double>(n));
for (size_t i = 0; i < n; i++) {
for (size_t j = 0; j < n; j++) {
inverse[i][j] = M[i][n + j];
}
}
return inverse;
}
vector<vector<double>> inversion_method(vector<vector<double>> A, vector<double> B)
{
vector<vector<double>> A_inv = matrix_inverse(A);
vector<vector<double>> B_col([Link](), vector<double>(1));
for (size_t i = 0; i < [Link](); i++) {
B_col[i][0] = B[i];
}
vector<vector<double>> X_col = matrix_multiply(A_inv, B_col);
vector<vector<double>> solution(1, vector<double>(X_col.size()));
for (size_t i = 0; i < X_col.size(); i++) {
solution[0][i] = X_col[i][0];
}
return solution;
}
int main() {
cout << "System of Linear Equations Solver" << endl;
cout << "========================================" << endl;
int n;
cout << "Enter the number of variables: ";
cin >> n;
vector<vector<double>> A(n, vector<double>(n));
cout << "\nEnter the coefficient matrix A (" << n << "x" << n << "):" << endl;
for (int i = 0; i < n; i++) {
cout << "Row " << i + 1 << ": ";
for (int j = 0; j < n; j++) {
cin >> A[i][j];
}
}
vector<double> B(n);
cout << "\nEnter the constant vector B (" << n << " elements):" << endl;
for (int i = 0; i < n; i++) {
cin >> B[i];
}
cout << "\n========================================" << endl;
try {
vector<vector<double>> gj_solution = gauss_jordan(A, B);
cout << "Gauss-Jordan Method Solution:" << endl;
for (size_t i = 0; i < gj_solution[0].size(); i++) {
cout << "x" << i + 1 << " = " << fixed << setprecision(6) <<
gj_solution[0][i] << endl;
}
} catch (const runtime_error& e) {
cout << "Gauss-Jordan Method failed: " << [Link]() << endl;
}
cout << "\n----------------------------------------" << endl;
try {
vector<vector<double>> inv_solution = inversion_method(A, B);
cout << "Inversion Method Solution:" << endl;
for (size_t i = 0; i < inv_solution[0].size(); i++) {
cout << "x" << i + 1 << " = " << fixed << setprecision(6) <<
inv_solution[0][i] << endl;
}
} catch (const runtime_error& e) {
cout << "Inversion Method failed: " << [Link]() << endl;
}
return 0;
}