Improved gaussian reducer (supports expressions like ln(pi^2), uses fractions of bignums)

mega-gaussian-elim.py

# Gaussian elimination
# -- fractional datatype version
# -- deluxe big integers and expressions edition
# 
# By blockeduser the wannabe programmer
# Original clown version: September 1-3, 2012
# Python Deluxe big integers version: September 14, 2013

import sys
import re
from math import *
from fractions import *

def n(k):
    while int(k) != k:
            k *= 10
    return int(k)

def d(k):
    d = 1
    while int(k) != k:
         k *= 10
         d *= 10
    return int(d)


# get a token from stdin
def symbol():
	def chkvalid(c):
		return c not in ['\t', '\n', '\r', 0] and c >= 0
	_str = ""
	while 1:
		c = sys.stdin.read(1)
		if chkvalid(c):
			break
	while 1:
		if chkvalid(c):
			_str += c
		else:
			break
		c = sys.stdin.read(1)
	return _str

# table of syntax conversions
def subs(string):
	a = ["ln", "\^"]
	b = ["log", "**"]
	for i in range(len(a)):
		string = re.sub(a[i], b[i], string)
	return string

matrix_n = [0 for i in range(500000)]
matrix_d = [0 for i in range(500000)]

print "Welcome to the deluxe matrix reducer"
print "Please separate entries by TABS !"
print "Expressions like 1 + 2 * 3 or ln(pi^2) + sin(3) are allowed"
print "=========================================="

cols = input("number of columns: ")
rows = input("number of rows: ")

print "matrix: "

k = 0
for i in range(rows):
	for j in range(cols):
		datum = eval(subs(symbol()))
		matrix_n[k] = n(datum)
		matrix_d[k] = d(datum)
		k += 1

print " "

mess_up = 0

# Column-by-column...
for c in range(cols):
	# We "clean up" this column row-by-row
	for r in range(rows):
		index = r * cols + c
		op = 0

		# 1. When r = c, we want a 1
		if r == c:
			# We only do work if the cell isn't
			# already a 1
			if not (matrix_d[index] != 0 and (matrix_n[index] == matrix_d[index])):
				# The way we get a 1 is that we divide this row by the
				# cell's value. We must make sure the cell is
				# a real nonzero number.
				if matrix_n[index] != 0 and matrix_d[index] != 0:
					# Say what we're going to do
					print "row "+str(r+1)+" divided by "+str(matrix_n[index] * 1.0 /matrix_d[index])+""

					# Do it
					tn = matrix_n[index]
					td = matrix_d[index]
					for c2 in range(cols):
						matrix_n[r * cols + c2] = int(matrix_n[r * cols + c2]) * int(td)
						matrix_d[r * cols + c2] = int(matrix_d[r * cols + c2]) * int(tn)

					# Print the new matrix
					op = 1

				# Otherwise
				if matrix_n[index] * matrix_d[index] == 0:
					mess_up = 1

		# 2. When r != c, we want a zero
		if r != c:
			# Needn't do anything more to this row
			# if it's already zero
			if matrix_n[index] != 0:
				# Look for a non-zero row
				# other than this one
				# to cancel this one off with
				chosen = -1
				for r2 in range(rows):
					different = (r2 != r)
					nonzero = (matrix_n[r2 * cols + c] != 0)
					defined = (matrix_d[r2 * cols + c] != 0)

					if different and nonzero and defined:
						# Count zeroes in the candidate row
						# in columns prior to this one
						tz = 0
						for c2 in range(c):
							if matrix_n[r2 * cols + c2] == 0 and (matrix_d[r2 * cols + c2] != 0):
								tz += 1

						# We choose the highest-found (lowest index)
						# row with trailing zeroes in all the
						# previous columns.
						if (tz == c) and (chosen == -1):
							chosen = r2

				if chosen == -1:
					mess_up = 1

				# Found a cancelling row
				if chosen >= 0:
					# This row - ratio * the other row
					# gives zero for this column of this row
					
					ratio_n = int(int(matrix_n[index]) * int(matrix_d[chosen * cols + c]))
					ratio_d = int(int(matrix_d[index]) * int(matrix_n[chosen * cols + c]))

					# Sign-normalize so things print out nicely
					if ((ratio_n > 0) and (ratio_d < 0)) or ((ratio_n < 0) and (ratio_d < 0)):
						ratio_n *= -1
						ratio_d *= -1
					
					# Do a quick simplification, again for niceness
					if (ratio_n * ratio_d != 0) and (ratio_n % ratio_d == 0):
						ratio_n /= ratio_d
						ratio_d = 1

					# Finally show the row-op
					sys.stdout.write("row "+str(r+1)+"");
					if ratio_n > 0:
						sys.stdout.write(" - ")
					if ratio_n < 0:
						sys.stdout.write(" + ")
					sys.stdout.write(""+str(abs(ratio_n) * 1.0 / abs(ratio_d))+"")
					sys.stdout.write(" * ")
					print " row "+str(chosen + 1)+""

					# Do the subtraction row-op
					for c2 in range(cols):
						# plain algebra serving linear algebra, eh
						#
						#	 a     c     e	   a	 ce     adf - bce	
						#	--- - --- * --- = --- - ---- = -----------
						#	 b     d     f	   b	 df 	   bdf
						
						aa = int(matrix_n[r * cols + c2])
						ab = int(matrix_d[r * cols + c2])
						ac = int(ratio_n)
						ad = int(ratio_d)
						ae = int(matrix_n[chosen * cols + c2])
						af = int(matrix_d[chosen * cols + c2])
						
						# come on python, don't use fucking floats !!!!
						new_num   = int(int(aa) * int(ad) * int(af)) - int(int(ab) * int(ac) * int(ae))
						new_denom = int(int(ab) * int(ad) * int(af))
						
						matrix_n[r * cols + c2] = int(new_num)
						matrix_d[r * cols + c2] = int(new_denom)
					op = 1

		if op:
			# If an operation was carried out, 
			# print out the new, operated-open
			# matrix, reducing fractions along
			# the way for the sake of prettiness
			# and integer-overflow prevention.
			print "--> "
			print ""
			for i in range(rows):
				sys.stdout.write(" | ")
				for j in range(cols):
					# Fraction reduction
					an = matrix_n[i * cols + j]
					ad = matrix_d[i * cols + j]
					_min = abs(ad)
					if abs(an) < abs(ad):
						_min = abs(an)

					p = gcd(int(an), int(ad))
					an = int(int(an) / int(p))
					ad = int(int(ad) / int(p))

					if (an != 0) and (ad != 0) and (an % ad == 0):
						an /= ad
						ad = 1	

					# Sign normalization
					if ((an < 0) and (ad < 0)) or ((an > 0) and (ad < 0)):
						an *= -1
						ad *= -1

					# Store the reduced fraction in its cell
					matrix_n[i * cols + j] = int(an)
					matrix_d[i * cols + j] = int(ad)

					# Print the cell
					sys.stdout.write(str((an * 1.0) / ad))
					#sys.stdout.write(str(an) + " / " + str(ad))

					if j < cols:
						sys.stdout.write(" | ")
				print " "
			print " "