Python (2) solver for real roots of a polynomial

algebra.py

## python algebra module

import math

def func_in(func,var='x'): ## lowercase var to split on
	# returns a list of the coefficients
	func = func.replace('-','+-').lower().replace(' ','')
	lst2 = {}
	try:
		for term in func.split('+'):
			if term:
				if var in term:
					coef,exp = term.split(var)
					if not exp:
						exp = 1
					if coef == '':
						coef = 1
					elif coef == '-':
						coef = -1
					try:
						exp = exp.replace('^','')
					except:
						pass
					coef = float(coef)
					exp = int(exp)
				else:
					coef = float(term)
					exp = 0
				lst2[exp] = coef
	except:
		print "arithmetic error, check equation?"
		return []
	maxpow = max(lst2.keys()) + 1
	lst = [0]*maxpow
	for i in range(maxpow):
		lst[i] = lst2.get(i,0)
	lst.reverse()
	return lst

def func_out(lst):
	# returns a function from a list of coefficients
	func = ''
	exp = len(lst) - 1
	for i in lst:
		if i != 0:
			if int(i) == i:
				i = int(i) ## remove decimals if not needed
			if i < 0:
				func += ' - ' + str(abs(i))
			else:
				func += ' + ' + str(i)
			if exp == 1:
				func += 'x'
			elif exp > 1:
				func += 'x^'+str(exp)
		exp-=1
	if func[:3] == ' + ':
		func = func[3:]
	return func

def num_factors(n): ## borrowed from http://stackoverflow.com/a/6800214/2635662
	return set(reduce(list.__add__, ([i, n//i] for i in range(1, int(n**0.5) + 1) if n % i == 0)))

def possible_zeros(lst):
	## possible rational factors
	p = num_factors(abs(lst[-1]))
	q = num_factors(abs(lst[0]))
	out = []
	for i in p:
		for j in q:
			yield (float(i)/float(j))
			yield -(float(i)/float(j))

def synthetic_div(lst,n,pr=True): #pr: print out work
	## synthetically divide
	lst2 = [0]
	lst3 = []
	for i in range(len(lst)):
		lst3.append(lst[i] + lst2[i])
		lst2.append(lst3[i] * n)
	del lst2[-1]
	if pr:
		print str(n) + '|\t\t' + '\t'.join(str(i) for i in lst)
		print '-'*len(str(n)) + '\t\t' + '\t'.join(str(i) for i in lst2)
		print '\t\t' + '\t'.join(str(i) for i in lst3)
	return lst3

def find_ratio_factors(func):
	## putting it all together
	lst = func_in(func)
	fact = []
	pz = set(possible_zeros(lst))
	for i in pz:
		l = synthetic_div(lst,i,pr=False)
		if l[-1] == 0: ## divisible by i
			fact.append(i)
	return fact

def quadratic(a,b,c):
	return ((-b + math.sqrt(b**2 - 4*a*c)) / 2*a,
			(-b - math.sqrt(b**2 - 4*a*c)) / 2*a)

'''
def cubic(a,b,c,d):
	q = math.sqrt((2 * b**3 - 9*a*b*c + 27 * a**2 * d)**2 - 4*(b**2 - 3*a*c)**3)
	z = (0.5 * (q + 2 * b**2 - 9*a*b*c + 27 * a**2 * d)) ** (1/3.0) ## cube root
	return (-b / (3*a)) - () ##todo:finish
'''

def find_real_factors(func):
	## putting it all together
	lst = func_in(func)
	tmp = lst
	fact = []
	pz = find_ratio_factors(func)*3 ## three times for up to multiplicity 3
	for i in pz:
		l = synthetic_div(tmp,i,pr=False)
		if l[-1] == 0: ## divisible by i
			tmp = l[:-1]
			fact.append(i)
	if tmp != [1]: ## means x = 1, all the way factored
		## uh oh
		if len(tmp) == 3: ## quadratics can be solved...
			for i in quadratic(tmp[0],tmp[1],tmp[2]):
				fact.append(i)
	return fact

print 'factors: ' + ', '.join(find_real_factors(raw_input('>> ')))