added findvalues function and better handling with non polynomial expressions and irrational numbers

shiroindev.py

@@ -4,10 +4,10 @@ import warnings,operator
 warnings.filterwarnings("ignore")	
 #Seed is needed to select the weights in linprog function.
 #None means that the seed is random.
-class Empty:
+class Vars:
 	pass
-sVars=Empty()
-sVars.print=print
+sVars=Vars()
+sVars.display=print
 sVars.seed=None
 sVars.translation={}
 translationList=['numerator:','denominator:','status:',
@@ -24,9 +24,10 @@ translationList=['numerator:','denominator:','status:',
 #Initialize english-english dictionary.
 for phrase in translationList:
 	sVars.translation[phrase]=phrase
-from scipy.optimize import linprog
+from scipy.optimize import linprog,fmin
 import random
-from sympy import S,cancel,fraction,Pow,expand,solve,latex,oo
+from sympy import S,cancel,fraction,Pow,expand,solve,latex,oo,Poly,lambdify
+from collections import Counter
 import re
 def _remzero(coef,fun): 
 #coef, fun represents an expression.
@@ -72,6 +73,28 @@ def ssolve(formula,variables):
 	return result
 def sstr(formula):
 	return str(formula).replace('**','^').replace('*','').replace(' ','')
+def reducegens(formula):
+	pol=Poly(formula)
+	newgens={}
+	ind={}
+	for gen in pol.gens:
+		base,pw=_powr(gen)
+		coef,_=pw.as_coeff_mul()
+		ml=pw/coef
+		if base**ml in newgens:
+			newgens[base**ml]=gcd(newgens[base**ml],coef)
+		else:
+			newgens[base**ml]=coef
+			ind[base**ml]=S('tmp'+str(len(ind)))
+	for gen in pol.gens:
+		base,pw=_powr(gen)
+		coef,_=pw.as_coeff_mul()
+		ml=pw/coef
+		pol=pol.replace(gen,ind[base**ml]**(coef/newgens[base**ml]))
+	newpol=Poly(pol.as_expr())
+	for gen in newgens:
+		newpol=newpol.replace(ind[gen],gen**newgens[gen])
+	return newpol
 def Sm(formula):
 #Adds multiplication signs and sympifies a formula.
 #For example, Sm('(2x+y)(7+5xz)') -> S('(2*x+y)*(7+5*x*z)')
@@ -88,14 +111,14 @@ def _input2fraction(formula,variables,values):
 	subst=[]
 	for x,y in zip(variables,values):
 		if y!=1:
-			sVars.print(sVars.translation['Substitute']+' $'+str(x)+'\\to '+slatex(S(y)*S(x))+'$')
+			sVars.display(sVars.translation['Substitute']+' $'+str(x)+'\\to '+slatex(S(y)*S(x))+'$')
 			subst+=[(x,x*y)]
 	formula=formula.subs(subst)
 	numerator,denominator=fractioncancel(formula)
-	sVars.print(sVars.translation['numerator:']+' $'+slatex(numerator)+'$')
-	sVars.print(sVars.translation['denominator:']+' $'+slatex(denominator)+'$')
+	sVars.display(sVars.translation['numerator:']+' $'+slatex(numerator)+'$')
+	sVars.display(sVars.translation['denominator:']+' $'+slatex(denominator)+'$')
 	return (numerator,denominator)
-def _formula2list(formula,variables):
+def _formula2list(formula):
 #Splits a polynomial to a difference of two polynomials with positive
 #coefficients and extracts coefficients and powers of both polynomials.
 #'variables' is used to set order of powers
@@ -103,30 +126,13 @@ def _formula2list(formula,variables):
 #the program tries to prove that 
 #0<=5x^2-4xy+8y^3
 #4xy<=5x^2+8y^3
-#lcoef=[4]
-#lfun=[[1,1]]
-#rcoef=[5,8]
-#rfun=[[2,0],[0,3]]
-	lfun=[]
-	lcoef=[]
-	rfun=[]
-	rcoef=[]
-	varorder=dict(zip(variables,range(len(variables))))
-	for addend in formula.as_ordered_terms():
-		coef,facts=addend.as_coeff_mul()
-		powers=[0]*len(variables)
-		for var in variables:
-			powers[varorder[var]]=0
-		for fact in facts:
-			var,pw=_powr(fact)
-			powers[varorder[var]]=int(pw)
-		if(coef<0):
-			lcoef+=[-coef]
-			lfun+=[powers]
-		else:
-			rcoef+=[coef]
-			rfun+=[powers]
-	return(lcoef,lfun,rcoef,rfun)
+#returns [4],[(1,1)], [5,8],[(2,0),(0,3)], (x,y)
+	formula=reducegens(formula)
+	neg=(formula.abs()-formula)/2
+	pos=(formula.abs()+formula)/2
+	neg=Poly(neg,gens=formula.gens)
+	pos=Poly(pos,gens=formula.gens)
+	return neg.coeffs(),neg.monoms(),pos.coeffs(),pos.monoms(),Poly(formula).gens
 def _list2proof(lcoef,lfun,rcoef,rfun,variables,itermax,linprogiter,_writ2=_writ2,theorem="From weighted AM-GM inequality:"):
 #Now the formula is splitted on two polynomials with positive coefficients.
 #we will call them LHS and RHS and our inequality to prove would
@@ -157,13 +163,15 @@ def _list2proof(lcoef,lfun,rcoef,rfun,variables,itermax,linprogiter,_writ2=_writ
 		#If LHS is empty, then break. 
 	localseed=sVars.seed
 	bufer=[]
+	lcoef,lfun=_remzero(lcoef,lfun)
+	rcoef,rfun=_remzero(rcoef,rfun)
 	itern=0
 	if len(lcoef)==0: #if LHS is empty
-		sVars.print(sVars.translation['status:']+' 0')
+		sVars.display(sVars.translation['status:']+' 0')
 		status=0
 	elif len(rcoef)==0: 
 	#if RHS is empty, but LHS is not
-		sVars.print(sVars.translation['status:']+' 2')
+		sVars.display(sVars.translation['status:']+' 2')
 		status=2
 		itermax=0
 	foundreal=0
@@ -200,9 +208,9 @@ def _list2proof(lcoef,lfun,rcoef,rfun,variables,itermax,linprogiter,_writ2=_writ
 		res=linprog(vecc,A_eq=A,b_eq=b,A_ub=A_ub,b_ub=b_ub,options={'maxiter':linprogiter})
 		status=res.status
 		if itern==1:
-			sVars.print(sVars.translation['status:']+' '+str(status))
+			sVars.display(sVars.translation['status:']+' '+str(status))
 			if status==0:
-				sVars.print(sVars.translation[theorem])
+				sVars.display(sVars.translation[theorem])
 		if status==2: #if real solution of current inequality doesn't exist
 			if foundreal==0: #if this is the first inequality, then break 
 				break
@@ -214,7 +222,7 @@ def _list2proof(lcoef,lfun,rcoef,rfun,variables,itermax,linprogiter,_writ2=_writ
 				continue
 		if status==0:#if found a solution with real coefficients
 			for ineq in bufer:
-				sVars.print(ineq)
+				sVars.display(ineq)
 			foundreal=1
 			bufer=[]
 			oldlfun,oldrfun=lfun,rfun
@@ -246,25 +254,25 @@ def _list2proof(lcoef,lfun,rcoef,rfun,variables,itermax,linprogiter,_writ2=_writ
 			lcoef,lfun=_remzero(lcoef,lfun)
 			rcoef,rfun=_remzero(rcoef,rfun)
 	for ineq in bufer:
-		sVars.print(ineq)
+		sVars.display(ineq)
 	lhs='+'.join([_writ2(c,f,variables) for c,f in zip(lcoef,lfun)])
 	if lhs=='':
 		lhs='0'
 	elif status==0:
-		sVars.print(sVars.translation[
+		sVars.display(sVars.translation[
 		"Program couldn't find a solution with integer coefficients. Try "+
 		"to multiple the formula by some integer and run this function again."])
 	elif(status==2):
-		sVars.print(sVars.translation["Program couldn't find any proof."])
+		sVars.display(sVars.translation["Program couldn't find any proof."])
 		#return res.status
 	elif status==1:
-		sVars.print(sVars.translation["Try to set higher linprogiter parameter."])
+		sVars.display(sVars.translation["Try to set higher linprogiter parameter."])
 	rhs='+'.join([_writ2(c,f,variables) for c,f in zip(rcoef,rfun)])
 	if rhs=='':
 		rhs='0'
-	sVars.print('$$ '+slatex(lhs)+' \\le '+slatex(rhs)+' $$')
+	sVars.display('$$ '+slatex(lhs)+' \\le '+slatex(rhs)+' $$')
 	if lhs=='0':
-		sVars.print(sVars.translation['The sum of all inequalities gives us a proof of the inequality.'])
+		sVars.display(sVars.translation['The sum of all inequalities gives us a proof of the inequality.'])
 	return status
 def _isiterable(obj):
 	try:
@@ -290,13 +298,13 @@ def prove(formula,values=None,variables=None,niter=200,linprogiter=10000):
 	if values: values=_smakeiterable(values)
 	else: values=[1]*len(variables)
 	num,den=_input2fraction(formula,variables,values)
-	st=_list2proof(*(_formula2list(num,variables)+(variables,niter,linprogiter)))
+	st=_list2proof(*(_formula2list(num)+(niter,linprogiter)))
 	if st==2 and issymetric(num):
 		fs=sorted([str(x) for x in num.free_symbols])
-		sVars.print(sVars.translation["It looks like the formula is symmetric. "+
+		sVars.display(sVars.translation["It looks like the formula is symmetric. "+
 		"You can assume without loss of generality that "]+
 		' >= '.join([str(x) for x in fs])+'. '+sVars.translation['Try'])
-		sVars.print('prove(makesubs(S("'+str(num)+'"),'+
+		sVars.display('prove(makesubs(S("'+str(num)+'"),'+
 		str([(str(x),'oo') for x in variables[1:]])+')')
 	return st
 def powerprove(formula,values=None,variables=None,niter=200,linprogiter=10000):
@@ -309,10 +317,11 @@ def powerprove(formula,values=None,variables=None,niter=200,linprogiter=10000):
 	else: values=[1]*len(variables)
 	num,den=_input2fraction(formula,variables,values)
 	subst2=[]
+	statusses=[]
 	for j in range(len(variables)):
 		subst2+=[(variables[j],1+variables[j])]
 	for i in range(1<<len(variables)): #tricky substitutions to improve speed
-		sVars.print('\n\\hline\n')
+		sVars.display('\n\\hline\n')
 		subst1=[]
 		substout=[]
 		for j in range(len(variables)):
@@ -321,11 +330,12 @@ def powerprove(formula,values=None,variables=None,niter=200,linprogiter=10000):
 				substout+=[str(variables[j])+'\\to 1/(1+'+str(variables[j])+')']
 			else:
 				substout+=[str(variables[j])+'\\to 1+'+str(variables[j])]
-		sVars.print(sVars.translation['Substitute']+ ' $'+','.join(substout)+'$')
+		sVars.display(sVars.translation['Substitute']+ ' $'+','.join(substout)+'$')
 		num1=fractioncancel(num.subs(subst1))[0]
 		num2=expand(num1.subs(subst2))
-		sVars.print(sVars.translation["Numerator after substitutions:"]+' $'+slatex(num2)+'$')
-		_list2proof(*(_formula2list(num2,variables)+(variables,niter,linprogiter)))
+		sVars.display(sVars.translation["Numerator after substitutions:"]+' $'+slatex(num2)+'$')
+		statusses+=[_list2proof(*(_formula2list(num2)+(niter,linprogiter)))]
+	return Counter(statusses)
 def makesubs(formula,intervals,values=None,variables=None,numden=False):
 #This function generates a new formula which satisfies this condition:
 #for all positive variables new formula is nonnegative iff
@@ -346,18 +356,18 @@ def makesubs(formula,intervals,values=None,variables=None,numden=False):
 		if {end1,end2}=={S('-oo'),S('oo')}:
 			formula=formula.subs(var,var-1/var)
 			equations=[equation.subs(var,var-1/var) for equation in equations]
-			sVars.print(sVars.translation['Substitute']+' $'+str(var)+'\\to '+var-1/var+'$')
+			sVars.display(sVars.translation['Substitute']+' $'+str(var)+'\\to '+sstr(var-1/var)+'$')
 		elif end2==S('oo'):
 			formula=formula.subs(var,end1+var)
 			equations=[equation.subs(var,end1+var) for equation in equations]
-			sVars.print(sVars.translation['Substitute']+' $'+str(var)+'\\to '+sstr(end1+var)+'$')
+			sVars.display(sVars.translation['Substitute']+' $'+str(var)+'\\to '+sstr(end1+var)+'$')
 		elif end2==S('-oo'):
 			formula=formula.subs(var,end1-var)
 			equations=[equation.subs(var,end1-var) for equation in equations]
-			sVars.print(sVars.translation['Substitute']+" $"+str(var)+'\\to '+sstr(end1-var)+'$')
+			sVars.display(sVars.translation['Substitute']+" $"+str(var)+'\\to '+sstr(end1-var)+'$')
 		else:
 			formula=formula.subs(var,end2+(end1-end2)/var)
-			sVars.print(sVars.translation['Substitute']+" $"+str(var)+'\\to '+sstr(end2+(end1-end2)/(1+var))+'$')
+			sVars.display(sVars.translation['Substitute']+" $"+str(var)+'\\to '+sstr(end2+(end1-end2)/(1+var))+'$')
 			equations=[equation.subs(var,end2+(end1-end2)/var) for equation in equations]
 	num,den=fractioncancel(formula)
 	for var,interval in zip(variables,intervals):
@@ -370,7 +380,7 @@ def makesubs(formula,intervals,values=None,variables=None,numden=False):
 		if len(values):
 			values=values[0]
 	num,den=expand(num),expand(den)
-	#sVars.print(sVars.translation["Formula after substitution:"],"$$",slatex(num/den),'$$')
+	#sVars.display(sVars.translation["Formula after substitution:"],"$$",slatex(num/den),'$$')
 	if values and numden:
 		return num,den,values
 	elif values:
@@ -395,7 +405,7 @@ def _formula2listf(formula):
 		else:
 			rcoef+=[coef]
 			rfun+=[facts[0].args]
-	return(lcoef,lfun,rcoef,rfun)
+	return(lcoef,lfun,rcoef,rfun,None)
 def provef(formula,niter=200,linprogiter=10000):
 #this function is similar to prove, formula is a linear combination of
 #values of f:R^k->R instead of a polynomial. provef checks if a formula
@@ -403,7 +413,7 @@ def provef(formula,niter=200,linprogiter=10000):
 #provides a proof of nonnegativity.
 	formula=S(formula)
 	num,den=_input2fraction(formula,[],[])
-	_list2proof(*(_formula2listf(num)+(None,niter,linprogiter,_writ,'From Jensen inequality:')))
+	return _list2proof(*(_formula2listf(num)+(niter,linprogiter,_writ,'From Jensen inequality:')))
 def issymetric(formula): #checks if formula is symmetric
 #and has at least two variables	
 	if len(formula.free_symbols)<2:
@@ -447,3 +457,19 @@ def symmetrize(formula,oper=operator.add,variables=None,init=None):
 	for i in range(1,len(variables)):
 		formula=cyclize(formula,oper,variables[:i+1])
 	return formula
+def findvalues(formula,values=None,variables=None):
+	formula=S(formula)
+	num,den=fractioncancel(formula)
+	if variables==None:
+		variables=sorted(num.free_symbols,key=str)
+	num=num.subs(zip(variables,list(map(lambda x:x**2,variables))))
+	num=Poly(num)
+	newformula=S((num.abs()+num)/(num.abs()-num))
+	f=lambdify(variables,newformula)
+	f2=lambda x:f(*x)
+	if values==None:
+		values=[1.0]*len(variables)
+	else:
+		values=S(values)
+	tup=tuple(fmin(f2,values))
+	return tuple([x*x for x in tup])