import numpy

eps = 2**-52

def LU(A):
	U = numpy.matrix(A)
	n = U.shape[0]
	assert len(U.shape)==2
	assert n==U.shape[1]
	L = numpy.matrix(numpy.identity(n))
	for j in range(n-1):
		assert abs(U[j,j])>=eps
		for i in range(j+1, n):
			mult=U[i,j]/U[j,j]
			for k in range(n):
				U[i,k] -= mult*U[j,k]
			L[i,j]=mult
	return L, U
	
def backsolve(A, b):
	b = numpy.matrix(b).transpose()
	L, U = LU(A)
	n = L.shape[0]
	assert b.shape[0]==n
	assert b.shape[1]==1
	c = numpy.matrix(numpy.zeros((n, 1)))
	x = c.copy()
	for i in range(n):
		mult = c[i] = b[i]
		for k in range(i,n):
			b[k]-= mult*L[k,i]
	for i in range(n-1, -1, -1):
		mult = x[i] = c[i]/U[i,i]
		for k in range(i):
			c[k]-=mult*U[k,i]
	return x