include "primitives.alh"
include "random.alh"

Float function v2f(i : Element):
	return cast_s2f(cast_v2s(i))!

Void function dict_overwrite(d : Element, key : Element, value : Element):
	if (dict_in(d, key)):
		dict_delete(d, key)
	if (dict_in_node(d, key)):
		dict_delete_node(d, key)
	dict_add(d, key, value)

	return !

Element function create_random_matrix(n : Integer):
	Element m
	Integer i
	Integer j
	Element t

	// Construct the matrix first, with as many rows as there are variables
	// Number of columns is 1 higher
	i = 0
	m = create_node()
	while (i < n):
		j = 0
		t = create_node()
		while (j < (n + 1)):
			list_append(t, random())
			j = j + 1
		list_append(m, t)
		i = i + 1

	return m!

Void function eliminateGaussJordan(m : Element):
	Integer i
	Integer j
	Integer f
	Integer g
	Boolean searching
	Element t
	Float divisor

	i = 0
	j = 0

	while (i < read_nr_out(m)):
		// Make sure pivot m[i][j] != 0, swapping if necessary
		while (v2f(m[i][j]) == 0.0):
			// Is zero, so find row which is not zero
			f = i + 1
			searching = True
			while (searching):
				if (f >= read_nr_out(m)):
					// No longer any rows left, so just increase column counter
					searching = False
					j = j + 1
				else:
					if (v2f(m[f][j]) == 0.0):
						// Also zero, so continue
						f = f + 1
					else:
						// Found non-zero, so swap row
						t = v2f(m[f])
						dict_overwrite(m, f, v2f(m[i]))
						dict_overwrite(m, i, t)
						searching = False
			// If we have increased j, we will just start the loop again (possibly), as m[i][j] might be zero again

		// Pivot in m[i][j] guaranteed to not be 0
		// Now divide complete row by value of m[i][j] to make it equal 1
		f = j
		divisor = v2f(m[i][j])
		while (f < read_nr_out(m[i])):
			dict_overwrite(m[i], f, float_division(v2f(m[i][f]), divisor))
			f = f + 1

		// Eliminate all rows in the j-th column, except the i-th row
		f = 0
		while (f < read_nr_out(m)):
			if (bool_not(f == i)):
				g = j
				divisor = v2f(m[f][j])
				while (g < read_nr_out(m[f])):
					dict_overwrite(m[f], g, v2f(m[f][g]) - (divisor * v2f(m[i][g])))
					g = g + 1
			f = f + 1

		// Increase row and column
		i = i + 1
		j = j + 1

	return !