import copy
import pprint
import numpy

def gaussian_with_pivot(A, b, round = lambda x: x):
	A = numpy.array(A, dtype = numpy.float64)
	b = numpy.array(b, dtype = numpy.float64)
	L = numpy.identity(len(A), dtype = numpy.float64)
	s = numpy.zeros(len(A), dtype = numpy.float64)

	for k in range(len(A)):
		max_a = numpy.absolute(A[k:,k]).argmax()
		A[[k, k + max_a]] = A[[k + max_a, k]]
		b[[k, k + max_a]] = b[[k + max_a, k]]
		s[k] = k + max_a

		for i in range(k + 1, len(A)):
			L[i][k] = round(A[i][k]/A[k][k])
			b[i] = round(b[i] -  round(L[i][k] * b[k]))

			for j in range(k + 1, len(A)):
				A[i][j] = round(A[i][j] -  round(L[i][k] * A[k][j]))

			A[i][k] = 0
	return A, b, L, s

def gaussian_no_pivot(A, b, round = lambda x: x):
	A = numpy.array(A, dtype = numpy.float64)
	b = numpy.array(b, dtype = numpy.float64)
	L = numpy.identity(len(A), dtype = numpy.float64)

	for k in range(len(A)):

		for i in range(k + 1, len(A)):
			L[i][k] = round(A[i][k]/A[k][k])
			b[i] = round(b[i] -  round(L[i][k] * b[k]))

			for j in range(k + 1, len(A)):
				A[i][j] = round(A[i][j] -  round(L[i][k] * A[k][j]))

			A[i][k] = 0
	return A, b, L


if( __name__ == "__main__"):
	# For Exercise 1
	A = [[10, 1, 0], [15, 2, 1], [5, 0, 1]]
	b = [12, 18, 4]
	pprint.pprint(gaussian_with_pivot(A, b))
	
	
	# For Exercise 2
	A = [[7, 1, 1]
		, [10, 1, 1]
		, [1000, 0, 1]]
	b = [10, 13, 1001]
	pprint.pprint(A)
	
	def round4(x):
		return float("%.4g" % x)
	
	pprint.pprint(gaussian_with_pivot(A, b, round4))
	pprint.pprint(gaussian_no_pivot(A, b, round4))
	
	
	A = numpy.zeros((10, 10))

	for i in range(10):
		if(i - 1 > 0):
			A[i][i - 1] = -1
		A[i][i] = 2
		if(i + 1 < 10):
			A[i][i + 1] = -1

	pprint.pprint(gaussian_no_pivot(A, numpy.ones(10)))