Skip to content

Instantly share code, notes, and snippets.

@vampjaz
Created March 19, 2018 18:38
Show Gist options
  • Save vampjaz/ec987fcb22beb301ad54013f1dfb7112 to your computer and use it in GitHub Desktop.
Save vampjaz/ec987fcb22beb301ad54013f1dfb7112 to your computer and use it in GitHub Desktop.
Python (2) solver for real roots of a polynomial
## 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('>> ')))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment