%ALIANJECON OTERO WHAHAHAH
clc, clear, close all;

function test(A,x,B)
if norm(A*x - B) < 1e-9
    disp("True!")
else
    disp("False!")
end
end

function x = gauss_elimination(A,B)
Ab = [A B];
[n,m]=size(Ab);
if m ~= n + 1
    error('[A B] must be augmented matrix of size n x (n+1)');
end

x = zeros(n,1);
for i = 1:n-1
    [~, maxRow] = max(abs(Ab(i:n, i)));
    maxRow = maxRow + i - 1;
    if Ab(maxRow, i) == 0
        error('singular matrix or zero column');
    end
    if maxRow ~= i
        Ab([i, maxRow], :) = Ab([maxRow, i], :);
    end

    for j = i+1:n
        %additive inverse then subtract to make zero
        m = Ab(j,i)/Ab(i,i);
        Ab(j,:)=Ab(j,:) - m * Ab(i,:);
    end
end

for i = n:-1:1
    if Ab(i,i)==0
        error('No unique sol');
    end
    %xi = bi-sum_j=i+1^n aijxj/aii xi = bi-knownterms/aii
    x(i) = (Ab(i,end) - Ab(i,i+1:n) * x(i+1:n)) / Ab(i,i);
end
end

A = [1 2 3; 2 1 1; 3 4 2];
B = [14; 10; 19];

x = gauss_elimination(A,B);

test(A,x,B)

A = [0 1 1; 1 1 1; 2 3 4];
B = [2; 3; 7];

x = gauss_elimination(A,B);

test(A,x,B)

A = [1 0 0; 0 1 0; 0 0 1];
B = [2; 3; 7];

x = gauss_elimination(A,B);

test(A,x,B)

A = [2 3 1; 5 2 2; 1 6 8];
B = [1; 2; 3];

x = gauss_elimination(A,B);

test(A,x,B)