Linear equations Python ver

linear.py

#!/usr/bin/python2

class vec(object):
    def __init__(self, val):
        self.val = list(val)
    
    def __add__(self, other):
        return vec(a+b for a, b in zip(self.val, other.val))

    def __sub__(self, other):
        return vec(a-b for a, b, in zip(self.val, other.val))

    def __mul__(self, other):
        return vec(a*other for a in self.val)

    def __div__(self, other):
        return vec(a*1.0/other for a in self.val)

    def __getitem__(self, i):
        return self.val[i]

    def __setitem__(self, i, v):
        self.val[i] = v

    def __repr__(self):
        return 'vec(%s)' % repr(self.val)


print 'Input n of unknowns:',
n = int(raw_input())

mat = []

print 'Input coefficients:'

for i in xrange(n):
    mat.append(vec(
        float(raw_input()) for j in range(n+1)
    ))

# begin
for i in range(n):
    if not mat[i][i]: # if cur coef is 0
        for j in range(i+1, n):
            if mat[j][i]: # if some other is not
                mat[i], mat[j] = mat[j], mat[i] # swap them
                break
        else:
            continue # all the coefs are 0
    co = mat[i][i]
    mat[i] /= co
    for j in range(i+1, n):
        mat[j] -= mat[i] * mat[j][i]

for i in reversed(range(n)):
    val = mat[i].val
    if not any(val[:-1]): # coef all zero
        if val[-1]: # but b is not?!
            raise Exception('No solution!')
    else:
        rank = i + 1
        break

print rank

for i in range(1, rank):
    for j in range(i):
        mat[j] -= mat[i] * mat[j][i]

print '(%s)' %  ', '.join(str(mat[i][n]) for i in range(n)),

for i in range(n - rank):
    line = [-mat[j][rank+i] for j in range(rank)]
    v1 = [0] * (n-rank)
    v1[i] = 1
    line.extend(v1)
    print '+ c%d * (%s)' % (i+1, ', '.join(str(v) for v in line)),