shiroindev/.ipynb_checkpoints/sandbox-checkpoint.ipynb
2020-05-17 23:46:53 +02:00

276 KiB
Raw Blame History

def extrema(fr):
    t=[]
    for x in fr.free_symbols:
        t+=[fr.diff(x)]
    print(t)
    return solve(t)
from shiroindev import *
from sympy import *
from itertools import permutations, combinations
shiro.seed=1
from IPython.display import Latex
shiro.display=lambda x:display(Latex(x))
Latex(r'__________________________')
__________________________
prove(makesubs(Sm('-a(s-b)-b(s-c)-c(s-a)+ s^2'),'[0,s],[0,s],[0,s]'))
Substitute $a\to s-\frac{s}{d+1}$
Substitute $b\to s-\frac{s}{e+1}$
Substitute $c\to s-\frac{s}{f+1}$
numerator: $defs^2+s^2$
denominator: $def+de+df+d+ef+e+f+1$
status: 0
$$ 0 \le defs^2+s^2 $$
The sum of all inequalities gives us a proof of the inequality.
0
def nonneg(formula):
    formula=expand(formula)
    for addend in formula.as_ordered_terms():
        coef,facts=addend.as_coeff_mul()
        if coef<0:
            return False
    return True
def symprove(formula,n):
    formula=S(formula)
    if n==0:
        return
    ls=list(formula.free_symbols)
    for i in range(len(ls)):
        a=ls[i]
        for j in range(i+1,len(ls)):
            b=ls[j]
            if expand(formula-formula.subs({a:b, b:a}, simultaneous=True))==S(0):
                formula=makesubs(formula,[[b,S('oo')]],variables=[a,b])
                sVars.display('$$'+latex(formula)+'$$')
                symprove(formula,n-1)
symprove('x^2-2*x*y+y^2',4)
Substitute $x\to g+y$
$$g^{2}$$
def provesym(formula,n):
    formula=S(formula)
    if n==0:
        return
    fs=list(formula.free_symbols)
    print 
    for i in range(2,len(fs)+1):
        for fs2 in combinations(fs,i):
            for fsp in permutations(fs2[1:]):
                if expand(formula-formula.subs(zip((fs2[0],)+fsp,fsp+(fs2[0],)), simultaneous=True))==S(0):
                    newformula=makesubs(formula,[[fs2[0],oo]]*(len(fsp)),variables=fsp)
                    sVars.display(str(n)+' $$'+latex(newformula)+'$$')
                    provesym(newformula,n-1)
provesym(Sm('x^t(x-y)(x-z) + y^t (y-z)(y-x) + z^t (z-x)(z-y)'),4)
Substitute $y\to h+x$
4 $$h^{2} \left(h + x\right)^{t} - h x x^{t} + h x z^{t} + h x \left(h + x\right)^{t} + h x^{t} z - h z z^{t} - h z \left(h + x\right)^{t} + x^{2} z^{t} - 2 x z z^{t} + z^{2} z^{t}$$
Substitute $z\to i+x$
4 $$i^{2} \left(i + x\right)^{t} - i x x^{t} + i x y^{t} + i x \left(i + x\right)^{t} + i x^{t} y - i y y^{t} - i y \left(i + x\right)^{t} + x^{2} y^{t} - 2 x y y^{t} + y^{2} y^{t}$$
Substitute $z\to j+y$
4 $$j^{2} \left(j + y\right)^{t} - j x x^{t} + j x y^{t} - j x \left(j + y\right)^{t} + j x^{t} y - j y y^{t} + j y \left(j + y\right)^{t} + x^{2} x^{t} - 2 x x^{t} y + x^{t} y^{2}$$
Substitute $y\to k+x$
Substitute $z\to l+x$
4 $$k^{2} \left(k + x\right)^{t} + k l x^{t} - k l \left(k + x\right)^{t} - k l \left(l + x\right)^{t} + l^{2} \left(l + x\right)^{t}$$
Substitute $l\to k+m$
3 $$k^{2} x^{t} + k m x^{t} - k m \left(k + x\right)^{t} + k m \left(k + m + x\right)^{t} + m^{2} \left(k + m + x\right)^{t}$$
Substitute $z\to n+x$
Substitute $y\to o+x$
4 $$n^{2} \left(n + x\right)^{t} + n o x^{t} - n o \left(n + x\right)^{t} - n o \left(o + x\right)^{t} + o^{2} \left(o + x\right)^{t}$$
Substitute $o\to n+p$
3 $$n^{2} x^{t} + n p x^{t} - n p \left(n + x\right)^{t} + n p \left(n + p + x\right)^{t} + p^{2} \left(n + p + x\right)^{t}$$
makesubs(Sm('x^t(x-y)(x-z) + y^t (y-z)(y-x) + z^t (z-x)(z-y)'),'[y,oo]',variables='x')
Substitute $x\to q+y$
$\displaystyle q^{2} \left(q + y\right)^{t} - q y y^{t} + q y z^{t} + q y \left(q + y\right)^{t} + q y^{t} z - q z z^{t} - q z \left(q + y\right)^{t} + y^{2} z^{t} - 2 y z z^{t} + z^{2} z^{t}$
makesubs(Sm('x^t*(x-y)(x-z) + y^t*(y-z)(y-x) + z^t*(z-x)(z-y)'),'[y,oo]',variables='x')
Substitute $x\to r+y$
$\displaystyle r^{2} \left(r + y\right)^{t} - r y y^{t} + r y z^{t} + r y \left(r + y\right)^{t} + r y^{t} z - r z z^{t} - r z \left(r + y\right)^{t} + y^{2} z^{t} - 2 y z z^{t} + z^{2} z^{t}$
provesym(Sm('x^t(x-y)(x-z) + y^t (y-z)(y-x) + z^t (z-x)(z-y)'),4)
Substitute $y\to u+x$
4 $$u^{2} \left(u + x\right)^{t} - u x x^{t} + u x z^{t} + u x \left(u + x\right)^{t} + u x^{t} z - u z z^{t} - u z \left(u + x\right)^{t} + x^{2} z^{t} - 2 x z z^{t} + z^{2} z^{t}$$
Substitute $z\to v+x$
4 $$v^{2} \left(v + x\right)^{t} - v x x^{t} + v x y^{t} + v x \left(v + x\right)^{t} + v x^{t} y - v y y^{t} - v y \left(v + x\right)^{t} + x^{2} y^{t} - 2 x y y^{t} + y^{2} y^{t}$$
Substitute $z\to w+y$
4 $$w^{2} \left(w + y\right)^{t} - w x x^{t} + w x y^{t} - w x \left(w + y\right)^{t} + w x^{t} y - w y y^{t} + w y \left(w + y\right)^{t} + x^{2} x^{t} - 2 x x^{t} y + x^{t} y^{2}$$
Substitute $y\to a_{1}+x$
Substitute $z\to b_{1}+x$
4 $$a_{1}^{2} \left(a_{1} + x\right)^{t} + a_{1} b_{1} x^{t} - a_{1} b_{1} \left(a_{1} + x\right)^{t} - a_{1} b_{1} \left(b_{1} + x\right)^{t} + b_{1}^{2} \left(b_{1} + x\right)^{t}$$
Substitute $b_{1}\to a_{1}+c_{1}$
3 $$a_{1}^{2} x^{t} + a_{1} c_{1} x^{t} - a_{1} c_{1} \left(a_{1} + x\right)^{t} + a_{1} c_{1} \left(a_{1} + c_{1} + x\right)^{t} + c_{1}^{2} \left(a_{1} + c_{1} + x\right)^{t}$$
Substitute $z\to d_{1}+x$
Substitute $y\to e_{1}+x$
4 $$d_{1}^{2} \left(d_{1} + x\right)^{t} + d_{1} e_{1} x^{t} - d_{1} e_{1} \left(d_{1} + x\right)^{t} - d_{1} e_{1} \left(e_{1} + x\right)^{t} + e_{1}^{2} \left(e_{1} + x\right)^{t}$$
Substitute $e_{1}\to d_{1}+f_{1}$
3 $$d_{1}^{2} x^{t} + d_{1} f_{1} x^{t} - d_{1} f_{1} \left(d_{1} + x\right)^{t} + d_{1} f_{1} \left(d_{1} + f_{1} + x\right)^{t} + f_{1}^{2} \left(d_{1} + f_{1} + x\right)^{t}$$
Sm('x^t(x-y)(x-z) + y^t (y-z)(y-x) + z^t (z-x)(z-y)')
$\displaystyle x^{t} \left(x - y\right) \left(x - z\right) + y^{t} \left(- x + y\right) \left(y - z\right) + z^{t} \left(- x + z\right) \left(- y + z\right)$
S('x+2*y+3*z').subs(zip(S('[x,y,z]'),S('[y,z,x]')),simultaneous=True)
$\displaystyle 3 x + y + 2 z$
S('[x,y,z]')
[x, y, z]
formula=Sm('x^2(x+y)^t-xyy^t+xyz^t+xy(x+y)^t+xy^tz-xzzt-xz(x+y)^t+y^2z^t-2yzz^t+z^2z^t')
expand(formula.subs(S('[y,z],[z,y]'),simultaneous=True)-formula)
$\displaystyle - t x y^{2} + t x z^{2} - x^{2} \left(x + y\right)^{t} + x^{2} \left(x + z\right)^{t} + x y y^{t} - x y \left(x + y\right)^{t} - x y \left(x + z\right)^{t} - x z z^{t} + x z \left(x + y\right)^{t} + x z \left(x + z\right)^{t} + y^{2} y^{t} - y^{2} z^{t} - 2 y y^{t} z + 2 y z z^{t} + y^{t} z^{2} - z^{2} z^{t}$
def point(formula):
    fs=list(formula.free_symbols)
    il=1
    for s in fs[1:]:
        il*=s
    fr=formula.subs(fs[0],1/il)
    print(fr)
    return(extrema(fr))
#point(Sm('(a^2+b^2+c^2+d^2)-a*(b+c+d)'))
prove(Sm('(a^2+b^2+c^2+d^2)-a*(b+c+d)'),'11/8,4/5,4/5,4/5')
Substitute $a\to 11g_{1}/8$
Substitute $b\to 4h_{1}/5$
Substitute $c\to 4i_{1}/5$
Substitute $d\to 4j_{1}/5$
numerator: $3025g_{1}^2-1760g_{1}h_{1}-1760g_{1}i_{1}-1760g_{1}j_{1}+1024h_{1}^2+1024i_{1}^2+1024j_{1}^2$
denominator: $1600$
status: 0
From weighted AM-GM inequality:
$$1760g_{1}h_{1} \le 880g_{1}^2+880h_{1}^2$$
$$1760g_{1}i_{1} \le 880g_{1}^2+880i_{1}^2$$
$$1760g_{1}j_{1} \le 880g_{1}^2+880j_{1}^2$$
$$ 0 \le 385g_{1}^2+144h_{1}^2+144i_{1}^2+144j_{1}^2 $$
The sum of all inequalities gives us a proof of the inequality.
0
216/125-5/3
0.06133333333333324
"""formula=cyclize(Sm('((a-b)/c)^2-8**(1/2)*(a-b)/c'))
display(formula)
from scipy.optimize import fmin
import numpy as np
def f(x):
    num,den=fraction(cancel(newformula))
    fs=sorted(newformula.free_symbols,key=str)
    return float(num.subs(list(zip(fs,x))))
newformula=(makesubs(formula,'[b,oo],[c,oo]'))
print(f([2,1,1]))
fmin(f,[2,1,1])
print(fmin(f,np.array([2,1,1])))
display(simplify(newformula))"""
#prove(newformula)
#prove(makesubs(formula,'[b,oo],[a,oo]',variables='c,b,a'))
"formula=cyclize(Sm('((a-b)/c)^2-8**(1/2)*(a-b)/c'))\ndisplay(formula)\nfrom scipy.optimize import fmin\nimport numpy as np\ndef f(x):\n    num,den=fraction(cancel(newformula))\n    fs=sorted(newformula.free_symbols,key=str)\n    return float(num.subs(list(zip(fs,x))))\nnewformula=(makesubs(formula,'[b,oo],[c,oo]'))\nprint(f([2,1,1]))\nfmin(f,[2,1,1])\nprint(fmin(f,np.array([2,1,1])))\ndisplay(simplify(newformula))"
def g(x):
    return x[0]+x[1]
fmin(g,[0,0])
Warning: Maximum number of function evaluations has been exceeded.
array([-1.62831265e+41, -1.93382892e+41])
formula=Sm('(a^2+b^2+c^2+d^2)/(a*(b+c+d))')
display(formula)
def f(x):
    fs=sorted(formula.free_symbols,key=str)
    return formula.subs(zip(fs,x))
nsimplify(tuple(fmin(f,(1,1,1,1))),tolerance=0.1,rational=True)
$\displaystyle \frac{a^{2} + b^{2} + c^{2} + d^{2}}{a \left(b + c + d\right)}$
Optimization terminated successfully.
         Current function value: 1.154701
         Iterations: 80
         Function evaluations: 144
$\displaystyle \left( \frac{11}{8}, \ \frac{4}{5}, \ \frac{4}{5}, \ \frac{4}{5}\right)$
_[0]/_[1]
$\displaystyle \frac{55}{32}$
formula=cyclize(Sm('((a-b)/c)^2-8**(1/2)*(a-b)/c'))
formula=(makesubs(formula,'[b,oo],[c,oo]'))
num,den=fraction(cancel(formula))
num=num.subs(zip(fs,list(map(lambda x:x**2,fs))))
display(num)
numm=0
nump=0
for addend in num.as_ordered_terms():
    coef,facts=addend.as_coeff_mul()
    if coef<0:
        numm-=addend
    else:
        nump+=addend
num=nump/numm
fs=sorted(num.free_symbols,key=str)
numm,nump=Poly(numm),Poly(nump)
def f(x):
    return nump.eval(dict(zip(fs,x)))/numm.eval(dict(zip(fs,x)))  
print(dict(zip(fs,(2,2,2))))
display(f((2,2,2)))
print('x')
tup=tuple(fmin(f,(2,2,2)))
display(tuple([x*x for x in tup]))
nsimplify(tuple([x*x for x in tup]),tolerance=0.1,rational=True)
Substitute $a\to b+k_{1}$
Substitute $b\to c+l_{1}$
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-24-10599ede4c20> in <module>()
      2 formula=(makesubs(formula,'[b,oo],[c,oo]'))
      3 num,den=fraction(cancel(formula))
----> 4 num=num.subs(zip(fs,list(map(lambda x:x**2,fs))))
      5 display(num)
      6 numm=0

NameError: name 'fs' is not defined
nsimplify((2.1002573656763053, 10.340431462974655, 1.7661001788212371),tolerance=0.3,rational=True)
$\displaystyle \left( 2, \ \frac{31}{3}, \ \frac{7}{4}\right)$
(2.1002573656763053/1.7661001788212371)**2
1.414211498165453
formula=cyclize(Sm('((a-b)/c)^2-8**(1/2)*(a-b)/c'))
display(Latex('Case $a\ge c\ge b$'))
formula1=makesubs(formula,'[c,oo],[b,oo]',variables='a,c,b')
prove(formula1)
display(Latex('Case $a\ge b\ge c$'))
formula2=makesubs(formula,'[b,oo],[c,oo]')
prove(formula2*4,values='1,sqrt(2),1')
formula
Case $a\ge c\ge b$
Substitute $a\to c+r_{1}$
Substitute $c\to b+s_{1}$
numerator: $2b^4r_{1}^2+2b^4r_{1}s_{1}+2b^4s_{1}^2+4b^3r_{1}^3+2\sqrt{2}b^3r_{1}^2s_{1}+10b^3r_{1}^2s_{1}+2\sqrt{2}b^3r_{1}s_{1}^2+6b^3r_{1}s_{1}^2+4b^3s_{1}^3+2b^2r_{1}^4+2\sqrt{2}b^2r_{1}^3s_{1}+10b^2r_{1}^3s_{1}+6\sqrt{2}b^2r_{1}^2s_{1}^2+12b^2r_{1}^2s_{1}^2+4b^2r_{1}s_{1}^3+4\sqrt{2}b^2r_{1}s_{1}^3+2b^2s_{1}^4+2br_{1}^4s_{1}+2\sqrt{2}br_{1}^3s_{1}^2+6br_{1}^3s_{1}^2+4br_{1}^2s_{1}^3+4\sqrt{2}br_{1}^2s_{1}^3+2\sqrt{2}br_{1}s_{1}^4+r_{1}^4s_{1}^2+2r_{1}^3s_{1}^3+r_{1}^2s_{1}^4$
denominator: $b^6+2b^5r_{1}+4b^5s_{1}+b^4r_{1}^2+6b^4r_{1}s_{1}+6b^4s_{1}^2+2b^3r_{1}^2s_{1}+6b^3r_{1}s_{1}^2+4b^3s_{1}^3+b^2r_{1}^2s_{1}^2+2b^2r_{1}s_{1}^3+b^2s_{1}^4$
status: 0
$$ 0 \le r_{1}^4s_{1}^2+2br_{1}^4s_{1}+2b^2r_{1}^4+2r_{1}^3s_{1}^3+2\sqrt{2}br_{1}^3s_{1}^2+6br_{1}^3s_{1}^2+2\sqrt{2}b^2r_{1}^3s_{1}+10b^2r_{1}^3s_{1}+4b^3r_{1}^3+r_{1}^2s_{1}^4+4\sqrt{2}br_{1}^2s_{1}^3+4br_{1}^2s_{1}^3+6\sqrt{2}b^2r_{1}^2s_{1}^2+12b^2r_{1}^2s_{1}^2+2\sqrt{2}b^3r_{1}^2s_{1}+10b^3r_{1}^2s_{1}+2b^4r_{1}^2+2\sqrt{2}br_{1}s_{1}^4+4\sqrt{2}b^2r_{1}s_{1}^3+4b^2r_{1}s_{1}^3+2\sqrt{2}b^3r_{1}s_{1}^2+6b^3r_{1}s_{1}^2+2b^4r_{1}s_{1}+2b^2s_{1}^4+4b^3s_{1}^3+2b^4s_{1}^2 $$
The sum of all inequalities gives us a proof of the inequality.
Case $a\ge b\ge c$
Substitute $a\to b+t_{1}$
Substitute $b\to c+u_{1}$
Substitute $t_{1}\to \sqrt{2}v_{1}$
numerator: $8c^4u_{1}^2+8\sqrt{2}c^4u_{1}v_{1}+16c^4v_{1}^2+16c^3u_{1}^3-16c^3u_{1}^2v_{1}+24\sqrt{2}c^3u_{1}^2v_{1}-16\sqrt{2}c^3u_{1}v_{1}^2+80c^3u_{1}v_{1}^2+32\sqrt{2}c^3v_{1}^3+8c^2u_{1}^4-32c^2u_{1}^3v_{1}+16\sqrt{2}c^2u_{1}^3v_{1}-48\sqrt{2}c^2u_{1}^2v_{1}^2+96c^2u_{1}^2v_{1}^2-32c^2u_{1}v_{1}^3+80\sqrt{2}c^2u_{1}v_{1}^3+32c^2v_{1}^4-16cu_{1}^4v_{1}-32\sqrt{2}cu_{1}^3v_{1}^2+32cu_{1}^3v_{1}^2-32cu_{1}^2v_{1}^3+48\sqrt{2}cu_{1}^2v_{1}^3+32cu_{1}v_{1}^4+8u_{1}^4v_{1}^2+16\sqrt{2}u_{1}^3v_{1}^3+16u_{1}^2v_{1}^4$
denominator: $c^6+4c^5u_{1}+2\sqrt{2}c^5v_{1}+6c^4u_{1}^2+6\sqrt{2}c^4u_{1}v_{1}+2c^4v_{1}^2+4c^3u_{1}^3+6\sqrt{2}c^3u_{1}^2v_{1}+4c^3u_{1}v_{1}^2+c^2u_{1}^4+2\sqrt{2}c^2u_{1}^3v_{1}+2c^2u_{1}^2v_{1}^2$
status: 0
From weighted AM-GM inequality:
$$16cu_{1}^4v_{1} \le 8u_{1}^4v_{1}^2+8c^2u_{1}^4$$
$$32\sqrt{2}cu_{1}^3v_{1}^2 \le 16\sqrt{2}u_{1}^3v_{1}^3+16\sqrt{2}c^2u_{1}^3v_{1}$$
$$32c^2u_{1}^3v_{1} \le 16cu_{1}^3v_{1}^2+16c^3u_{1}^3$$
$$32cu_{1}^2v_{1}^3 \le 16cu_{1}^3v_{1}^2+16cu_{1}v_{1}^4$$
$$48\sqrt{2}c^2u_{1}^2v_{1}^2 \le 24\sqrt{2}cu_{1}^2v_{1}^3+24\sqrt{2}c^3u_{1}^2v_{1}$$
$$16c^3u_{1}^2v_{1} \le 8c^2u_{1}^2v_{1}^2+8c^4u_{1}^2$$
$$32c^2u_{1}v_{1}^3 \le 16u_{1}^2v_{1}^4+16c^4v_{1}^2$$
$$16\sqrt{2}c^3u_{1}v_{1}^2 \le 8\sqrt{2}c^2u_{1}v_{1}^3+8\sqrt{2}c^4u_{1}v_{1}$$
$$ 0 \le 24\sqrt{2}cu_{1}^2v_{1}^3+88c^2u_{1}^2v_{1}^2+16cu_{1}v_{1}^4+72\sqrt{2}c^2u_{1}v_{1}^3+80c^3u_{1}v_{1}^2+32c^2v_{1}^4+32\sqrt{2}c^3v_{1}^3 $$
The sum of all inequalities gives us a proof of the inequality.
$\displaystyle - \frac{2 \sqrt{2} \left(a - b\right)}{c} + \frac{\left(a - b\right)^{2}}{c^{2}} - \frac{2 \sqrt{2} \left(- a + c\right)}{b} + \frac{\left(- a + c\right)^{2}}{b^{2}} - \frac{2 \sqrt{2} \left(b - c\right)}{a} + \frac{\left(b - c\right)^{2}}{a^{2}}$
formula=cyclize(Sm('((a-b)/c)^2-8**(1/2)*(a-b)/c'))
formula=(makesubs(formula,'[b,oo],[c,oo]'))
formula,_=fraction(cancel(formula))
formula=formula.subs('a','2**(1/2)*a')
sVars.display(r'Substitute $a\to\sqrt{2}a$')
prove(formula)
Substitute $a\to b+w_{1}$
Substitute $b\to c+x_{1}$
Substitute $a\to\sqrt{2}a$
numerator: $2c^4w_{1}^2+2c^4w_{1}x_{1}+2c^4x_{1}^2+4c^3w_{1}^3-2\sqrt{2}c^3w_{1}^2x_{1}+10c^3w_{1}^2x_{1}-2\sqrt{2}c^3w_{1}x_{1}^2+6c^3w_{1}x_{1}^2+4c^3x_{1}^3+2c^2w_{1}^4-2\sqrt{2}c^2w_{1}^3x_{1}+10c^2w_{1}^3x_{1}-6\sqrt{2}c^2w_{1}^2x_{1}^2+12c^2w_{1}^2x_{1}^2-4\sqrt{2}c^2w_{1}x_{1}^3+4c^2w_{1}x_{1}^3+2c^2x_{1}^4+2cw_{1}^4x_{1}-2\sqrt{2}cw_{1}^3x_{1}^2+6cw_{1}^3x_{1}^2-4\sqrt{2}cw_{1}^2x_{1}^3+4cw_{1}^2x_{1}^3-2\sqrt{2}cw_{1}x_{1}^4+w_{1}^4x_{1}^2+2w_{1}^3x_{1}^3+w_{1}^2x_{1}^4$
denominator: $1$
status: 2
Program couldn't find any proof.
$$ 2\sqrt{2}cw_{1}x_{1}^4+4\sqrt{2}cw_{1}^2x_{1}^3+4\sqrt{2}c^2w_{1}x_{1}^3+2\sqrt{2}cw_{1}^3x_{1}^2+6\sqrt{2}c^2w_{1}^2x_{1}^2+2\sqrt{2}c^3w_{1}x_{1}^2+2\sqrt{2}c^2w_{1}^3x_{1}+2\sqrt{2}c^3w_{1}^2x_{1} \le w_{1}^2x_{1}^4+2c^2x_{1}^4+2w_{1}^3x_{1}^3+4cw_{1}^2x_{1}^3+4c^2w_{1}x_{1}^3+4c^3x_{1}^3+w_{1}^4x_{1}^2+6cw_{1}^3x_{1}^2+12c^2w_{1}^2x_{1}^2+6c^3w_{1}x_{1}^2+2c^4x_{1}^2+2cw_{1}^4x_{1}+10c^2w_{1}^3x_{1}+10c^3w_{1}^2x_{1}+2c^4w_{1}x_{1}+2c^2w_{1}^4+4c^3w_{1}^3+2c^4w_{1}^2 $$
2
S('2**(1/2)*x*y').as_coeff_mul()
(1, (x, y, sqrt(2)))
S('sqrt(2)*x*y').as_coeff_mul()
(1, (x, y, sqrt(2)))
Poly('sqrt(2)*x*y+8').monoms()
[(1, 1, 1), (0, 0, 0)]
Poly('sqrt(2)*x*y+58').coeffs()
[1, 58]
Poly('sqrt(2)*x*y+sqrt(3)+sqrt(2)*sqrt(3)').gens
(x, y, sqrt(2), sqrt(3), sqrt(6))
def _formula2list(formula):
    neg=pos=0
    for addend in formula.as_ordered_terms():
        coef,facts=addend.as_coeff_mul()
        if coef<0:
            neg-=addend
        else:
            pos+=addend
    neg=Poly(neg,gens=Poly(formula).gens)
    pos=Poly(pos,gens=Poly(formula).gens)
    return neg.coeffs(),neg.monoms(),pos.coeffs(),pos.monoms(),Poly(formula).gens
_formula2list(S('x^2+7*y'))
([0], [(0, 0)], [1, 7], [(2, 0), (0, 1)], (x, y))
prove(S('sqrt(2)*x^2-sqrt(8)*x*y+sqrt(2)*y^2'))
numerator: $\sqrt{2}x^2-2\sqrt{2}xy+\sqrt{2}y^2$
denominator: $1$
status: 0
From weighted AM-GM inequality:
$$2\sqrt{2}xy \le \sqrt{2}x^2+\sqrt{2}y^2$$
$$ 0 \le 0 $$
The sum of all inequalities gives us a proof of the inequality.
0
[Poly('1+2**(1/3)+4**(1/3)+x').monoms(),Poly('1+2**(1/3)+4**(1/3)+x').coeffs()]
[[(1, 0), (0, 2), (0, 1), (0, 0)], [1, 1, 1, 1]]
x=Symbol('x', positive=True)
Poly(x+sqrt(x))
$\displaystyle \operatorname{Poly}{\left( x + \sqrt{x}, x, \sqrt{x}, domain=\mathbb{Z} \right)}$
prove?
Poly('x^2-1').abs()
$\displaystyle \operatorname{Poly}{\left( x^{2} + 1, x, domain=\mathbb{Z} \right)}$
Poly(Poly('x')/Poly('x+y'))
$\displaystyle \operatorname{Poly}{\left( x\frac{1}{x + y}, x, \frac{1}{x + y}, domain=\mathbb{Z} \right)}$
Poly('sqrt(2)+x+sqrt(x)+x**(1/4)+x**(1/3)+x**(2/3)',extension=1)
$\displaystyle \operatorname{Poly}{\left( x + \sqrt{x} + x^{\frac{2}{3}} + \sqrt[3]{x} + \sqrt[4]{x} + \sqrt{2}, x, \sqrt{x}, \sqrt[3]{x}, \sqrt[4]{x}, domain=QQ \right)}$
Poly('x+sqrt(6)',S('x'),S('sqrt(2)'),S('sqrt(3)')).coeffs()
[1, sqrt(6)]
Poly('x+sqrt(x)',S('sqrt(x)')).coeffs()
[1, x]
source(Poly)
/home/grzegorz/Pobrane/SageMath/local/lib/python3.7/site-packages/sympy/core/decorators.py:38: SymPyDeprecationWarning: 

source has been deprecated since SymPy 1.3. Use ?? in IPython/Jupyter
or inspect.getsource instead. See
https://github.com/sympy/sympy/issues/14905 for more info.

  _warn_deprecation(wrapped, 3)
In file: /home/grzegorz/Pobrane/SageMath/local/lib/python3.7/site-packages/sympy/polys/polytools.py
class Poly(Expr):
    """
    Generic class for representing and operating on polynomial expressions.
    Subclasses Expr class.

    Examples
    ========

    >>> from sympy import Poly
    >>> from sympy.abc import x, y

    Create a univariate polynomial:

    >>> Poly(x*(x**2 + x - 1)**2)
    Poly(x**5 + 2*x**4 - x**3 - 2*x**2 + x, x, domain='ZZ')

    Create a univariate polynomial with specific domain:

    >>> from sympy import sqrt
    >>> Poly(x**2 + 2*x + sqrt(3), domain='R')
    Poly(1.0*x**2 + 2.0*x + 1.73205080756888, x, domain='RR')

    Create a multivariate polynomial:

    >>> Poly(y*x**2 + x*y + 1)
    Poly(x**2*y + x*y + 1, x, y, domain='ZZ')

    Create a univariate polynomial, where y is a constant:

    >>> Poly(y*x**2 + x*y + 1,x)
    Poly(y*x**2 + y*x + 1, x, domain='ZZ[y]')

    You can evaluate the above polynomial as a function of y:

    >>> Poly(y*x**2 + x*y + 1,x).eval(2)
    6*y + 1

    See Also
    ========

    sympy.core.expr.Expr

    """

    __slots__ = ['rep', 'gens']

    is_commutative = True
    is_Poly = True
    _op_priority = 10.001

    def __new__(cls, rep, *gens, **args):
        """Create a new polynomial instance out of something useful. """
        opt = options.build_options(gens, args)

        if 'order' in opt:
            raise NotImplementedError("'order' keyword is not implemented yet")

        if iterable(rep, exclude=str):
            if isinstance(rep, dict):
                return cls._from_dict(rep, opt)
            else:
                return cls._from_list(list(rep), opt)
        else:
            rep = sympify(rep)

            if rep.is_Poly:
                return cls._from_poly(rep, opt)
            else:
                return cls._from_expr(rep, opt)

    @classmethod
    def new(cls, rep, *gens):
        """Construct :class:`Poly` instance from raw representation. """
        if not isinstance(rep, DMP):
            raise PolynomialError(
                "invalid polynomial representation: %s" % rep)
        elif rep.lev != len(gens) - 1:
            raise PolynomialError("invalid arguments: %s, %s" % (rep, gens))

        obj = Basic.__new__(cls)

        obj.rep = rep
        obj.gens = gens

        return obj

    @classmethod
    def from_dict(cls, rep, *gens, **args):
        """Construct a polynomial from a ``dict``. """
        opt = options.build_options(gens, args)
        return cls._from_dict(rep, opt)

    @classmethod
    def from_list(cls, rep, *gens, **args):
        """Construct a polynomial from a ``list``. """
        opt = options.build_options(gens, args)
        return cls._from_list(rep, opt)

    @classmethod
    def from_poly(cls, rep, *gens, **args):
        """Construct a polynomial from a polynomial. """
        opt = options.build_options(gens, args)
        return cls._from_poly(rep, opt)

    @classmethod
    def from_expr(cls, rep, *gens, **args):
        """Construct a polynomial from an expression. """
        opt = options.build_options(gens, args)
        return cls._from_expr(rep, opt)

    @classmethod
    def _from_dict(cls, rep, opt):
        """Construct a polynomial from a ``dict``. """
        gens = opt.gens

        if not gens:
            raise GeneratorsNeeded(
                "can't initialize from 'dict' without generators")

        level = len(gens) - 1
        domain = opt.domain

        if domain is None:
            domain, rep = construct_domain(rep, opt=opt)
        else:
            for monom, coeff in rep.items():
                rep[monom] = domain.convert(coeff)

        return cls.new(DMP.from_dict(rep, level, domain), *gens)

    @classmethod
    def _from_list(cls, rep, opt):
        """Construct a polynomial from a ``list``. """
        gens = opt.gens

        if not gens:
            raise GeneratorsNeeded(
                "can't initialize from 'list' without generators")
        elif len(gens) != 1:
            raise MultivariatePolynomialError(
                "'list' representation not supported")

        level = len(gens) - 1
        domain = opt.domain

        if domain is None:
            domain, rep = construct_domain(rep, opt=opt)
        else:
            rep = list(map(domain.convert, rep))

        return cls.new(DMP.from_list(rep, level, domain), *gens)

    @classmethod
    def _from_poly(cls, rep, opt):
        """Construct a polynomial from a polynomial. """
        if cls != rep.__class__:
            rep = cls.new(rep.rep, *rep.gens)

        gens = opt.gens
        field = opt.field
        domain = opt.domain

        if gens and rep.gens != gens:
            if set(rep.gens) != set(gens):
                return cls._from_expr(rep.as_expr(), opt)
            else:
                rep = rep.reorder(*gens)

        if 'domain' in opt and domain:
            rep = rep.set_domain(domain)
        elif field is True:
            rep = rep.to_field()

        return rep

    @classmethod
    def _from_expr(cls, rep, opt):
        """Construct a polynomial from an expression. """
        rep, opt = _dict_from_expr(rep, opt)
        return cls._from_dict(rep, opt)

    def _hashable_content(self):
        """Allow SymPy to hash Poly instances. """
        return (self.rep, self.gens)

    def __hash__(self):
        return super(Poly, self).__hash__()

    @property
    def free_symbols(self):
        """
        Free symbols of a polynomial expression.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y, z

        >>> Poly(x**2 + 1).free_symbols
        {x}
        >>> Poly(x**2 + y).free_symbols
        {x, y}
        >>> Poly(x**2 + y, x).free_symbols
        {x, y}
        >>> Poly(x**2 + y, x, z).free_symbols
        {x, y}

        """
        symbols = set()
        gens = self.gens
        for i in range(len(gens)):
            for monom in self.monoms():
                if monom[i]:
                    symbols |= gens[i].free_symbols
                    break

        return symbols | self.free_symbols_in_domain

    @property
    def free_symbols_in_domain(self):
        """
        Free symbols of the domain of ``self``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x**2 + 1).free_symbols_in_domain
        set()
        >>> Poly(x**2 + y).free_symbols_in_domain
        set()
        >>> Poly(x**2 + y, x).free_symbols_in_domain
        {y}

        """
        domain, symbols = self.rep.dom, set()

        if domain.is_Composite:
            for gen in domain.symbols:
                symbols |= gen.free_symbols
        elif domain.is_EX:
            for coeff in self.coeffs():
                symbols |= coeff.free_symbols

        return symbols

    @property
    def args(self):
        """
        Don't mess up with the core.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 + 1, x).args
        (x**2 + 1,)

        """
        return (self.as_expr(),)

    @property
    def gen(self):
        """
        Return the principal generator.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 + 1, x).gen
        x

        """
        return self.gens[0]

    @property
    def domain(self):
        """Get the ground domain of ``self``. """
        return self.get_domain()

    @property
    def zero(self):
        """Return zero polynomial with ``self``'s properties. """
        return self.new(self.rep.zero(self.rep.lev, self.rep.dom), *self.gens)

    @property
    def one(self):
        """Return one polynomial with ``self``'s properties. """
        return self.new(self.rep.one(self.rep.lev, self.rep.dom), *self.gens)

    @property
    def unit(self):
        """Return unit polynomial with ``self``'s properties. """
        return self.new(self.rep.unit(self.rep.lev, self.rep.dom), *self.gens)

    def unify(f, g):
        """
        Make ``f`` and ``g`` belong to the same domain.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> f, g = Poly(x/2 + 1), Poly(2*x + 1)

        >>> f
        Poly(1/2*x + 1, x, domain='QQ')
        >>> g
        Poly(2*x + 1, x, domain='ZZ')

        >>> F, G = f.unify(g)

        >>> F
        Poly(1/2*x + 1, x, domain='QQ')
        >>> G
        Poly(2*x + 1, x, domain='QQ')

        """
        _, per, F, G = f._unify(g)
        return per(F), per(G)

    def _unify(f, g):
        g = sympify(g)

        if not g.is_Poly:
            try:
                return f.rep.dom, f.per, f.rep, f.rep.per(f.rep.dom.from_sympy(g))
            except CoercionFailed:
                raise UnificationFailed("can't unify %s with %s" % (f, g))

        if isinstance(f.rep, DMP) and isinstance(g.rep, DMP):
            gens = _unify_gens(f.gens, g.gens)

            dom, lev = f.rep.dom.unify(g.rep.dom, gens), len(gens) - 1

            if f.gens != gens:
                f_monoms, f_coeffs = _dict_reorder(
                    f.rep.to_dict(), f.gens, gens)

                if f.rep.dom != dom:
                    f_coeffs = [dom.convert(c, f.rep.dom) for c in f_coeffs]

                F = DMP(dict(list(zip(f_monoms, f_coeffs))), dom, lev)
            else:
                F = f.rep.convert(dom)

            if g.gens != gens:
                g_monoms, g_coeffs = _dict_reorder(
                    g.rep.to_dict(), g.gens, gens)

                if g.rep.dom != dom:
                    g_coeffs = [dom.convert(c, g.rep.dom) for c in g_coeffs]

                G = DMP(dict(list(zip(g_monoms, g_coeffs))), dom, lev)
            else:
                G = g.rep.convert(dom)
        else:
            raise UnificationFailed("can't unify %s with %s" % (f, g))

        cls = f.__class__

        def per(rep, dom=dom, gens=gens, remove=None):
            if remove is not None:
                gens = gens[:remove] + gens[remove + 1:]

                if not gens:
                    return dom.to_sympy(rep)

            return cls.new(rep, *gens)

        return dom, per, F, G

    def per(f, rep, gens=None, remove=None):
        """
        Create a Poly out of the given representation.

        Examples
        ========

        >>> from sympy import Poly, ZZ
        >>> from sympy.abc import x, y

        >>> from sympy.polys.polyclasses import DMP

        >>> a = Poly(x**2 + 1)

        >>> a.per(DMP([ZZ(1), ZZ(1)], ZZ), gens=[y])
        Poly(y + 1, y, domain='ZZ')

        """
        if gens is None:
            gens = f.gens

        if remove is not None:
            gens = gens[:remove] + gens[remove + 1:]

            if not gens:
                return f.rep.dom.to_sympy(rep)

        return f.__class__.new(rep, *gens)

    def set_domain(f, domain):
        """Set the ground domain of ``f``. """
        opt = options.build_options(f.gens, {'domain': domain})
        return f.per(f.rep.convert(opt.domain))

    def get_domain(f):
        """Get the ground domain of ``f``. """
        return f.rep.dom

    def set_modulus(f, modulus):
        """
        Set the modulus of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(5*x**2 + 2*x - 1, x).set_modulus(2)
        Poly(x**2 + 1, x, modulus=2)

        """
        modulus = options.Modulus.preprocess(modulus)
        return f.set_domain(FF(modulus))

    def get_modulus(f):
        """
        Get the modulus of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 + 1, modulus=2).get_modulus()
        2

        """
        domain = f.get_domain()

        if domain.is_FiniteField:
            return Integer(domain.characteristic())
        else:
            raise PolynomialError("not a polynomial over a Galois field")

    def _eval_subs(f, old, new):
        """Internal implementation of :func:`subs`. """
        if old in f.gens:
            if new.is_number:
                return f.eval(old, new)
            else:
                try:
                    return f.replace(old, new)
                except PolynomialError:
                    pass

        return f.as_expr().subs(old, new)

    def exclude(f):
        """
        Remove unnecessary generators from ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import a, b, c, d, x

        >>> Poly(a + x, a, b, c, d, x).exclude()
        Poly(a + x, a, x, domain='ZZ')

        """
        J, new = f.rep.exclude()
        gens = []

        for j in range(len(f.gens)):
            if j not in J:
                gens.append(f.gens[j])

        return f.per(new, gens=gens)

    def replace(f, x, y=None, *_ignore):
        # XXX this does not match Basic's signature
        """
        Replace ``x`` with ``y`` in generators list.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x**2 + 1, x).replace(x, y)
        Poly(y**2 + 1, y, domain='ZZ')

        """
        if y is None:
            if f.is_univariate:
                x, y = f.gen, x
            else:
                raise PolynomialError(
                    "syntax supported only in univariate case")

        if x == y or x not in f.gens:
            return f

        if x in f.gens and y not in f.gens:
            dom = f.get_domain()

            if not dom.is_Composite or y not in dom.symbols:
                gens = list(f.gens)
                gens[gens.index(x)] = y
                return f.per(f.rep, gens=gens)

        raise PolynomialError("can't replace %s with %s in %s" % (x, y, f))

    def reorder(f, *gens, **args):
        """
        Efficiently apply new order of generators.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x**2 + x*y**2, x, y).reorder(y, x)
        Poly(y**2*x + x**2, y, x, domain='ZZ')

        """
        opt = options.Options((), args)

        if not gens:
            gens = _sort_gens(f.gens, opt=opt)
        elif set(f.gens) != set(gens):
            raise PolynomialError(
                "generators list can differ only up to order of elements")

        rep = dict(list(zip(*_dict_reorder(f.rep.to_dict(), f.gens, gens))))

        return f.per(DMP(rep, f.rep.dom, len(gens) - 1), gens=gens)

    def ltrim(f, gen):
        """
        Remove dummy generators from ``f`` that are to the left of
        specified ``gen`` in the generators as ordered. When ``gen``
        is an integer, it refers to the generator located at that
        position within the tuple of generators of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y, z

        >>> Poly(y**2 + y*z**2, x, y, z).ltrim(y)
        Poly(y**2 + y*z**2, y, z, domain='ZZ')
        >>> Poly(z, x, y, z).ltrim(-1)
        Poly(z, z, domain='ZZ')

        """
        rep = f.as_dict(native=True)
        j = f._gen_to_level(gen)

        terms = {}

        for monom, coeff in rep.items():

            if any(i for i in monom[:j]):
                # some generator is used in the portion to be trimmed
                raise PolynomialError("can't left trim %s" % f)

            terms[monom[j:]] = coeff

        gens = f.gens[j:]

        return f.new(DMP.from_dict(terms, len(gens) - 1, f.rep.dom), *gens)

    def has_only_gens(f, *gens):
        """
        Return ``True`` if ``Poly(f, *gens)`` retains ground domain.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y, z

        >>> Poly(x*y + 1, x, y, z).has_only_gens(x, y)
        True
        >>> Poly(x*y + z, x, y, z).has_only_gens(x, y)
        False

        """
        indices = set()

        for gen in gens:
            try:
                index = f.gens.index(gen)
            except ValueError:
                raise GeneratorsError(
                    "%s doesn't have %s as generator" % (f, gen))
            else:
                indices.add(index)

        for monom in f.monoms():
            for i, elt in enumerate(monom):
                if i not in indices and elt:
                    return False

        return True

    def to_ring(f):
        """
        Make the ground domain a ring.

        Examples
        ========

        >>> from sympy import Poly, QQ
        >>> from sympy.abc import x

        >>> Poly(x**2 + 1, domain=QQ).to_ring()
        Poly(x**2 + 1, x, domain='ZZ')

        """
        if hasattr(f.rep, 'to_ring'):
            result = f.rep.to_ring()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'to_ring')

        return f.per(result)

    def to_field(f):
        """
        Make the ground domain a field.

        Examples
        ========

        >>> from sympy import Poly, ZZ
        >>> from sympy.abc import x

        >>> Poly(x**2 + 1, x, domain=ZZ).to_field()
        Poly(x**2 + 1, x, domain='QQ')

        """
        if hasattr(f.rep, 'to_field'):
            result = f.rep.to_field()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'to_field')

        return f.per(result)

    def to_exact(f):
        """
        Make the ground domain exact.

        Examples
        ========

        >>> from sympy import Poly, RR
        >>> from sympy.abc import x

        >>> Poly(x**2 + 1.0, x, domain=RR).to_exact()
        Poly(x**2 + 1, x, domain='QQ')

        """
        if hasattr(f.rep, 'to_exact'):
            result = f.rep.to_exact()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'to_exact')

        return f.per(result)

    def retract(f, field=None):
        """
        Recalculate the ground domain of a polynomial.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> f = Poly(x**2 + 1, x, domain='QQ[y]')
        >>> f
        Poly(x**2 + 1, x, domain='QQ[y]')

        >>> f.retract()
        Poly(x**2 + 1, x, domain='ZZ')
        >>> f.retract(field=True)
        Poly(x**2 + 1, x, domain='QQ')

        """
        dom, rep = construct_domain(f.as_dict(zero=True),
            field=field, composite=f.domain.is_Composite or None)
        return f.from_dict(rep, f.gens, domain=dom)

    def slice(f, x, m, n=None):
        """Take a continuous subsequence of terms of ``f``. """
        if n is None:
            j, m, n = 0, x, m
        else:
            j = f._gen_to_level(x)

        m, n = int(m), int(n)

        if hasattr(f.rep, 'slice'):
            result = f.rep.slice(m, n, j)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'slice')

        return f.per(result)

    def coeffs(f, order=None):
        """
        Returns all non-zero coefficients from ``f`` in lex order.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**3 + 2*x + 3, x).coeffs()
        [1, 2, 3]

        See Also
        ========
        all_coeffs
        coeff_monomial
        nth

        """
        return [f.rep.dom.to_sympy(c) for c in f.rep.coeffs(order=order)]

    def monoms(f, order=None):
        """
        Returns all non-zero monomials from ``f`` in lex order.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x**2 + 2*x*y**2 + x*y + 3*y, x, y).monoms()
        [(2, 0), (1, 2), (1, 1), (0, 1)]

        See Also
        ========
        all_monoms

        """
        return f.rep.monoms(order=order)

    def terms(f, order=None):
        """
        Returns all non-zero terms from ``f`` in lex order.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x**2 + 2*x*y**2 + x*y + 3*y, x, y).terms()
        [((2, 0), 1), ((1, 2), 2), ((1, 1), 1), ((0, 1), 3)]

        See Also
        ========
        all_terms

        """
        return [(m, f.rep.dom.to_sympy(c)) for m, c in f.rep.terms(order=order)]

    def all_coeffs(f):
        """
        Returns all coefficients from a univariate polynomial ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**3 + 2*x - 1, x).all_coeffs()
        [1, 0, 2, -1]

        """
        return [f.rep.dom.to_sympy(c) for c in f.rep.all_coeffs()]

    def all_monoms(f):
        """
        Returns all monomials from a univariate polynomial ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**3 + 2*x - 1, x).all_monoms()
        [(3,), (2,), (1,), (0,)]

        See Also
        ========
        all_terms

        """
        return f.rep.all_monoms()

    def all_terms(f):
        """
        Returns all terms from a univariate polynomial ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**3 + 2*x - 1, x).all_terms()
        [((3,), 1), ((2,), 0), ((1,), 2), ((0,), -1)]

        """
        return [(m, f.rep.dom.to_sympy(c)) for m, c in f.rep.all_terms()]

    def termwise(f, func, *gens, **args):
        """
        Apply a function to all terms of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> def func(k, coeff):
        ...     k = k[0]
        ...     return coeff//10**(2-k)

        >>> Poly(x**2 + 20*x + 400).termwise(func)
        Poly(x**2 + 2*x + 4, x, domain='ZZ')

        """
        terms = {}

        for monom, coeff in f.terms():
            result = func(monom, coeff)

            if isinstance(result, tuple):
                monom, coeff = result
            else:
                coeff = result

            if coeff:
                if monom not in terms:
                    terms[monom] = coeff
                else:
                    raise PolynomialError(
                        "%s monomial was generated twice" % monom)

        return f.from_dict(terms, *(gens or f.gens), **args)

    def length(f):
        """
        Returns the number of non-zero terms in ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 + 2*x - 1).length()
        3

        """
        return len(f.as_dict())

    def as_dict(f, native=False, zero=False):
        """
        Switch to a ``dict`` representation.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x**2 + 2*x*y**2 - y, x, y).as_dict()
        {(0, 1): -1, (1, 2): 2, (2, 0): 1}

        """
        if native:
            return f.rep.to_dict(zero=zero)
        else:
            return f.rep.to_sympy_dict(zero=zero)

    def as_list(f, native=False):
        """Switch to a ``list`` representation. """
        if native:
            return f.rep.to_list()
        else:
            return f.rep.to_sympy_list()

    def as_expr(f, *gens):
        """
        Convert a Poly instance to an Expr instance.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> f = Poly(x**2 + 2*x*y**2 - y, x, y)

        >>> f.as_expr()
        x**2 + 2*x*y**2 - y
        >>> f.as_expr({x: 5})
        10*y**2 - y + 25
        >>> f.as_expr(5, 6)
        379

        """
        if not gens:
            gens = f.gens
        elif len(gens) == 1 and isinstance(gens[0], dict):
            mapping = gens[0]
            gens = list(f.gens)

            for gen, value in mapping.items():
                try:
                    index = gens.index(gen)
                except ValueError:
                    raise GeneratorsError(
                        "%s doesn't have %s as generator" % (f, gen))
                else:
                    gens[index] = value

        return basic_from_dict(f.rep.to_sympy_dict(), *gens)

    def lift(f):
        """
        Convert algebraic coefficients to rationals.

        Examples
        ========

        >>> from sympy import Poly, I
        >>> from sympy.abc import x

        >>> Poly(x**2 + I*x + 1, x, extension=I).lift()
        Poly(x**4 + 3*x**2 + 1, x, domain='QQ')

        """
        if hasattr(f.rep, 'lift'):
            result = f.rep.lift()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'lift')

        return f.per(result)

    def deflate(f):
        """
        Reduce degree of ``f`` by mapping ``x_i**m`` to ``y_i``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x**6*y**2 + x**3 + 1, x, y).deflate()
        ((3, 2), Poly(x**2*y + x + 1, x, y, domain='ZZ'))

        """
        if hasattr(f.rep, 'deflate'):
            J, result = f.rep.deflate()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'deflate')

        return J, f.per(result)

    def inject(f, front=False):
        """
        Inject ground domain generators into ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> f = Poly(x**2*y + x*y**3 + x*y + 1, x)

        >>> f.inject()
        Poly(x**2*y + x*y**3 + x*y + 1, x, y, domain='ZZ')
        >>> f.inject(front=True)
        Poly(y**3*x + y*x**2 + y*x + 1, y, x, domain='ZZ')

        """
        dom = f.rep.dom

        if dom.is_Numerical:
            return f
        elif not dom.is_Poly:
            raise DomainError("can't inject generators over %s" % dom)

        if hasattr(f.rep, 'inject'):
            result = f.rep.inject(front=front)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'inject')

        if front:
            gens = dom.symbols + f.gens
        else:
            gens = f.gens + dom.symbols

        return f.new(result, *gens)

    def eject(f, *gens):
        """
        Eject selected generators into the ground domain.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> f = Poly(x**2*y + x*y**3 + x*y + 1, x, y)

        >>> f.eject(x)
        Poly(x*y**3 + (x**2 + x)*y + 1, y, domain='ZZ[x]')
        >>> f.eject(y)
        Poly(y*x**2 + (y**3 + y)*x + 1, x, domain='ZZ[y]')

        """
        dom = f.rep.dom

        if not dom.is_Numerical:
            raise DomainError("can't eject generators over %s" % dom)

        k = len(gens)

        if f.gens[:k] == gens:
            _gens, front = f.gens[k:], True
        elif f.gens[-k:] == gens:
            _gens, front = f.gens[:-k], False
        else:
            raise NotImplementedError(
                "can only eject front or back generators")

        dom = dom.inject(*gens)

        if hasattr(f.rep, 'eject'):
            result = f.rep.eject(dom, front=front)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'eject')

        return f.new(result, *_gens)

    def terms_gcd(f):
        """
        Remove GCD of terms from the polynomial ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x**6*y**2 + x**3*y, x, y).terms_gcd()
        ((3, 1), Poly(x**3*y + 1, x, y, domain='ZZ'))

        """
        if hasattr(f.rep, 'terms_gcd'):
            J, result = f.rep.terms_gcd()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'terms_gcd')

        return J, f.per(result)

    def add_ground(f, coeff):
        """
        Add an element of the ground domain to ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x + 1).add_ground(2)
        Poly(x + 3, x, domain='ZZ')

        """
        if hasattr(f.rep, 'add_ground'):
            result = f.rep.add_ground(coeff)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'add_ground')

        return f.per(result)

    def sub_ground(f, coeff):
        """
        Subtract an element of the ground domain from ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x + 1).sub_ground(2)
        Poly(x - 1, x, domain='ZZ')

        """
        if hasattr(f.rep, 'sub_ground'):
            result = f.rep.sub_ground(coeff)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'sub_ground')

        return f.per(result)

    def mul_ground(f, coeff):
        """
        Multiply ``f`` by a an element of the ground domain.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x + 1).mul_ground(2)
        Poly(2*x + 2, x, domain='ZZ')

        """
        if hasattr(f.rep, 'mul_ground'):
            result = f.rep.mul_ground(coeff)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'mul_ground')

        return f.per(result)

    def quo_ground(f, coeff):
        """
        Quotient of ``f`` by a an element of the ground domain.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(2*x + 4).quo_ground(2)
        Poly(x + 2, x, domain='ZZ')

        >>> Poly(2*x + 3).quo_ground(2)
        Poly(x + 1, x, domain='ZZ')

        """
        if hasattr(f.rep, 'quo_ground'):
            result = f.rep.quo_ground(coeff)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'quo_ground')

        return f.per(result)

    def exquo_ground(f, coeff):
        """
        Exact quotient of ``f`` by a an element of the ground domain.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(2*x + 4).exquo_ground(2)
        Poly(x + 2, x, domain='ZZ')

        >>> Poly(2*x + 3).exquo_ground(2)
        Traceback (most recent call last):
        ...
        ExactQuotientFailed: 2 does not divide 3 in ZZ

        """
        if hasattr(f.rep, 'exquo_ground'):
            result = f.rep.exquo_ground(coeff)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'exquo_ground')

        return f.per(result)

    def abs(f):
        """
        Make all coefficients in ``f`` positive.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 - 1, x).abs()
        Poly(x**2 + 1, x, domain='ZZ')

        """
        if hasattr(f.rep, 'abs'):
            result = f.rep.abs()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'abs')

        return f.per(result)

    def neg(f):
        """
        Negate all coefficients in ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 - 1, x).neg()
        Poly(-x**2 + 1, x, domain='ZZ')

        >>> -Poly(x**2 - 1, x)
        Poly(-x**2 + 1, x, domain='ZZ')

        """
        if hasattr(f.rep, 'neg'):
            result = f.rep.neg()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'neg')

        return f.per(result)

    def add(f, g):
        """
        Add two polynomials ``f`` and ``g``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 + 1, x).add(Poly(x - 2, x))
        Poly(x**2 + x - 1, x, domain='ZZ')

        >>> Poly(x**2 + 1, x) + Poly(x - 2, x)
        Poly(x**2 + x - 1, x, domain='ZZ')

        """
        g = sympify(g)

        if not g.is_Poly:
            return f.add_ground(g)

        _, per, F, G = f._unify(g)

        if hasattr(f.rep, 'add'):
            result = F.add(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'add')

        return per(result)

    def sub(f, g):
        """
        Subtract two polynomials ``f`` and ``g``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 + 1, x).sub(Poly(x - 2, x))
        Poly(x**2 - x + 3, x, domain='ZZ')

        >>> Poly(x**2 + 1, x) - Poly(x - 2, x)
        Poly(x**2 - x + 3, x, domain='ZZ')

        """
        g = sympify(g)

        if not g.is_Poly:
            return f.sub_ground(g)

        _, per, F, G = f._unify(g)

        if hasattr(f.rep, 'sub'):
            result = F.sub(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'sub')

        return per(result)

    def mul(f, g):
        """
        Multiply two polynomials ``f`` and ``g``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 + 1, x).mul(Poly(x - 2, x))
        Poly(x**3 - 2*x**2 + x - 2, x, domain='ZZ')

        >>> Poly(x**2 + 1, x)*Poly(x - 2, x)
        Poly(x**3 - 2*x**2 + x - 2, x, domain='ZZ')

        """
        g = sympify(g)

        if not g.is_Poly:
            return f.mul_ground(g)

        _, per, F, G = f._unify(g)

        if hasattr(f.rep, 'mul'):
            result = F.mul(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'mul')

        return per(result)

    def sqr(f):
        """
        Square a polynomial ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x - 2, x).sqr()
        Poly(x**2 - 4*x + 4, x, domain='ZZ')

        >>> Poly(x - 2, x)**2
        Poly(x**2 - 4*x + 4, x, domain='ZZ')

        """
        if hasattr(f.rep, 'sqr'):
            result = f.rep.sqr()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'sqr')

        return f.per(result)

    def pow(f, n):
        """
        Raise ``f`` to a non-negative power ``n``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x - 2, x).pow(3)
        Poly(x**3 - 6*x**2 + 12*x - 8, x, domain='ZZ')

        >>> Poly(x - 2, x)**3
        Poly(x**3 - 6*x**2 + 12*x - 8, x, domain='ZZ')

        """
        n = int(n)

        if hasattr(f.rep, 'pow'):
            result = f.rep.pow(n)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'pow')

        return f.per(result)

    def pdiv(f, g):
        """
        Polynomial pseudo-division of ``f`` by ``g``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 + 1, x).pdiv(Poly(2*x - 4, x))
        (Poly(2*x + 4, x, domain='ZZ'), Poly(20, x, domain='ZZ'))

        """
        _, per, F, G = f._unify(g)

        if hasattr(f.rep, 'pdiv'):
            q, r = F.pdiv(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'pdiv')

        return per(q), per(r)

    def prem(f, g):
        """
        Polynomial pseudo-remainder of ``f`` by ``g``.

        Caveat: The function prem(f, g, x) can be safely used to compute
          in Z[x] _only_ subresultant polynomial remainder sequences (prs's).

          To safely compute Euclidean and Sturmian prs's in Z[x]
          employ anyone of the corresponding functions found in
          the module sympy.polys.subresultants_qq_zz. The functions
          in the module with suffix _pg compute prs's in Z[x] employing
          rem(f, g, x), whereas the functions with suffix _amv
          compute prs's in Z[x] employing rem_z(f, g, x).

          The function rem_z(f, g, x) differs from prem(f, g, x) in that
          to compute the remainder polynomials in Z[x] it premultiplies
          the divident times the absolute value of the leading coefficient
          of the divisor raised to the power degree(f, x) - degree(g, x) + 1.


        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 + 1, x).prem(Poly(2*x - 4, x))
        Poly(20, x, domain='ZZ')

        """
        _, per, F, G = f._unify(g)

        if hasattr(f.rep, 'prem'):
            result = F.prem(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'prem')

        return per(result)

    def pquo(f, g):
        """
        Polynomial pseudo-quotient of ``f`` by ``g``.

        See the Caveat note in the function prem(f, g).

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 + 1, x).pquo(Poly(2*x - 4, x))
        Poly(2*x + 4, x, domain='ZZ')

        >>> Poly(x**2 - 1, x).pquo(Poly(2*x - 2, x))
        Poly(2*x + 2, x, domain='ZZ')

        """
        _, per, F, G = f._unify(g)

        if hasattr(f.rep, 'pquo'):
            result = F.pquo(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'pquo')

        return per(result)

    def pexquo(f, g):
        """
        Polynomial exact pseudo-quotient of ``f`` by ``g``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 - 1, x).pexquo(Poly(2*x - 2, x))
        Poly(2*x + 2, x, domain='ZZ')

        >>> Poly(x**2 + 1, x).pexquo(Poly(2*x - 4, x))
        Traceback (most recent call last):
        ...
        ExactQuotientFailed: 2*x - 4 does not divide x**2 + 1

        """
        _, per, F, G = f._unify(g)

        if hasattr(f.rep, 'pexquo'):
            try:
                result = F.pexquo(G)
            except ExactQuotientFailed as exc:
                raise exc.new(f.as_expr(), g.as_expr())
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'pexquo')

        return per(result)

    def div(f, g, auto=True):
        """
        Polynomial division with remainder of ``f`` by ``g``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 + 1, x).div(Poly(2*x - 4, x))
        (Poly(1/2*x + 1, x, domain='QQ'), Poly(5, x, domain='QQ'))

        >>> Poly(x**2 + 1, x).div(Poly(2*x - 4, x), auto=False)
        (Poly(0, x, domain='ZZ'), Poly(x**2 + 1, x, domain='ZZ'))

        """
        dom, per, F, G = f._unify(g)
        retract = False

        if auto and dom.is_Ring and not dom.is_Field:
            F, G = F.to_field(), G.to_field()
            retract = True

        if hasattr(f.rep, 'div'):
            q, r = F.div(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'div')

        if retract:
            try:
                Q, R = q.to_ring(), r.to_ring()
            except CoercionFailed:
                pass
            else:
                q, r = Q, R

        return per(q), per(r)

    def rem(f, g, auto=True):
        """
        Computes the polynomial remainder of ``f`` by ``g``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 + 1, x).rem(Poly(2*x - 4, x))
        Poly(5, x, domain='ZZ')

        >>> Poly(x**2 + 1, x).rem(Poly(2*x - 4, x), auto=False)
        Poly(x**2 + 1, x, domain='ZZ')

        """
        dom, per, F, G = f._unify(g)
        retract = False

        if auto and dom.is_Ring and not dom.is_Field:
            F, G = F.to_field(), G.to_field()
            retract = True

        if hasattr(f.rep, 'rem'):
            r = F.rem(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'rem')

        if retract:
            try:
                r = r.to_ring()
            except CoercionFailed:
                pass

        return per(r)

    def quo(f, g, auto=True):
        """
        Computes polynomial quotient of ``f`` by ``g``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 + 1, x).quo(Poly(2*x - 4, x))
        Poly(1/2*x + 1, x, domain='QQ')

        >>> Poly(x**2 - 1, x).quo(Poly(x - 1, x))
        Poly(x + 1, x, domain='ZZ')

        """
        dom, per, F, G = f._unify(g)
        retract = False

        if auto and dom.is_Ring and not dom.is_Field:
            F, G = F.to_field(), G.to_field()
            retract = True

        if hasattr(f.rep, 'quo'):
            q = F.quo(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'quo')

        if retract:
            try:
                q = q.to_ring()
            except CoercionFailed:
                pass

        return per(q)

    def exquo(f, g, auto=True):
        """
        Computes polynomial exact quotient of ``f`` by ``g``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 - 1, x).exquo(Poly(x - 1, x))
        Poly(x + 1, x, domain='ZZ')

        >>> Poly(x**2 + 1, x).exquo(Poly(2*x - 4, x))
        Traceback (most recent call last):
        ...
        ExactQuotientFailed: 2*x - 4 does not divide x**2 + 1

        """
        dom, per, F, G = f._unify(g)
        retract = False

        if auto and dom.is_Ring and not dom.is_Field:
            F, G = F.to_field(), G.to_field()
            retract = True

        if hasattr(f.rep, 'exquo'):
            try:
                q = F.exquo(G)
            except ExactQuotientFailed as exc:
                raise exc.new(f.as_expr(), g.as_expr())
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'exquo')

        if retract:
            try:
                q = q.to_ring()
            except CoercionFailed:
                pass

        return per(q)

    def _gen_to_level(f, gen):
        """Returns level associated with the given generator. """
        if isinstance(gen, int):
            length = len(f.gens)

            if -length <= gen < length:
                if gen < 0:
                    return length + gen
                else:
                    return gen
            else:
                raise PolynomialError("-%s <= gen < %s expected, got %s" %
                                      (length, length, gen))
        else:
            try:
                return f.gens.index(sympify(gen))
            except ValueError:
                raise PolynomialError(
                    "a valid generator expected, got %s" % gen)

    def degree(f, gen=0):
        """
        Returns degree of ``f`` in ``x_j``.

        The degree of 0 is negative infinity.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x**2 + y*x + 1, x, y).degree()
        2
        >>> Poly(x**2 + y*x + y, x, y).degree(y)
        1
        >>> Poly(0, x).degree()
        -oo

        """
        j = f._gen_to_level(gen)

        if hasattr(f.rep, 'degree'):
            return f.rep.degree(j)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'degree')

    def degree_list(f):
        """
        Returns a list of degrees of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x**2 + y*x + 1, x, y).degree_list()
        (2, 1)

        """
        if hasattr(f.rep, 'degree_list'):
            return f.rep.degree_list()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'degree_list')

    def total_degree(f):
        """
        Returns the total degree of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x**2 + y*x + 1, x, y).total_degree()
        2
        >>> Poly(x + y**5, x, y).total_degree()
        5

        """
        if hasattr(f.rep, 'total_degree'):
            return f.rep.total_degree()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'total_degree')

    def homogenize(f, s):
        """
        Returns the homogeneous polynomial of ``f``.

        A homogeneous polynomial is a polynomial whose all monomials with
        non-zero coefficients have the same total degree. If you only
        want to check if a polynomial is homogeneous, then use
        :func:`Poly.is_homogeneous`. If you want not only to check if a
        polynomial is homogeneous but also compute its homogeneous order,
        then use :func:`Poly.homogeneous_order`.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y, z

        >>> f = Poly(x**5 + 2*x**2*y**2 + 9*x*y**3)
        >>> f.homogenize(z)
        Poly(x**5 + 2*x**2*y**2*z + 9*x*y**3*z, x, y, z, domain='ZZ')

        """
        if not isinstance(s, Symbol):
            raise TypeError("``Symbol`` expected, got %s" % type(s))
        if s in f.gens:
            i = f.gens.index(s)
            gens = f.gens
        else:
            i = len(f.gens)
            gens = f.gens + (s,)
        if hasattr(f.rep, 'homogenize'):
            return f.per(f.rep.homogenize(i), gens=gens)
        raise OperationNotSupported(f, 'homogeneous_order')

    def homogeneous_order(f):
        """
        Returns the homogeneous order of ``f``.

        A homogeneous polynomial is a polynomial whose all monomials with
        non-zero coefficients have the same total degree. This degree is
        the homogeneous order of ``f``. If you only want to check if a
        polynomial is homogeneous, then use :func:`Poly.is_homogeneous`.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> f = Poly(x**5 + 2*x**3*y**2 + 9*x*y**4)
        >>> f.homogeneous_order()
        5

        """
        if hasattr(f.rep, 'homogeneous_order'):
            return f.rep.homogeneous_order()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'homogeneous_order')

    def LC(f, order=None):
        """
        Returns the leading coefficient of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(4*x**3 + 2*x**2 + 3*x, x).LC()
        4

        """
        if order is not None:
            return f.coeffs(order)[0]

        if hasattr(f.rep, 'LC'):
            result = f.rep.LC()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'LC')

        return f.rep.dom.to_sympy(result)

    def TC(f):
        """
        Returns the trailing coefficient of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**3 + 2*x**2 + 3*x, x).TC()
        0

        """
        if hasattr(f.rep, 'TC'):
            result = f.rep.TC()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'TC')

        return f.rep.dom.to_sympy(result)

    def EC(f, order=None):
        """
        Returns the last non-zero coefficient of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**3 + 2*x**2 + 3*x, x).EC()
        3

        """
        if hasattr(f.rep, 'coeffs'):
            return f.coeffs(order)[-1]
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'EC')

    def coeff_monomial(f, monom):
        """
        Returns the coefficient of ``monom`` in ``f`` if there, else None.

        Examples
        ========

        >>> from sympy import Poly, exp
        >>> from sympy.abc import x, y

        >>> p = Poly(24*x*y*exp(8) + 23*x, x, y)

        >>> p.coeff_monomial(x)
        23
        >>> p.coeff_monomial(y)
        0
        >>> p.coeff_monomial(x*y)
        24*exp(8)

        Note that ``Expr.coeff()`` behaves differently, collecting terms
        if possible; the Poly must be converted to an Expr to use that
        method, however:

        >>> p.as_expr().coeff(x)
        24*y*exp(8) + 23
        >>> p.as_expr().coeff(y)
        24*x*exp(8)
        >>> p.as_expr().coeff(x*y)
        24*exp(8)

        See Also
        ========
        nth: more efficient query using exponents of the monomial's generators

        """
        return f.nth(*Monomial(monom, f.gens).exponents)

    def nth(f, *N):
        """
        Returns the ``n``-th coefficient of ``f`` where ``N`` are the
        exponents of the generators in the term of interest.

        Examples
        ========

        >>> from sympy import Poly, sqrt
        >>> from sympy.abc import x, y

        >>> Poly(x**3 + 2*x**2 + 3*x, x).nth(2)
        2
        >>> Poly(x**3 + 2*x*y**2 + y**2, x, y).nth(1, 2)
        2
        >>> Poly(4*sqrt(x)*y)
        Poly(4*y*(sqrt(x)), y, sqrt(x), domain='ZZ')
        >>> _.nth(1, 1)
        4

        See Also
        ========
        coeff_monomial

        """
        if hasattr(f.rep, 'nth'):
            if len(N) != len(f.gens):
                raise ValueError('exponent of each generator must be specified')
            result = f.rep.nth(*list(map(int, N)))
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'nth')

        return f.rep.dom.to_sympy(result)

    def coeff(f, x, n=1, right=False):
        # the semantics of coeff_monomial and Expr.coeff are different;
        # if someone is working with a Poly, they should be aware of the
        # differences and chose the method best suited for the query.
        # Alternatively, a pure-polys method could be written here but
        # at this time the ``right`` keyword would be ignored because Poly
        # doesn't work with non-commutatives.
        raise NotImplementedError(
            'Either convert to Expr with `as_expr` method '
            'to use Expr\'s coeff method or else use the '
            '`coeff_monomial` method of Polys.')

    def LM(f, order=None):
        """
        Returns the leading monomial of ``f``.

        The Leading monomial signifies the monomial having
        the highest power of the principal generator in the
        expression f.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(4*x**2 + 2*x*y**2 + x*y + 3*y, x, y).LM()
        x**2*y**0

        """
        return Monomial(f.monoms(order)[0], f.gens)

    def EM(f, order=None):
        """
        Returns the last non-zero monomial of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(4*x**2 + 2*x*y**2 + x*y + 3*y, x, y).EM()
        x**0*y**1

        """
        return Monomial(f.monoms(order)[-1], f.gens)

    def LT(f, order=None):
        """
        Returns the leading term of ``f``.

        The Leading term signifies the term having
        the highest power of the principal generator in the
        expression f along with its coefficient.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(4*x**2 + 2*x*y**2 + x*y + 3*y, x, y).LT()
        (x**2*y**0, 4)

        """
        monom, coeff = f.terms(order)[0]
        return Monomial(monom, f.gens), coeff

    def ET(f, order=None):
        """
        Returns the last non-zero term of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(4*x**2 + 2*x*y**2 + x*y + 3*y, x, y).ET()
        (x**0*y**1, 3)

        """
        monom, coeff = f.terms(order)[-1]
        return Monomial(monom, f.gens), coeff

    def max_norm(f):
        """
        Returns maximum norm of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(-x**2 + 2*x - 3, x).max_norm()
        3

        """
        if hasattr(f.rep, 'max_norm'):
            result = f.rep.max_norm()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'max_norm')

        return f.rep.dom.to_sympy(result)

    def l1_norm(f):
        """
        Returns l1 norm of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(-x**2 + 2*x - 3, x).l1_norm()
        6

        """
        if hasattr(f.rep, 'l1_norm'):
            result = f.rep.l1_norm()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'l1_norm')

        return f.rep.dom.to_sympy(result)

    def clear_denoms(self, convert=False):
        """
        Clear denominators, but keep the ground domain.

        Examples
        ========

        >>> from sympy import Poly, S, QQ
        >>> from sympy.abc import x

        >>> f = Poly(x/2 + S(1)/3, x, domain=QQ)

        >>> f.clear_denoms()
        (6, Poly(3*x + 2, x, domain='QQ'))
        >>> f.clear_denoms(convert=True)
        (6, Poly(3*x + 2, x, domain='ZZ'))

        """
        f = self

        if not f.rep.dom.is_Field:
            return S.One, f

        dom = f.get_domain()
        if dom.has_assoc_Ring:
            dom = f.rep.dom.get_ring()

        if hasattr(f.rep, 'clear_denoms'):
            coeff, result = f.rep.clear_denoms()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'clear_denoms')

        coeff, f = dom.to_sympy(coeff), f.per(result)

        if not convert or not dom.has_assoc_Ring:
            return coeff, f
        else:
            return coeff, f.to_ring()

    def rat_clear_denoms(self, g):
        """
        Clear denominators in a rational function ``f/g``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> f = Poly(x**2/y + 1, x)
        >>> g = Poly(x**3 + y, x)

        >>> p, q = f.rat_clear_denoms(g)

        >>> p
        Poly(x**2 + y, x, domain='ZZ[y]')
        >>> q
        Poly(y*x**3 + y**2, x, domain='ZZ[y]')

        """
        f = self

        dom, per, f, g = f._unify(g)

        f = per(f)
        g = per(g)

        if not (dom.is_Field and dom.has_assoc_Ring):
            return f, g

        a, f = f.clear_denoms(convert=True)
        b, g = g.clear_denoms(convert=True)

        f = f.mul_ground(b)
        g = g.mul_ground(a)

        return f, g

    def integrate(self, *specs, **args):
        """
        Computes indefinite integral of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x**2 + 2*x + 1, x).integrate()
        Poly(1/3*x**3 + x**2 + x, x, domain='QQ')

        >>> Poly(x*y**2 + x, x, y).integrate((0, 1), (1, 0))
        Poly(1/2*x**2*y**2 + 1/2*x**2, x, y, domain='QQ')

        """
        f = self

        if args.get('auto', True) and f.rep.dom.is_Ring:
            f = f.to_field()

        if hasattr(f.rep, 'integrate'):
            if not specs:
                return f.per(f.rep.integrate(m=1))

            rep = f.rep

            for spec in specs:
                if type(spec) is tuple:
                    gen, m = spec
                else:
                    gen, m = spec, 1

                rep = rep.integrate(int(m), f._gen_to_level(gen))

            return f.per(rep)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'integrate')

    def diff(f, *specs, **kwargs):
        """
        Computes partial derivative of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x**2 + 2*x + 1, x).diff()
        Poly(2*x + 2, x, domain='ZZ')

        >>> Poly(x*y**2 + x, x, y).diff((0, 0), (1, 1))
        Poly(2*x*y, x, y, domain='ZZ')

        """
        if not kwargs.get('evaluate', True):
            return Derivative(f, *specs, **kwargs)

        if hasattr(f.rep, 'diff'):
            if not specs:
                return f.per(f.rep.diff(m=1))

            rep = f.rep

            for spec in specs:
                if type(spec) is tuple:
                    gen, m = spec
                else:
                    gen, m = spec, 1

                rep = rep.diff(int(m), f._gen_to_level(gen))

            return f.per(rep)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'diff')

    _eval_derivative = diff

    def eval(self, x, a=None, auto=True):
        """
        Evaluate ``f`` at ``a`` in the given variable.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y, z

        >>> Poly(x**2 + 2*x + 3, x).eval(2)
        11

        >>> Poly(2*x*y + 3*x + y + 2, x, y).eval(x, 2)
        Poly(5*y + 8, y, domain='ZZ')

        >>> f = Poly(2*x*y + 3*x + y + 2*z, x, y, z)

        >>> f.eval({x: 2})
        Poly(5*y + 2*z + 6, y, z, domain='ZZ')
        >>> f.eval({x: 2, y: 5})
        Poly(2*z + 31, z, domain='ZZ')
        >>> f.eval({x: 2, y: 5, z: 7})
        45

        >>> f.eval((2, 5))
        Poly(2*z + 31, z, domain='ZZ')
        >>> f(2, 5)
        Poly(2*z + 31, z, domain='ZZ')

        """
        f = self

        if a is None:
            if isinstance(x, dict):
                mapping = x

                for gen, value in mapping.items():
                    f = f.eval(gen, value)

                return f
            elif isinstance(x, (tuple, list)):
                values = x

                if len(values) > len(f.gens):
                    raise ValueError("too many values provided")

                for gen, value in zip(f.gens, values):
                    f = f.eval(gen, value)

                return f
            else:
                j, a = 0, x
        else:
            j = f._gen_to_level(x)

        if not hasattr(f.rep, 'eval'):  # pragma: no cover
            raise OperationNotSupported(f, 'eval')

        try:
            result = f.rep.eval(a, j)
        except CoercionFailed:
            if not auto:
                raise DomainError("can't evaluate at %s in %s" % (a, f.rep.dom))
            else:
                a_domain, [a] = construct_domain([a])
                new_domain = f.get_domain().unify_with_symbols(a_domain, f.gens)

                f = f.set_domain(new_domain)
                a = new_domain.convert(a, a_domain)

                result = f.rep.eval(a, j)

        return f.per(result, remove=j)

    def __call__(f, *values):
        """
        Evaluate ``f`` at the give values.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y, z

        >>> f = Poly(2*x*y + 3*x + y + 2*z, x, y, z)

        >>> f(2)
        Poly(5*y + 2*z + 6, y, z, domain='ZZ')
        >>> f(2, 5)
        Poly(2*z + 31, z, domain='ZZ')
        >>> f(2, 5, 7)
        45

        """
        return f.eval(values)

    def half_gcdex(f, g, auto=True):
        """
        Half extended Euclidean algorithm of ``f`` and ``g``.

        Returns ``(s, h)`` such that ``h = gcd(f, g)`` and ``s*f = h (mod g)``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> f = x**4 - 2*x**3 - 6*x**2 + 12*x + 15
        >>> g = x**3 + x**2 - 4*x - 4

        >>> Poly(f).half_gcdex(Poly(g))
        (Poly(-1/5*x + 3/5, x, domain='QQ'), Poly(x + 1, x, domain='QQ'))

        """
        dom, per, F, G = f._unify(g)

        if auto and dom.is_Ring:
            F, G = F.to_field(), G.to_field()

        if hasattr(f.rep, 'half_gcdex'):
            s, h = F.half_gcdex(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'half_gcdex')

        return per(s), per(h)

    def gcdex(f, g, auto=True):
        """
        Extended Euclidean algorithm of ``f`` and ``g``.

        Returns ``(s, t, h)`` such that ``h = gcd(f, g)`` and ``s*f + t*g = h``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> f = x**4 - 2*x**3 - 6*x**2 + 12*x + 15
        >>> g = x**3 + x**2 - 4*x - 4

        >>> Poly(f).gcdex(Poly(g))
        (Poly(-1/5*x + 3/5, x, domain='QQ'),
         Poly(1/5*x**2 - 6/5*x + 2, x, domain='QQ'),
         Poly(x + 1, x, domain='QQ'))

        """
        dom, per, F, G = f._unify(g)

        if auto and dom.is_Ring:
            F, G = F.to_field(), G.to_field()

        if hasattr(f.rep, 'gcdex'):
            s, t, h = F.gcdex(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'gcdex')

        return per(s), per(t), per(h)

    def invert(f, g, auto=True):
        """
        Invert ``f`` modulo ``g`` when possible.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 - 1, x).invert(Poly(2*x - 1, x))
        Poly(-4/3, x, domain='QQ')

        >>> Poly(x**2 - 1, x).invert(Poly(x - 1, x))
        Traceback (most recent call last):
        ...
        NotInvertible: zero divisor

        """
        dom, per, F, G = f._unify(g)

        if auto and dom.is_Ring:
            F, G = F.to_field(), G.to_field()

        if hasattr(f.rep, 'invert'):
            result = F.invert(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'invert')

        return per(result)

    def revert(f, n):
        """
        Compute ``f**(-1)`` mod ``x**n``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(1, x).revert(2)
        Poly(1, x, domain='ZZ')

        >>> Poly(1 + x, x).revert(1)
        Poly(1, x, domain='ZZ')

        >>> Poly(x**2 - 1, x).revert(1)
        Traceback (most recent call last):
        ...
        NotReversible: only unity is reversible in a ring

        >>> Poly(1/x, x).revert(1)
        Traceback (most recent call last):
        ...
        PolynomialError: 1/x contains an element of the generators set

        """
        if hasattr(f.rep, 'revert'):
            result = f.rep.revert(int(n))
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'revert')

        return f.per(result)

    def subresultants(f, g):
        """
        Computes the subresultant PRS of ``f`` and ``g``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 + 1, x).subresultants(Poly(x**2 - 1, x))
        [Poly(x**2 + 1, x, domain='ZZ'),
         Poly(x**2 - 1, x, domain='ZZ'),
         Poly(-2, x, domain='ZZ')]

        """
        _, per, F, G = f._unify(g)

        if hasattr(f.rep, 'subresultants'):
            result = F.subresultants(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'subresultants')

        return list(map(per, result))

    def resultant(f, g, includePRS=False):
        """
        Computes the resultant of ``f`` and ``g`` via PRS.

        If includePRS=True, it includes the subresultant PRS in the result.
        Because the PRS is used to calculate the resultant, this is more
        efficient than calling :func:`subresultants` separately.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> f = Poly(x**2 + 1, x)

        >>> f.resultant(Poly(x**2 - 1, x))
        4
        >>> f.resultant(Poly(x**2 - 1, x), includePRS=True)
        (4, [Poly(x**2 + 1, x, domain='ZZ'), Poly(x**2 - 1, x, domain='ZZ'),
             Poly(-2, x, domain='ZZ')])

        """
        _, per, F, G = f._unify(g)

        if hasattr(f.rep, 'resultant'):
            if includePRS:
                result, R = F.resultant(G, includePRS=includePRS)
            else:
                result = F.resultant(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'resultant')

        if includePRS:
            return (per(result, remove=0), list(map(per, R)))
        return per(result, remove=0)

    def discriminant(f):
        """
        Computes the discriminant of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 + 2*x + 3, x).discriminant()
        -8

        """
        if hasattr(f.rep, 'discriminant'):
            result = f.rep.discriminant()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'discriminant')

        return f.per(result, remove=0)

    def dispersionset(f, g=None):
        r"""Compute the *dispersion set* of two polynomials.

        For two polynomials `f(x)` and `g(x)` with `\deg f > 0`
        and `\deg g > 0` the dispersion set `\operatorname{J}(f, g)` is defined as:

        .. math::
            \operatorname{J}(f, g)
            & := \\{a \in \mathbb{N}_0 | \gcd(f(x), g(x+a)) \neq 1\\} \\\\
            &  = \\{a \in \mathbb{N}_0 | \deg \gcd(f(x), g(x+a)) \geq 1\\}

        For a single polynomial one defines `\operatorname{J}(f) := \operatorname{J}(f, f)`.

        Examples
        ========

        >>> from sympy import poly
        >>> from sympy.polys.dispersion import dispersion, dispersionset
        >>> from sympy.abc import x

        Dispersion set and dispersion of a simple polynomial:

        >>> fp = poly((x - 3)*(x + 3), x)
        >>> sorted(dispersionset(fp))
        [0, 6]
        >>> dispersion(fp)
        6

        Note that the definition of the dispersion is not symmetric:

        >>> fp = poly(x**4 - 3*x**2 + 1, x)
        >>> gp = fp.shift(-3)
        >>> sorted(dispersionset(fp, gp))
        [2, 3, 4]
        >>> dispersion(fp, gp)
        4
        >>> sorted(dispersionset(gp, fp))
        []
        >>> dispersion(gp, fp)
        -oo

        Computing the dispersion also works over field extensions:

        >>> from sympy import sqrt
        >>> fp = poly(x**2 + sqrt(5)*x - 1, x, domain='QQ<sqrt(5)>')
        >>> gp = poly(x**2 + (2 + sqrt(5))*x + sqrt(5), x, domain='QQ<sqrt(5)>')
        >>> sorted(dispersionset(fp, gp))
        [2]
        >>> sorted(dispersionset(gp, fp))
        [1, 4]

        We can even perform the computations for polynomials
        having symbolic coefficients:

        >>> from sympy.abc import a
        >>> fp = poly(4*x**4 + (4*a + 8)*x**3 + (a**2 + 6*a + 4)*x**2 + (a**2 + 2*a)*x, x)
        >>> sorted(dispersionset(fp))
        [0, 1]

        See Also
        ========

        dispersion

        References
        ==========

        1. [ManWright94]_
        2. [Koepf98]_
        3. [Abramov71]_
        4. [Man93]_
        """
        from sympy.polys.dispersion import dispersionset
        return dispersionset(f, g)

    def dispersion(f, g=None):
        r"""Compute the *dispersion* of polynomials.

        For two polynomials `f(x)` and `g(x)` with `\deg f > 0`
        and `\deg g > 0` the dispersion `\operatorname{dis}(f, g)` is defined as:

        .. math::
            \operatorname{dis}(f, g)
            & := \max\\{ J(f,g) \cup \\{0\\} \\} \\\\
            &  = \max\\{ \\{a \in \mathbb{N} | \gcd(f(x), g(x+a)) \neq 1\\} \cup \\{0\\} \\}

        and for a single polynomial `\operatorname{dis}(f) := \operatorname{dis}(f, f)`.

        Examples
        ========

        >>> from sympy import poly
        >>> from sympy.polys.dispersion import dispersion, dispersionset
        >>> from sympy.abc import x

        Dispersion set and dispersion of a simple polynomial:

        >>> fp = poly((x - 3)*(x + 3), x)
        >>> sorted(dispersionset(fp))
        [0, 6]
        >>> dispersion(fp)
        6

        Note that the definition of the dispersion is not symmetric:

        >>> fp = poly(x**4 - 3*x**2 + 1, x)
        >>> gp = fp.shift(-3)
        >>> sorted(dispersionset(fp, gp))
        [2, 3, 4]
        >>> dispersion(fp, gp)
        4
        >>> sorted(dispersionset(gp, fp))
        []
        >>> dispersion(gp, fp)
        -oo

        Computing the dispersion also works over field extensions:

        >>> from sympy import sqrt
        >>> fp = poly(x**2 + sqrt(5)*x - 1, x, domain='QQ<sqrt(5)>')
        >>> gp = poly(x**2 + (2 + sqrt(5))*x + sqrt(5), x, domain='QQ<sqrt(5)>')
        >>> sorted(dispersionset(fp, gp))
        [2]
        >>> sorted(dispersionset(gp, fp))
        [1, 4]

        We can even perform the computations for polynomials
        having symbolic coefficients:

        >>> from sympy.abc import a
        >>> fp = poly(4*x**4 + (4*a + 8)*x**3 + (a**2 + 6*a + 4)*x**2 + (a**2 + 2*a)*x, x)
        >>> sorted(dispersionset(fp))
        [0, 1]

        See Also
        ========

        dispersionset

        References
        ==========

        1. [ManWright94]_
        2. [Koepf98]_
        3. [Abramov71]_
        4. [Man93]_
        """
        from sympy.polys.dispersion import dispersion
        return dispersion(f, g)

    def cofactors(f, g):
        """
        Returns the GCD of ``f`` and ``g`` and their cofactors.

        Returns polynomials ``(h, cff, cfg)`` such that ``h = gcd(f, g)``, and
        ``cff = quo(f, h)`` and ``cfg = quo(g, h)`` are, so called, cofactors
        of ``f`` and ``g``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 - 1, x).cofactors(Poly(x**2 - 3*x + 2, x))
        (Poly(x - 1, x, domain='ZZ'),
         Poly(x + 1, x, domain='ZZ'),
         Poly(x - 2, x, domain='ZZ'))

        """
        _, per, F, G = f._unify(g)

        if hasattr(f.rep, 'cofactors'):
            h, cff, cfg = F.cofactors(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'cofactors')

        return per(h), per(cff), per(cfg)

    def gcd(f, g):
        """
        Returns the polynomial GCD of ``f`` and ``g``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 - 1, x).gcd(Poly(x**2 - 3*x + 2, x))
        Poly(x - 1, x, domain='ZZ')

        """
        _, per, F, G = f._unify(g)

        if hasattr(f.rep, 'gcd'):
            result = F.gcd(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'gcd')

        return per(result)

    def lcm(f, g):
        """
        Returns polynomial LCM of ``f`` and ``g``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 - 1, x).lcm(Poly(x**2 - 3*x + 2, x))
        Poly(x**3 - 2*x**2 - x + 2, x, domain='ZZ')

        """
        _, per, F, G = f._unify(g)

        if hasattr(f.rep, 'lcm'):
            result = F.lcm(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'lcm')

        return per(result)

    def trunc(f, p):
        """
        Reduce ``f`` modulo a constant ``p``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(2*x**3 + 3*x**2 + 5*x + 7, x).trunc(3)
        Poly(-x**3 - x + 1, x, domain='ZZ')

        """
        p = f.rep.dom.convert(p)

        if hasattr(f.rep, 'trunc'):
            result = f.rep.trunc(p)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'trunc')

        return f.per(result)

    def monic(self, auto=True):
        """
        Divides all coefficients by ``LC(f)``.

        Examples
        ========

        >>> from sympy import Poly, ZZ
        >>> from sympy.abc import x

        >>> Poly(3*x**2 + 6*x + 9, x, domain=ZZ).monic()
        Poly(x**2 + 2*x + 3, x, domain='QQ')

        >>> Poly(3*x**2 + 4*x + 2, x, domain=ZZ).monic()
        Poly(x**2 + 4/3*x + 2/3, x, domain='QQ')

        """
        f = self

        if auto and f.rep.dom.is_Ring:
            f = f.to_field()

        if hasattr(f.rep, 'monic'):
            result = f.rep.monic()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'monic')

        return f.per(result)

    def content(f):
        """
        Returns the GCD of polynomial coefficients.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(6*x**2 + 8*x + 12, x).content()
        2

        """
        if hasattr(f.rep, 'content'):
            result = f.rep.content()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'content')

        return f.rep.dom.to_sympy(result)

    def primitive(f):
        """
        Returns the content and a primitive form of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(2*x**2 + 8*x + 12, x).primitive()
        (2, Poly(x**2 + 4*x + 6, x, domain='ZZ'))

        """
        if hasattr(f.rep, 'primitive'):
            cont, result = f.rep.primitive()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'primitive')

        return f.rep.dom.to_sympy(cont), f.per(result)

    def compose(f, g):
        """
        Computes the functional composition of ``f`` and ``g``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 + x, x).compose(Poly(x - 1, x))
        Poly(x**2 - x, x, domain='ZZ')

        """
        _, per, F, G = f._unify(g)

        if hasattr(f.rep, 'compose'):
            result = F.compose(G)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'compose')

        return per(result)

    def decompose(f):
        """
        Computes a functional decomposition of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**4 + 2*x**3 - x - 1, x, domain='ZZ').decompose()
        [Poly(x**2 - x - 1, x, domain='ZZ'), Poly(x**2 + x, x, domain='ZZ')]

        """
        if hasattr(f.rep, 'decompose'):
            result = f.rep.decompose()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'decompose')

        return list(map(f.per, result))

    def shift(f, a):
        """
        Efficiently compute Taylor shift ``f(x + a)``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 - 2*x + 1, x).shift(2)
        Poly(x**2 + 2*x + 1, x, domain='ZZ')

        """
        if hasattr(f.rep, 'shift'):
            result = f.rep.shift(a)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'shift')

        return f.per(result)

    def transform(f, p, q):
        """
        Efficiently evaluate the functional transformation ``q**n * f(p/q)``.


        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 - 2*x + 1, x).transform(Poly(x + 1, x), Poly(x - 1, x))
        Poly(4, x, domain='ZZ')

        """
        P, Q = p.unify(q)
        F, P = f.unify(P)
        F, Q = F.unify(Q)

        if hasattr(F.rep, 'transform'):
            result = F.rep.transform(P.rep, Q.rep)
        else:  # pragma: no cover
            raise OperationNotSupported(F, 'transform')

        return F.per(result)

    def sturm(self, auto=True):
        """
        Computes the Sturm sequence of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**3 - 2*x**2 + x - 3, x).sturm()
        [Poly(x**3 - 2*x**2 + x - 3, x, domain='QQ'),
         Poly(3*x**2 - 4*x + 1, x, domain='QQ'),
         Poly(2/9*x + 25/9, x, domain='QQ'),
         Poly(-2079/4, x, domain='QQ')]

        """
        f = self

        if auto and f.rep.dom.is_Ring:
            f = f.to_field()

        if hasattr(f.rep, 'sturm'):
            result = f.rep.sturm()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'sturm')

        return list(map(f.per, result))

    def gff_list(f):
        """
        Computes greatest factorial factorization of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> f = x**5 + 2*x**4 - x**3 - 2*x**2

        >>> Poly(f).gff_list()
        [(Poly(x, x, domain='ZZ'), 1), (Poly(x + 2, x, domain='ZZ'), 4)]

        """
        if hasattr(f.rep, 'gff_list'):
            result = f.rep.gff_list()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'gff_list')

        return [(f.per(g), k) for g, k in result]

    def norm(f):
        """
        Computes the product, ``Norm(f)``, of the conjugates of
        a polynomial ``f`` defined over a number field ``K``.

        Examples
        ========

        >>> from sympy import Poly, sqrt
        >>> from sympy.abc import x

        >>> a, b = sqrt(2), sqrt(3)

        A polynomial over a quadratic extension.
        Two conjugates x - a and x + a.

        >>> f = Poly(x - a, x, extension=a)
        >>> f.norm()
        Poly(x**2 - 2, x, domain='QQ')

        A polynomial over a quartic extension.
        Four conjugates x - a, x - a, x + a and x + a.

        >>> f = Poly(x - a, x, extension=(a, b))
        >>> f.norm()
        Poly(x**4 - 4*x**2 + 4, x, domain='QQ')

        """
        if hasattr(f.rep, 'norm'):
            r = f.rep.norm()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'norm')

        return f.per(r)

    def sqf_norm(f):
        """
        Computes square-free norm of ``f``.

        Returns ``s``, ``f``, ``r``, such that ``g(x) = f(x-sa)`` and
        ``r(x) = Norm(g(x))`` is a square-free polynomial over ``K``,
        where ``a`` is the algebraic extension of the ground domain.

        Examples
        ========

        >>> from sympy import Poly, sqrt
        >>> from sympy.abc import x

        >>> s, f, r = Poly(x**2 + 1, x, extension=[sqrt(3)]).sqf_norm()

        >>> s
        1
        >>> f
        Poly(x**2 - 2*sqrt(3)*x + 4, x, domain='QQ<sqrt(3)>')
        >>> r
        Poly(x**4 - 4*x**2 + 16, x, domain='QQ')

        """
        if hasattr(f.rep, 'sqf_norm'):
            s, g, r = f.rep.sqf_norm()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'sqf_norm')

        return s, f.per(g), f.per(r)

    def sqf_part(f):
        """
        Computes square-free part of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**3 - 3*x - 2, x).sqf_part()
        Poly(x**2 - x - 2, x, domain='ZZ')

        """
        if hasattr(f.rep, 'sqf_part'):
            result = f.rep.sqf_part()
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'sqf_part')

        return f.per(result)

    def sqf_list(f, all=False):
        """
        Returns a list of square-free factors of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> f = 2*x**5 + 16*x**4 + 50*x**3 + 76*x**2 + 56*x + 16

        >>> Poly(f).sqf_list()
        (2, [(Poly(x + 1, x, domain='ZZ'), 2),
             (Poly(x + 2, x, domain='ZZ'), 3)])

        >>> Poly(f).sqf_list(all=True)
        (2, [(Poly(1, x, domain='ZZ'), 1),
             (Poly(x + 1, x, domain='ZZ'), 2),
             (Poly(x + 2, x, domain='ZZ'), 3)])

        """
        if hasattr(f.rep, 'sqf_list'):
            coeff, factors = f.rep.sqf_list(all)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'sqf_list')

        return f.rep.dom.to_sympy(coeff), [(f.per(g), k) for g, k in factors]

    def sqf_list_include(f, all=False):
        """
        Returns a list of square-free factors of ``f``.

        Examples
        ========

        >>> from sympy import Poly, expand
        >>> from sympy.abc import x

        >>> f = expand(2*(x + 1)**3*x**4)
        >>> f
        2*x**7 + 6*x**6 + 6*x**5 + 2*x**4

        >>> Poly(f).sqf_list_include()
        [(Poly(2, x, domain='ZZ'), 1),
         (Poly(x + 1, x, domain='ZZ'), 3),
         (Poly(x, x, domain='ZZ'), 4)]

        >>> Poly(f).sqf_list_include(all=True)
        [(Poly(2, x, domain='ZZ'), 1),
         (Poly(1, x, domain='ZZ'), 2),
         (Poly(x + 1, x, domain='ZZ'), 3),
         (Poly(x, x, domain='ZZ'), 4)]

        """
        if hasattr(f.rep, 'sqf_list_include'):
            factors = f.rep.sqf_list_include(all)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'sqf_list_include')

        return [(f.per(g), k) for g, k in factors]

    def factor_list(f):
        """
        Returns a list of irreducible factors of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> f = 2*x**5 + 2*x**4*y + 4*x**3 + 4*x**2*y + 2*x + 2*y

        >>> Poly(f).factor_list()
        (2, [(Poly(x + y, x, y, domain='ZZ'), 1),
             (Poly(x**2 + 1, x, y, domain='ZZ'), 2)])

        """
        if hasattr(f.rep, 'factor_list'):
            try:
                coeff, factors = f.rep.factor_list()
            except DomainError:
                return S.One, [(f, 1)]
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'factor_list')

        return f.rep.dom.to_sympy(coeff), [(f.per(g), k) for g, k in factors]

    def factor_list_include(f):
        """
        Returns a list of irreducible factors of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> f = 2*x**5 + 2*x**4*y + 4*x**3 + 4*x**2*y + 2*x + 2*y

        >>> Poly(f).factor_list_include()
        [(Poly(2*x + 2*y, x, y, domain='ZZ'), 1),
         (Poly(x**2 + 1, x, y, domain='ZZ'), 2)]

        """
        if hasattr(f.rep, 'factor_list_include'):
            try:
                factors = f.rep.factor_list_include()
            except DomainError:
                return [(f, 1)]
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'factor_list_include')

        return [(f.per(g), k) for g, k in factors]

    def intervals(f, all=False, eps=None, inf=None, sup=None, fast=False, sqf=False):
        """
        Compute isolating intervals for roots of ``f``.

        For real roots the Vincent-Akritas-Strzebonski (VAS) continued fractions method is used.

        References
        ==========
        .. [#] Alkiviadis G. Akritas and Adam W. Strzebonski: A Comparative Study of Two Real Root
            Isolation Methods . Nonlinear Analysis: Modelling and Control, Vol. 10, No. 4, 297-304, 2005.
        .. [#] Alkiviadis G. Akritas, Adam W. Strzebonski and Panagiotis S. Vigklas: Improving the
            Performance of the Continued Fractions Method Using new Bounds of Positive Roots. Nonlinear
            Analysis: Modelling and Control, Vol. 13, No. 3, 265-279, 2008.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 - 3, x).intervals()
        [((-2, -1), 1), ((1, 2), 1)]
        >>> Poly(x**2 - 3, x).intervals(eps=1e-2)
        [((-26/15, -19/11), 1), ((19/11, 26/15), 1)]

        """
        if eps is not None:
            eps = QQ.convert(eps)

            if eps <= 0:
                raise ValueError("'eps' must be a positive rational")

        if inf is not None:
            inf = QQ.convert(inf)
        if sup is not None:
            sup = QQ.convert(sup)

        if hasattr(f.rep, 'intervals'):
            result = f.rep.intervals(
                all=all, eps=eps, inf=inf, sup=sup, fast=fast, sqf=sqf)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'intervals')

        if sqf:
            def _real(interval):
                s, t = interval
                return (QQ.to_sympy(s), QQ.to_sympy(t))

            if not all:
                return list(map(_real, result))

            def _complex(rectangle):
                (u, v), (s, t) = rectangle
                return (QQ.to_sympy(u) + I*QQ.to_sympy(v),
                        QQ.to_sympy(s) + I*QQ.to_sympy(t))

            real_part, complex_part = result

            return list(map(_real, real_part)), list(map(_complex, complex_part))
        else:
            def _real(interval):
                (s, t), k = interval
                return ((QQ.to_sympy(s), QQ.to_sympy(t)), k)

            if not all:
                return list(map(_real, result))

            def _complex(rectangle):
                ((u, v), (s, t)), k = rectangle
                return ((QQ.to_sympy(u) + I*QQ.to_sympy(v),
                         QQ.to_sympy(s) + I*QQ.to_sympy(t)), k)

            real_part, complex_part = result

            return list(map(_real, real_part)), list(map(_complex, complex_part))

    def refine_root(f, s, t, eps=None, steps=None, fast=False, check_sqf=False):
        """
        Refine an isolating interval of a root to the given precision.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 - 3, x).refine_root(1, 2, eps=1e-2)
        (19/11, 26/15)

        """
        if check_sqf and not f.is_sqf:
            raise PolynomialError("only square-free polynomials supported")

        s, t = QQ.convert(s), QQ.convert(t)

        if eps is not None:
            eps = QQ.convert(eps)

            if eps <= 0:
                raise ValueError("'eps' must be a positive rational")

        if steps is not None:
            steps = int(steps)
        elif eps is None:
            steps = 1

        if hasattr(f.rep, 'refine_root'):
            S, T = f.rep.refine_root(s, t, eps=eps, steps=steps, fast=fast)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'refine_root')

        return QQ.to_sympy(S), QQ.to_sympy(T)

    def count_roots(f, inf=None, sup=None):
        """
        Return the number of roots of ``f`` in ``[inf, sup]`` interval.

        Examples
        ========

        >>> from sympy import Poly, I
        >>> from sympy.abc import x

        >>> Poly(x**4 - 4, x).count_roots(-3, 3)
        2
        >>> Poly(x**4 - 4, x).count_roots(0, 1 + 3*I)
        1

        """
        inf_real, sup_real = True, True

        if inf is not None:
            inf = sympify(inf)

            if inf is S.NegativeInfinity:
                inf = None
            else:
                re, im = inf.as_real_imag()

                if not im:
                    inf = QQ.convert(inf)
                else:
                    inf, inf_real = list(map(QQ.convert, (re, im))), False

        if sup is not None:
            sup = sympify(sup)

            if sup is S.Infinity:
                sup = None
            else:
                re, im = sup.as_real_imag()

                if not im:
                    sup = QQ.convert(sup)
                else:
                    sup, sup_real = list(map(QQ.convert, (re, im))), False

        if inf_real and sup_real:
            if hasattr(f.rep, 'count_real_roots'):
                count = f.rep.count_real_roots(inf=inf, sup=sup)
            else:  # pragma: no cover
                raise OperationNotSupported(f, 'count_real_roots')
        else:
            if inf_real and inf is not None:
                inf = (inf, QQ.zero)

            if sup_real and sup is not None:
                sup = (sup, QQ.zero)

            if hasattr(f.rep, 'count_complex_roots'):
                count = f.rep.count_complex_roots(inf=inf, sup=sup)
            else:  # pragma: no cover
                raise OperationNotSupported(f, 'count_complex_roots')

        return Integer(count)

    def root(f, index, radicals=True):
        """
        Get an indexed root of a polynomial.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> f = Poly(2*x**3 - 7*x**2 + 4*x + 4)

        >>> f.root(0)
        -1/2
        >>> f.root(1)
        2
        >>> f.root(2)
        2
        >>> f.root(3)
        Traceback (most recent call last):
        ...
        IndexError: root index out of [-3, 2] range, got 3

        >>> Poly(x**5 + x + 1).root(0)
        CRootOf(x**3 - x**2 + 1, 0)

        """
        return sympy.polys.rootoftools.rootof(f, index, radicals=radicals)

    def real_roots(f, multiple=True, radicals=True):
        """
        Return a list of real roots with multiplicities.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(2*x**3 - 7*x**2 + 4*x + 4).real_roots()
        [-1/2, 2, 2]
        >>> Poly(x**3 + x + 1).real_roots()
        [CRootOf(x**3 + x + 1, 0)]

        """
        reals = sympy.polys.rootoftools.CRootOf.real_roots(f, radicals=radicals)

        if multiple:
            return reals
        else:
            return group(reals, multiple=False)

    def all_roots(f, multiple=True, radicals=True):
        """
        Return a list of real and complex roots with multiplicities.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(2*x**3 - 7*x**2 + 4*x + 4).all_roots()
        [-1/2, 2, 2]
        >>> Poly(x**3 + x + 1).all_roots()
        [CRootOf(x**3 + x + 1, 0),
         CRootOf(x**3 + x + 1, 1),
         CRootOf(x**3 + x + 1, 2)]

        """
        roots = sympy.polys.rootoftools.CRootOf.all_roots(f, radicals=radicals)

        if multiple:
            return roots
        else:
            return group(roots, multiple=False)

    def nroots(f, n=15, maxsteps=50, cleanup=True):
        """
        Compute numerical approximations of roots of ``f``.

        Parameters
        ==========

        n ... the number of digits to calculate
        maxsteps ... the maximum number of iterations to do

        If the accuracy `n` cannot be reached in `maxsteps`, it will raise an
        exception. You need to rerun with higher maxsteps.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 - 3).nroots(n=15)
        [-1.73205080756888, 1.73205080756888]
        >>> Poly(x**2 - 3).nroots(n=30)
        [-1.73205080756887729352744634151, 1.73205080756887729352744634151]

        """
        from sympy.functions.elementary.complexes import sign
        if f.is_multivariate:
            raise MultivariatePolynomialError(
                "can't compute numerical roots of %s" % f)

        if f.degree() <= 0:
            return []

        # For integer and rational coefficients, convert them to integers only
        # (for accuracy). Otherwise just try to convert the coefficients to
        # mpmath.mpc and raise an exception if the conversion fails.
        if f.rep.dom is ZZ:
            coeffs = [int(coeff) for coeff in f.all_coeffs()]
        elif f.rep.dom is QQ:
            denoms = [coeff.q for coeff in f.all_coeffs()]
            from sympy.core.numbers import ilcm
            fac = ilcm(*denoms)
            coeffs = [int(coeff*fac) for coeff in f.all_coeffs()]
        else:
            coeffs = [coeff.evalf(n=n).as_real_imag()
                    for coeff in f.all_coeffs()]
            try:
                coeffs = [mpmath.mpc(*coeff) for coeff in coeffs]
            except TypeError:
                raise DomainError("Numerical domain expected, got %s" % \
                        f.rep.dom)

        dps = mpmath.mp.dps
        mpmath.mp.dps = n

        try:
            # We need to add extra precision to guard against losing accuracy.
            # 10 times the degree of the polynomial seems to work well.
            roots = mpmath.polyroots(coeffs, maxsteps=maxsteps,
                    cleanup=cleanup, error=False, extraprec=f.degree()*10)

            # Mpmath puts real roots first, then complex ones (as does all_roots)
            # so we make sure this convention holds here, too.
            roots = list(map(sympify,
                sorted(roots, key=lambda r: (1 if r.imag else 0, r.real, abs(r.imag), sign(r.imag)))))
        except NoConvergence:
            raise NoConvergence(
                'convergence to root failed; try n < %s or maxsteps > %s' % (
                n, maxsteps))
        finally:
            mpmath.mp.dps = dps

        return roots

    def ground_roots(f):
        """
        Compute roots of ``f`` by factorization in the ground domain.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**6 - 4*x**4 + 4*x**3 - x**2).ground_roots()
        {0: 2, 1: 2}

        """
        if f.is_multivariate:
            raise MultivariatePolynomialError(
                "can't compute ground roots of %s" % f)

        roots = {}

        for factor, k in f.factor_list()[1]:
            if factor.is_linear:
                a, b = factor.all_coeffs()
                roots[-b/a] = k

        return roots

    def nth_power_roots_poly(f, n):
        """
        Construct a polynomial with n-th powers of roots of ``f``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> f = Poly(x**4 - x**2 + 1)

        >>> f.nth_power_roots_poly(2)
        Poly(x**4 - 2*x**3 + 3*x**2 - 2*x + 1, x, domain='ZZ')
        >>> f.nth_power_roots_poly(3)
        Poly(x**4 + 2*x**2 + 1, x, domain='ZZ')
        >>> f.nth_power_roots_poly(4)
        Poly(x**4 + 2*x**3 + 3*x**2 + 2*x + 1, x, domain='ZZ')
        >>> f.nth_power_roots_poly(12)
        Poly(x**4 - 4*x**3 + 6*x**2 - 4*x + 1, x, domain='ZZ')

        """
        if f.is_multivariate:
            raise MultivariatePolynomialError(
                "must be a univariate polynomial")

        N = sympify(n)

        if N.is_Integer and N >= 1:
            n = int(N)
        else:
            raise ValueError("'n' must an integer and n >= 1, got %s" % n)

        x = f.gen
        t = Dummy('t')

        r = f.resultant(f.__class__.from_expr(x**n - t, x, t))

        return r.replace(t, x)

    def cancel(f, g, include=False):
        """
        Cancel common factors in a rational function ``f/g``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(2*x**2 - 2, x).cancel(Poly(x**2 - 2*x + 1, x))
        (1, Poly(2*x + 2, x, domain='ZZ'), Poly(x - 1, x, domain='ZZ'))

        >>> Poly(2*x**2 - 2, x).cancel(Poly(x**2 - 2*x + 1, x), include=True)
        (Poly(2*x + 2, x, domain='ZZ'), Poly(x - 1, x, domain='ZZ'))

        """
        dom, per, F, G = f._unify(g)

        if hasattr(F, 'cancel'):
            result = F.cancel(G, include=include)
        else:  # pragma: no cover
            raise OperationNotSupported(f, 'cancel')

        if not include:
            if dom.has_assoc_Ring:
                dom = dom.get_ring()

            cp, cq, p, q = result

            cp = dom.to_sympy(cp)
            cq = dom.to_sympy(cq)

            return cp/cq, per(p), per(q)
        else:
            return tuple(map(per, result))

    @property
    def is_zero(f):
        """
        Returns ``True`` if ``f`` is a zero polynomial.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(0, x).is_zero
        True
        >>> Poly(1, x).is_zero
        False

        """
        return f.rep.is_zero

    @property
    def is_one(f):
        """
        Returns ``True`` if ``f`` is a unit polynomial.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(0, x).is_one
        False
        >>> Poly(1, x).is_one
        True

        """
        return f.rep.is_one

    @property
    def is_sqf(f):
        """
        Returns ``True`` if ``f`` is a square-free polynomial.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 - 2*x + 1, x).is_sqf
        False
        >>> Poly(x**2 - 1, x).is_sqf
        True

        """
        return f.rep.is_sqf

    @property
    def is_monic(f):
        """
        Returns ``True`` if the leading coefficient of ``f`` is one.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x + 2, x).is_monic
        True
        >>> Poly(2*x + 2, x).is_monic
        False

        """
        return f.rep.is_monic

    @property
    def is_primitive(f):
        """
        Returns ``True`` if GCD of the coefficients of ``f`` is one.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(2*x**2 + 6*x + 12, x).is_primitive
        False
        >>> Poly(x**2 + 3*x + 6, x).is_primitive
        True

        """
        return f.rep.is_primitive

    @property
    def is_ground(f):
        """
        Returns ``True`` if ``f`` is an element of the ground domain.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x, x).is_ground
        False
        >>> Poly(2, x).is_ground
        True
        >>> Poly(y, x).is_ground
        True

        """
        return f.rep.is_ground

    @property
    def is_linear(f):
        """
        Returns ``True`` if ``f`` is linear in all its variables.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x + y + 2, x, y).is_linear
        True
        >>> Poly(x*y + 2, x, y).is_linear
        False

        """
        return f.rep.is_linear

    @property
    def is_quadratic(f):
        """
        Returns ``True`` if ``f`` is quadratic in all its variables.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x*y + 2, x, y).is_quadratic
        True
        >>> Poly(x*y**2 + 2, x, y).is_quadratic
        False

        """
        return f.rep.is_quadratic

    @property
    def is_monomial(f):
        """
        Returns ``True`` if ``f`` is zero or has only one term.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(3*x**2, x).is_monomial
        True
        >>> Poly(3*x**2 + 1, x).is_monomial
        False

        """
        return f.rep.is_monomial

    @property
    def is_homogeneous(f):
        """
        Returns ``True`` if ``f`` is a homogeneous polynomial.

        A homogeneous polynomial is a polynomial whose all monomials with
        non-zero coefficients have the same total degree. If you want not
        only to check if a polynomial is homogeneous but also compute its
        homogeneous order, then use :func:`Poly.homogeneous_order`.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x**2 + x*y, x, y).is_homogeneous
        True
        >>> Poly(x**3 + x*y, x, y).is_homogeneous
        False

        """
        return f.rep.is_homogeneous

    @property
    def is_irreducible(f):
        """
        Returns ``True`` if ``f`` has no factors over its domain.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> Poly(x**2 + x + 1, x, modulus=2).is_irreducible
        True
        >>> Poly(x**2 + 1, x, modulus=2).is_irreducible
        False

        """
        return f.rep.is_irreducible

    @property
    def is_univariate(f):
        """
        Returns ``True`` if ``f`` is a univariate polynomial.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x**2 + x + 1, x).is_univariate
        True
        >>> Poly(x*y**2 + x*y + 1, x, y).is_univariate
        False
        >>> Poly(x*y**2 + x*y + 1, x).is_univariate
        True
        >>> Poly(x**2 + x + 1, x, y).is_univariate
        False

        """
        return len(f.gens) == 1

    @property
    def is_multivariate(f):
        """
        Returns ``True`` if ``f`` is a multivariate polynomial.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x**2 + x + 1, x).is_multivariate
        False
        >>> Poly(x*y**2 + x*y + 1, x, y).is_multivariate
        True
        >>> Poly(x*y**2 + x*y + 1, x).is_multivariate
        False
        >>> Poly(x**2 + x + 1, x, y).is_multivariate
        True

        """
        return len(f.gens) != 1

    @property
    def is_cyclotomic(f):
        """
        Returns ``True`` if ``f`` is a cyclotomic polnomial.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x

        >>> f = x**16 + x**14 - x**10 + x**8 - x**6 + x**2 + 1

        >>> Poly(f).is_cyclotomic
        False

        >>> g = x**16 + x**14 - x**10 - x**8 - x**6 + x**2 + 1

        >>> Poly(g).is_cyclotomic
        True

        """
        return f.rep.is_cyclotomic

    def __abs__(f):
        return f.abs()

    def __neg__(f):
        return f.neg()

    @_sympifyit('g', NotImplemented)
    def __add__(f, g):
        if not g.is_Poly:
            try:
                g = f.__class__(g, *f.gens)
            except PolynomialError:
                return f.as_expr() + g

        return f.add(g)

    @_sympifyit('g', NotImplemented)
    def __radd__(f, g):
        if not g.is_Poly:
            try:
                g = f.__class__(g, *f.gens)
            except PolynomialError:
                return g + f.as_expr()

        return g.add(f)

    @_sympifyit('g', NotImplemented)
    def __sub__(f, g):
        if not g.is_Poly:
            try:
                g = f.__class__(g, *f.gens)
            except PolynomialError:
                return f.as_expr() - g

        return f.sub(g)

    @_sympifyit('g', NotImplemented)
    def __rsub__(f, g):
        if not g.is_Poly:
            try:
                g = f.__class__(g, *f.gens)
            except PolynomialError:
                return g - f.as_expr()

        return g.sub(f)

    @_sympifyit('g', NotImplemented)
    def __mul__(f, g):
        if not g.is_Poly:
            try:
                g = f.__class__(g, *f.gens)
            except PolynomialError:
                return f.as_expr()*g

        return f.mul(g)

    @_sympifyit('g', NotImplemented)
    def __rmul__(f, g):
        if not g.is_Poly:
            try:
                g = f.__class__(g, *f.gens)
            except PolynomialError:
                return g*f.as_expr()

        return g.mul(f)

    @_sympifyit('n', NotImplemented)
    def __pow__(f, n):
        if n.is_Integer and n >= 0:
            return f.pow(n)
        else:
            return f.as_expr()**n

    @_sympifyit('g', NotImplemented)
    def __divmod__(f, g):
        if not g.is_Poly:
            g = f.__class__(g, *f.gens)

        return f.div(g)

    @_sympifyit('g', NotImplemented)
    def __rdivmod__(f, g):
        if not g.is_Poly:
            g = f.__class__(g, *f.gens)

        return g.div(f)

    @_sympifyit('g', NotImplemented)
    def __mod__(f, g):
        if not g.is_Poly:
            g = f.__class__(g, *f.gens)

        return f.rem(g)

    @_sympifyit('g', NotImplemented)
    def __rmod__(f, g):
        if not g.is_Poly:
            g = f.__class__(g, *f.gens)

        return g.rem(f)

    @_sympifyit('g', NotImplemented)
    def __floordiv__(f, g):
        if not g.is_Poly:
            g = f.__class__(g, *f.gens)

        return f.quo(g)

    @_sympifyit('g', NotImplemented)
    def __rfloordiv__(f, g):
        if not g.is_Poly:
            g = f.__class__(g, *f.gens)

        return g.quo(f)

    @_sympifyit('g', NotImplemented)
    def __div__(f, g):
        return f.as_expr()/g.as_expr()

    @_sympifyit('g', NotImplemented)
    def __rdiv__(f, g):
        return g.as_expr()/f.as_expr()

    __truediv__ = __div__
    __rtruediv__ = __rdiv__

    @_sympifyit('other', NotImplemented)
    def __eq__(self, other):
        f, g = self, other

        if not g.is_Poly:
            try:
                g = f.__class__(g, f.gens, domain=f.get_domain())
            except (PolynomialError, DomainError, CoercionFailed):
                return False

        if f.gens != g.gens:
            return False

        if f.rep.dom != g.rep.dom:
            try:
                dom = f.rep.dom.unify(g.rep.dom, f.gens)
            except UnificationFailed:
                return False

            f = f.set_domain(dom)
            g = g.set_domain(dom)

        return f.rep == g.rep

    @_sympifyit('g', NotImplemented)
    def __ne__(f, g):
        return not f == g

    def __nonzero__(f):
        return not f.is_zero

    __bool__ = __nonzero__

    def eq(f, g, strict=False):
        if not strict:
            return f == g
        else:
            return f._strict_eq(sympify(g))

    def ne(f, g, strict=False):
        return not f.eq(g, strict=strict)

    def _strict_eq(f, g):
        return isinstance(g, f.__class__) and f.gens == g.gens and f.rep.eq(g.rep, strict=True)

Poly({1:1,2:1},gens=S('sqrt(x)'))
$\displaystyle \operatorname{Poly}{\left( x + \sqrt{x}, \sqrt{x}, domain=\mathbb{Z} \right)}$
x=symbols('x', positive=True)
S('sqrt(x^2)-x')
$\displaystyle - x + \sqrt{x^{2}}$
str(Poly({1:1,2:1,4:1},gens=S('sqrt(x)')))
"Poly((sqrt(x))**4 + (sqrt(x))**2 + (sqrt(x)), sqrt(x), domain='ZZ')"
str(Poly('x^2+x+sqrt(x)'))
"Poly(x**2 + x + (sqrt(x)), x, sqrt(x), domain='ZZ')"
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)
    tup=tuple(fmin(f2,values))
    return tuple([x*x for x in tup])
formula=cyclize(Sm('((a-b)/c)^2-8**(1/2)*(a-b)/c'))
formula=(makesubs(formula,'[b,oo],[c,oo]'))
values=findvalues(formula)
display(values)
nsimplify(values,tolerance=0.1,rational=True)
values[0]/values[2]
Substitute $a\to b+y_{1}$
Substitute $b\to c+z_{1}$
Optimization terminated successfully.
         Current function value: 1.000000
         Iterations: 137
         Function evaluations: 249
(1.7908873553542452e-10, 2.5326984818340415e-10, 7.129450063690368)
2.5119572187973787e-11
from random import random
formula=cyclize(Sm('((a-b)/c)^2-8**(1/2)*(a-b)/c'))
formula=(makesubs(formula,'[b,oo],[c,oo]'))
num,den=fractioncancel(formula)
fs=sorted(num.free_symbols,key=str)
num=num.subs(zip(fs,list(map(lambda x:x**2,fs))))
num=Poly(num,domain='RR')
num1=(num.abs()+num)
num2=(num.abs()-num)
newformula=S(num1/num2)
def evaluate(x):
    return newformula.evalf(subs=dict(zip(fs,(1,1,1))))
def evaluate2(x):
    return num1.eval(dict(zip(fs,x)))/num2.eval(dict(zip(fs,x)))
#display(num2.eval(dict(zip(fs,(1,1,1)))))
print('ok')
%timeit evaluate((random(),random(),random()))
%timeit evaluate2((random(),random(),random()))
Poly('x^2+sqrt(x)').eval
Substitute $a\to b+e_{2}$
Substitute $b\to c+f_{2}$
ok
100 loops, best of 5: 8.66 ms per loop
1000 loops, best of 5: 1.31 ms per loop
<bound method Poly.eval of Poly(x**2 + (sqrt(x)), x, sqrt(x), domain='ZZ')>
evaluate4=lambdify(fs,newformula)
evaluate5=lambda x:evaluate4(*x)
%timeit evaluate4(*(random(),random(),random()))
%timeit evaluate5((random(),random(),random()))
from sympy import *
print('first attempt:',Poly('x^2+x+sqrt(x)'))
x=Symbol('x', positive=True)
print('second attempt:',Poly(x**2+x+sqrt(x)))
print('third attempt:',str(Poly({1:1,2:1,4:1},gens=S('sqrt(x)'))))
(Poly({1:1,2:1,4:1},gens=S('sqrt(x)'))-Poly('x^2')).as_expr()
$\displaystyle \sqrt{x} + x$
from sympy.solvers.solvers import unrad
unrad('sqrt(x)+x')
(x*(1 - x), [])
formula=('x^(2/3)+x^(3/4)+sqrt(x+y)')
type((Poly(formula).gens[1]).args[0])==Symbol
True
S('sqrt(x)').args
(x, 1/2)
def _powr(formula):
	if formula.func==Pow:
		return formula.args
	else:
		return [formula,S('1')]
pol=Poly('x**2+y**2+sqrt(x+y)+x**(1/2)+x**(1/3)')
newgens={}
for gen in pol.gens:
    base,pw=_powr(gen)
    num,den=fraction(pw)
    if (gen[0],num) in newgens:
        newgens[gen[0]]=
    else:
        newgens[gen[0]]=gen[1]
  File "<ipython-input-58-578d4067f943>", line 12
    newgens[gen[0]]=
                    ^
SyntaxError: invalid syntax
isinstance(S('1'),Rational)
True
6 in {3:4,5:6}
False
from math import *
from sympy import *
gcd(S('a/b'),S('c/d'))
$\displaystyle 1$
Poly('x^(2*pi)')
$\displaystyle \operatorname{Poly}{\left( x^{2 \pi}, x^{\pi}, domain=\mathbb{Z} \right)}$
gcd(S('1/8'),S('1/4'))
$\displaystyle \frac{1}{8}$
S('sqrt(2)*3*x').as_coeff_mul()
(3, (x, sqrt(2)))
def _powr(formula):
	if formula.func==Pow:
		return formula.args
	else:
		return [formula,S('1')]
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
print(reducegens(S('x^2+x+sqrt(x)')))
print(reducegens(S('x**2+y**2+x**(1/2)+x**(1/3)')))
print(reducegens(S('sqrt(x)+sqrt(z)+sqrt(x*z)')))
print(reducegens(S('sqrt(x)+sqrt(z)+sqrt(x)*sqrt(z)')))
print(reducegens(S('sqrt(x+y)+sqrt(z+t)+sqrt((x+y)*(z+t))')))
print(reducegens(S('sqrt(2)+sqrt(3)+sqrt(6)')))
Poly((sqrt(x))**4 + (sqrt(x))**2 + (sqrt(x)), sqrt(x), domain='ZZ')
Poly((x**(1/6))**12 + (x**(1/6))**3 + (x**(1/6))**2 + y**2, x**(1/6), y, domain='ZZ')
Poly((sqrt(x)) + (sqrt(x*z)) + (sqrt(z)), sqrt(x), sqrt(x*z), sqrt(z), domain='ZZ')
Poly((sqrt(x))*(sqrt(z)) + (sqrt(x)) + (sqrt(z)), sqrt(x), sqrt(z), domain='ZZ')
Poly((sqrt(t + z)) + (sqrt(t*x + t*y + x*z + y*z)) + (sqrt(x + y)), sqrt(t + z), sqrt(t*x + t*y + x*z + y*z), sqrt(x + y), domain='ZZ')
Poly((sqrt(2)) + (sqrt(3)) + (sqrt(6)), sqrt(2), sqrt(3), sqrt(6), domain='ZZ')
S('-2').as_coeff_mul()
x=Symbol('x', positive=True)
y=Symbol('y', positive=True)
Poly(sqrt(x)+sqrt(y)+sqrt(x*y))
sqrt(x*y)
simplify(S('sqrt(x*y)-sqrt(x)*sqrt(y)'))
x=symbols('x',positive=True)
ask(Q.positive(x))
Poly('x^2+x')+Poly('-x+7')
print(reducegens(S('sqrt(8)')))
sqrt(-1*-1)
Poly('sqrt(2*x)+sqrt(3*x)+sqrt(6*x)')
Poly('sqrt(2*pi)+sqrt(3*pi)')
solve(S('sqrt(6*x)')-S('sqrt(2*y)'))
Poly('sqrt(6)+sqrt(10)+sqrt(14)+sqrt(15)+sqrt(21)+sqrt(35)+sqrt()')
from shiroindev import *
formula=S('(a+b-c)^2')
prove(makesubs(formula,'[c,oo]'))
prove(makesubs(formula,'[a,oo]',variables='c'))
formula=S('(a+b-d)^2+(a+b-c)^2')
prove(makesubs(formula,'[a+b,oo],[a+b,oo]',variables='c,d'))
prove(makesubs(formula,'[0,a+b],[a+b,oo]',variables='c,d'))
prove(makesubs(formula,'[0,a+b],[0,a+b]',variables='c,d'))
findvalues(formula)
display(cyclize('((a-b)/c)^2'))
display(S('sqrt(8)')*cyclize('(a-b)/c'))
formula=cyclize('((a-b)/c)^2-8**(1/2)*(a-b)/c')
display(Latex('Case $a\ge c\ge b$'))
formula1=makesubs(formula,'[c,oo],[b,oo]',variables='a,c,b')
prove(formula1)
display(Latex('Case $a\ge b\ge c$'))
formula2=makesubs(formula,'[b,oo],[c,oo]')
prove(formula2*4,values='2**(1/2),1,1')
formula=cyclize('((a-b)/c)^2-8**(1/2)*(a-b)/c')
formula2=makesubs(formula,'[c,oo],[c,oo]')
findvalues(formula2,values=[5,2,7])
x,y=symbols('x,y')
ask(Q.positive(x**2+1),Q.real(x))
formula=Sm('(a^2+b^2+c^2)^2- 3(a^3b+b^3c+c^3a)')
formula1=makesubs(formula,'[b,oo],[c,oo]',variables='a,b')
formula2=makesubs(formula,'[c,oo],[b,oo]',variables='a,c')
prove(formula2)
print(findvalues(formula1))
#prove(formula1,values='4*b*(1-b)-1,4*c*(1-c)-1,1')
makesubs(Sm('a^2-2ab+bc-c^2+ca'),'[c,oo],[a,oo]',variables='b,c')

formula=S('sqrt(x*y)-sqrt(x)*sqrt(y)')
display(formula)
def assumeall(formula,**kwargs):
    formula=S(formula)
    fs=formula.free_symbols
    for x in fs:
        y=Symbol(str(x),**kwargs)
        formula=formula.subs(x,y)
    return formula
display(assumeall(formula,positive=True))
$\displaystyle - \sqrt{x} \sqrt{y} + \sqrt{x y}$
$\displaystyle 0$
Poly(assumeall('x+sqrt(x)',positive=True))
str(type(S('f1(x,y,z)')))=='f1'
ReprPrinter?
srepr(S('x^2+f(x)+f(x,y)+f()+sqrt(y)'))
from sympy import Function,srepr,S
import re
f=Function('f')
def allsymbols(formula):
    formula=S(formula)
    funcsymbols=[x[10:-2] for x in re.findall(r"Function\(\'.*?\'\)",srepr(formula))]
    return set(funcsymbols)|set(map(str,formula.free_symbols))
def vargen(n):
    x='abcdefghijklmnopqrstuvwxyz'[n%26]
    if n>=26:
        x+=str(n//26)
    return x
allsymbols(S('f(5,7)+f(5)+g(8)+sqrt(x*y)'))
vargen(36)
from sympy import *
from sympy.printing.repr import *
ReprPrinter??
{5,7} | {8,7}
[x[10:-2] for x in re.findall(r"Function\(\'.*?\'\)",srepr(S('x^2')))]
S('f(5,7)+f(5)+g(8)+sqrt(x*y)')
srepr(S('f(5,7)+f(5)+g(8)+sqrt(x*y)'))
display(Latex('$'+str(S('z12345'))+'$'))
'xyz'[-2::]
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-68-247357fa9171> in <module>()
      1 'xyz'[-2::]
----> 2 int('')

ValueError: invalid literal for int() with base 10: ''
def sortkey(x):
    x=str(x)
    i=len(x)-1
    while i>=0:
        if x[i] not in '0123456789':
            break
        i-=1
    if i==len(x)-1:
        return (-1,x)
    return (int(x[i+1:]),x[:i+1])
sorted(S('x,y,z1,t2,qq1,x0'),key=sortkey)
[x, y, x0, qq1, z1, t2]
S('x_1^2+x_2^2+x_3^2+x_4^2+x_5^2')
formula=S('x_1^2+x_2^2+x_3^2+x_4^2+x_5^2- (x_1*x_2+x_2*x_3+x_3*x_4+x_4*x_5)')
values=findvalues(formula)
values=[value/values[0] for value in values]
[value**2 for value in values]
Optimization terminated successfully.
         Current function value: 1.154701
         Iterations: 85
         Function evaluations: 152
[1.0,
 3.000209557356011,
 3.9998018341331227,
 2.9997236704279775,
 0.999903278712487]
formula1=S('(x_1^2+x_2^2+x_3^2+x_4^2+x_5^2)/(x_1*x_2+x_2*x_3+x_3*x_4+x_4*x_5)')
formula1.subs(S('[[x_1,1],[x_2,sqrt(3)],[x_3,2],[x_4,sqrt(3)],[x_5,1]]'))
$\displaystyle \frac{2 \sqrt{3}}{3}$
formula2=S('x_1^2+x_2^2+x_3^2+x_4^2+x_5^2- (2/sqrt(3))*(x_1*x_2+x_2*x_3+x_3*x_4+x_4*x_5)')
prove(formula2,values='1,sqrt(3),2,sqrt(3),1')
Substitute $x_{2}\to \sqrt{3} g$
Substitute $x_{3}\to 2 h$
Substitute $x_{4}\to \sqrt{3} i$
numerator: $3 g^{2} - 4 g h - 2 g x_{1} + 4 h^{2} - 4 h i + 3 i^{2} - 2 i x_{5} + x_{1}^{2} + x_{5}^{2}$
denominator: $1$
status: 0
From weighted AM-GM inequality:
$$4 g h \le 2 g^{2}+2 h^{2}$$
$$2 g x_{1} \le g^{2}+x_{1}^{2}$$
$$4 h i \le 2 h^{2}+2 i^{2}$$
$$2 i x_{5} \le i^{2}+x_{5}^{2}$$
$$ 0 \le 0 $$
The sum of all inequalities gives us a proof of the inequality.
0
Poly('a^(n+k)')
$\displaystyle \operatorname{Poly}{\left( a^{k}a^{n}, a^{k}, a^{n}, domain=\mathbb{Z} \right)}$
S('a^{n+k}/{b^n}+{b^(n+k)/(c^n)+c^(n+k)/{a^n}-( a^k+b^k+c^k)')
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-13-ca88f614d264> in <module>()
----> 1 S('a^{n+k}/{b^n}+{b^{n+k}}/{c^n}+{c^{n+k}}/{a^n}-( a^k+b^k+c^k)')

/home/grzegorz/Pobrane/SageMath/local/lib/python3.7/site-packages/sympy/core/sympify.py in sympify(a, locals, convert_xor, strict, rational, evaluate)
    382     try:
    383         a = a.replace('\n', '')
--> 384         expr = parse_expr(a, local_dict=locals, transformations=transformations, evaluate=evaluate)
    385     except (TokenError, SyntaxError) as exc:
    386         raise SympifyError('could not parse %r' % a, exc)

/home/grzegorz/Pobrane/SageMath/local/lib/python3.7/site-packages/sympy/parsing/sympy_parser.py in parse_expr(s, local_dict, transformations, global_dict, evaluate)
    964         code = compile(evaluateFalse(code), '<string>', 'eval')
    965 
--> 966     return eval_expr(code, local_dict, global_dict)
    967 
    968 

/home/grzegorz/Pobrane/SageMath/local/lib/python3.7/site-packages/sympy/parsing/sympy_parser.py in eval_expr(code, local_dict, global_dict)
    877     """
    878     expr = eval(
--> 879         code, global_dict, local_dict)  # take local objects in preference
    880 
    881     return expr

<string> in <module>()

/home/grzegorz/Pobrane/SageMath/local/lib/python3.7/site-packages/sympy/core/expr.py in __pow__(self, other, mod)
    163     def __pow__(self, other, mod=None):
    164         if mod is None:
--> 165             return self._pow(other)
    166         try:
    167             _self, other, mod = as_int(self), as_int(other), as_int(mod)

/home/grzegorz/Pobrane/SageMath/local/lib/python3.7/site-packages/sympy/core/decorators.py in __sympifyit_wrapper(a, b)
     89                 if not hasattr(b, '_op_priority'):
     90                     b = sympify(b, strict=True)
---> 91                 return func(a, b)
     92             except SympifyError:
     93                 return retval

/home/grzegorz/Pobrane/SageMath/local/lib/python3.7/site-packages/sympy/core/decorators.py in binary_op_wrapper(self, other)
    127                     if f is not None:
    128                         return f(self)
--> 129             return func(self, other)
    130         return binary_op_wrapper
    131     return priority_decorator

/home/grzegorz/Pobrane/SageMath/local/lib/python3.7/site-packages/sympy/core/expr.py in _pow(self, other)
    159     @call_highest_priority('__rpow__')
    160     def _pow(self, other):
--> 161         return Pow(self, other)
    162 
    163     def __pow__(self, other, mod=None):

/home/grzegorz/Pobrane/SageMath/local/lib/python3.7/site-packages/sympy/core/cache.py in wrapper(*args, **kwargs)
     92             def wrapper(*args, **kwargs):
     93                 try:
---> 94                     retval = cfunc(*args, **kwargs)
     95                 except TypeError:
     96                     retval = func(*args, **kwargs)

/home/grzegorz/Pobrane/SageMath/local/lib/python3.7/site-packages/sympy/core/power.py in __new__(cls, b, e, evaluate)
    274                 if not e.is_Atom and b is not S.Exp1 and not isinstance(b, exp_polar):
    275                     from sympy import numer, denom, log, sign, im, factor_terms
--> 276                     c, ex = factor_terms(e, sign=False).as_coeff_Mul()
    277                     den = denom(ex)
    278                     if isinstance(den, log) and den.args[0] == b:

AttributeError: 'FiniteSet' object has no attribute 'as_coeff_Mul'
import re
s=r"""99-my-name-is-John-Smith-6376827-%^-1-2-767980716"""
re.compile(r"^(.*?)-").search(str(s[::-1])).group(1)[::-1]
'767980716'
from sympy import *
from sympy.parsing.latex import parse_latex
from shiroindev import *
shiro.seed=1
from IPython.display import Latex
shiro.display=lambda x:display(Latex(x))
formula=expand(parse_latex(r'(\frac{a^{n+k}}{b^n}+\frac{b^{n+k}}{c^n}+\frac{c^{n+k}}{a^n}-( a^k+b^k+c^k))*a^nb^nc^n'))
newproof()
prove(formula)
prove(makesubs(formula,['k',oo],variables='n'))
ANTLR runtime and generated code versions disagree: 4.8!=4.7.1
ANTLR runtime and generated code versions disagree: 4.8!=4.7.1
numerator: $a^{k} a^{2 n} c^{n} - a^{k} a^{n} b^{n} c^{n} + a^{n} b^{k} b^{2 n} - a^{n} b^{k} b^{n} c^{n} - a^{n} b^{n} c^{k} c^{n} + b^{n} c^{k} c^{2 n}$
denominator: $1$
status: 2
Program couldn't find any proof.
$$ a^{k} a^{n} b^{n} c^{n}+a^{n} b^{k} b^{n} c^{n}+a^{n} b^{n} c^{k} c^{n} \le a^{k} a^{2 n} c^{n}+a^{n} b^{k} b^{2 n}+b^{n} c^{k} c^{2 n} $$
Substitute $n\to d + k$
numerator: $a^{2 d} a^{3 k} c^{d} c^{k} - a^{d} a^{2 k} b^{d} b^{k} c^{d} c^{k} + a^{d} a^{k} b^{2 d} b^{3 k} - a^{d} a^{k} b^{d} b^{2 k} c^{d} c^{k} - a^{d} a^{k} b^{d} b^{k} c^{d} c^{2 k} + b^{d} b^{k} c^{2 d} c^{3 k}$
denominator: $1$
status: 2
Program couldn't find any proof.
$$ a^{d} a^{2 k} b^{d} b^{k} c^{d} c^{k}+a^{d} a^{k} b^{d} b^{2 k} c^{d} c^{k}+a^{d} a^{k} b^{d} b^{k} c^{d} c^{2 k} \le a^{2 d} a^{3 k} c^{d} c^{k}+a^{d} a^{k} b^{2 d} b^{3 k}+b^{d} b^{k} c^{2 d} c^{3 k} $$
2
Poly(dict(zip([(1,8),(2,7),(3,9)],[4,5,6])),gens=S('x,y'))
$\displaystyle \operatorname{Poly}{\left( 6 x^{3}y^{9} + 5 x^{2}y^{7} + 4 xy^{8}, x, y, domain=\mathbb{Z} \right)}$
Poly({1:1,2:1,4:1},gens=S('sqrt(x)')).as_expr()
$\displaystyle \sqrt{x} + x^{2} + x$
def _writ2(coef,fun,variables):
	return latex(Poly(dict(zip(fun,coef)),gens=variables).as_expr())
_writ2([7,8],[(1,2),(3,4)],S('x,y'))
'8 x^{3} y^{4} + 7 x y^{2}'
parse_latex(r' a^2+b^2\geq 2ab').rhs
ANTLR runtime and generated code versions disagree: 4.8!=4.7.1
ANTLR runtime and generated code versions disagree: 4.8!=4.7.1
$\displaystyle 2 a b$
parse_latex(r'2ab \geq a^2+b^2').rhs
ANTLR runtime and generated code versions disagree: 4.8!=4.7.1
ANTLR runtime and generated code versions disagree: 4.8!=4.7.1
$\displaystyle a^{2} + b^{2}$
fr=parse_latex(r'\sqrt{ab}+\sqrt{cd}+\sqrt{ef}\leq\sqrt{(a+c+e)(b+d+f)}')
prove(fr.rhs**2-fr.lhs**2)
numerator: $a d + a f + b c + b e + c f + d e - 2 \sqrt{a b} \sqrt{c d} - 2 \sqrt{a b} \sqrt{e f} - 2 \sqrt{c d} \sqrt{e f}$
denominator: $1$
Poly((sqrt(a))**2*(sqrt(d))**2 + (sqrt(a))**2*(sqrt(f))**2 - 2*(sqrt(a))*(sqrt(b))*(sqrt(c))*(sqrt(d)) - 2*(sqrt(a))*(sqrt(b))*(sqrt(e))*(sqrt(f)) + (sqrt(b))**2*(sqrt(c))**2 + (sqrt(b))**2*(sqrt(e))**2 + (sqrt(c))**2*(sqrt(f))**2 - 2*(sqrt(c))*(sqrt(d))*(sqrt(e))*(sqrt(f)) + (sqrt(d))**2*(sqrt(e))**2, sqrt(a), sqrt(b), sqrt(c), sqrt(d), sqrt(e), sqrt(f), domain='ZZ')
Poly((sqrt(a))**2*(sqrt(d))**2 + (sqrt(a))**2*(sqrt(f))**2 + (sqrt(b))**2*(sqrt(c))**2 + (sqrt(b))**2*(sqrt(e))**2 + (sqrt(c))**2*(sqrt(f))**2 + (sqrt(d))**2*(sqrt(e))**2, sqrt(a), sqrt(b), sqrt(c), sqrt(d), sqrt(e), sqrt(f), domain='QQ')
Poly(2*(sqrt(a))*(sqrt(b))*(sqrt(c))*(sqrt(d)) + 2*(sqrt(a))*(sqrt(b))*(sqrt(e))*(sqrt(f)) + 2*(sqrt(c))*(sqrt(d))*(sqrt(e))*(sqrt(f)), sqrt(a), sqrt(b), sqrt(c), sqrt(d), sqrt(e), sqrt(f), domain='QQ')
status: 0
From weighted AM-GM inequality:
$$2 \sqrt{a} \sqrt{b} \sqrt{c} \sqrt{d} \le a d+b c$$
$$2 \sqrt{a} \sqrt{b} \sqrt{e} \sqrt{f} \le a f+b e$$
$$2 \sqrt{c} \sqrt{d} \sqrt{e} \sqrt{f} \le c f+d e$$
$$ 0 \le 0 $$
The sum of all inequalities gives us a proof of the inequality.
0
assumeall(S('sqrt(x*y)-sqrt(x)*sqrt(y)'),positive=True)
$\displaystyle 0$
p=reducegens('x+sqrt(x)')
q=S('x+sqrt(x)')
print(Poly(p,gens=[S('sqrt(x)'),S('sqrt(x)')]))
print(Poly(q,gens=S('sqrt(x)')))
---------------------------------------------------------------------------
GeneratorsError                           Traceback (most recent call last)
<ipython-input-31-b2d05a6575cd> in <module>()
      1 p=reducegens('x+sqrt(x)')
      2 q=S('x+sqrt(x)')
----> 3 print(Poly(p,gens=[S('sqrt(x)'),S('sqrt(x)')]))
      4 print(Poly(q,gens=S('sqrt(x)')))

/home/grzegorz/Pobrane/SageMath/local/lib/python3.7/site-packages/sympy/polys/polytools.py in __new__(cls, rep, *gens, **args)
    107     def __new__(cls, rep, *gens, **args):
    108         """Create a new polynomial instance out of something useful. """
--> 109         opt = options.build_options(gens, args)
    110 
    111         if 'order' in opt:

/home/grzegorz/Pobrane/SageMath/local/lib/python3.7/site-packages/sympy/polys/polyoptions.py in build_options(gens, args)
    729 
    730     if len(args) != 1 or 'opt' not in args or gens:
--> 731         return Options(gens, args)
    732     else:
    733         return args['opt']

/home/grzegorz/Pobrane/SageMath/local/lib/python3.7/site-packages/sympy/polys/polyoptions.py in __init__(self, gens, args, flags, strict)
    152                     self[option] = cls.preprocess(value)
    153 
--> 154         preprocess_options(args)
    155 
    156         for key, value in dict(defaults).items():

/home/grzegorz/Pobrane/SageMath/local/lib/python3.7/site-packages/sympy/polys/polyoptions.py in preprocess_options(args)
    150 
    151                 if value is not None:
--> 152                     self[option] = cls.preprocess(value)
    153 
    154         preprocess_options(args)

/home/grzegorz/Pobrane/SageMath/local/lib/python3.7/site-packages/sympy/polys/polyoptions.py in preprocess(cls, gens)
    289             gens = ()
    290         elif has_dups(gens):
--> 291             raise GeneratorsError("duplicated generators: %s" % str(gens))
    292         elif any(gen.is_commutative is False for gen in gens):
    293             raise GeneratorsError("non-commutative generators: %s" % str(gens))

GeneratorsError: duplicated generators: [sqrt(x), sqrt(x)]
from importlib import reload
from sympy import *
import shiroindev
reload(shiroindev)
from shiroindev import *
from sympy.parsing.latex import parse_latex
from itertools import permutations, combinations
shiro.seed=1
from IPython.display import Latex
shiro.display=lambda x:display(Latex(x))
Poly('sqrt(x)')*Poly('sqrt(x)')
$\displaystyle \operatorname{Poly}{\left( x, \sqrt{x}, domain=\mathbb{Z} \right)}$
Poly('x')-Poly('x').abs()
$\displaystyle \operatorname{Poly}{\left( 0, x, domain=\mathbb{Z} \right)}$
fr=Poly('x-y')
(fr+fr.abs())*(S('1/2'))
$\displaystyle \operatorname{Poly}{\left( x, x, y, domain=\mathbb{Q} \right)}$
from itertools import permutations 
t5=[]
for ineq in ineqs2:
    u=0
    for vars in permutations(ineq.free_symbols):
        ineqp=makesubs(ineq,list(zip(vars[1:],[oo]*len(vars[1:]))),variables=vars[:-1])
        u=max(u,prove(ineqp))
    print(u,end=',')
    t5+=[u]
Counter(t5)
from statsmodels.stats.contingency_tables import mcnemar
print(mcnemar([[101,121],[59,33]]))
pvalue      4.43444926375551e-06
statistic   59.0
nsimplify?
shiro.display=lambda x:display(Latex(x))
newproof()
fr=makesubs(ineqs2[15],'[c,oo],[a,oo]',variables='b,c')
display(fr)
fr2=fractioncancel(fr.subs(S('[[e,(sqrt(5)-1)/2*f]]')))[0]
display(fr2)
prove(makesubs(fr2,'11/5,9/4',variables=[S('sqrt(5)')]))
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
#for all variables in corresponding intervals old formula is nonnegative
	formula=S(formula)
	addsymbols(formula)
	intervals=_smakeiterable2(intervals)
	if variables: variables=_smakeiterable(variables)
	else: variables=sorted(formula.free_symbols,key=str)
	if values!=None:
		values=_smakeiterable(values)
		equations=[var-value for var,value in zip(variables,values)]
	else:
		equations=[]
	newvars=[]
	usedvars=set()
	for var,interval in zip(variables,intervals):
		end1,end2=interval
		z=newvar()
		newvars+=[z]
		usedvars|={z}
		if (end1.free_symbols|end2.free_symbols)&usedvars:
			shiro.warning(shiro.translation[
			'Warning: intervals contain backwards dependencies. Consider changing order of variables'])
		if end1 in {S('-oo'),S('oo')}:
			end1,end2=end2,end1
		if {end1,end2}=={S('-oo'),S('oo')}:
			sub1=sub2=(z-1/z)
		elif end2==S('oo'):
			sub1=sub2=(end1+z)
		elif end2==S('-oo'):
			sub1=sub2=end1-z
		else:
			sub1=end2+(end1-end2)/z
			sub2=end2+(end1-end2)/(1+z)
		formula=formula.subs(var,sub1)
		shiro.display(shiro.translation['Substitute']+" $"+latex(var)+'\\\\to '+latex(sub2)+'$')
		equations=[equation.subs(var,sub1) for equation in equations]
	num,den=fractioncancel(formula)
	for var,interval in zip(newvars,intervals):
		if {interval[0],interval[1]} & {S('oo'),S('-oo')}==set():
			num=num.subs(var,var+1)
			den=den.subs(var,var+1)
			equations=[equation.subs(var,var+1) for equation in equations]
	if values:
		values=ssolve(equations,newvars)
		if len(values):
			values=values[0]
	num,den=expand(num),expand(den)
	#shiro.display(shiro.translation["Formula after substitution:"],"$$",latex(num/den),'$$')
	if values and numden:
		return num,den,values
	elif values:
		return num/den,values
	elif numden:
		return num,den
	else:
		return num/den
x={5,7,8}
x|={3}
x
{3, 5, 7, 8}
from shiroindev import *
S('x+7').free_symbols
{x}
bool(set())
False
from importlib import reload
newproof()
import shiroindev
reload(shiroindev)
from shiroindev import *
shiro.display=lambda x:display(Latex(x))
formula=Sm('xy/z+yz/x+zx/y-1').subs('z',S('1-x-y'))
newformula,values=makesubs(formula,'[0,1-y],[0,1]','1/3,1/3')
prove(newformula,values)
Substitute $x\to - y + 1 + \frac{y - 1}{a + 1}$
Substitute $y\to 1 - \frac{1}{b + 1}$
Substitute $b\to \frac{c}{2}$
numerator: $a^{4} c^{2} + a^{3} c^{2} - 2 a^{3} c - 4 a^{2} c + 4 a^{2} + a c^{2} - 2 a c + c^{2}$
denominator: $a^{3} c^{2} + 2 a^{3} c + 2 a^{2} c^{2} + 4 a^{2} c + a c^{2} + 2 a c$
status: 0
From weighted AM-GM inequality:
$$2 a^{3} c \le a^{4} c^{2}+a^{2}$$
$$4 a^{2} c \le a^{3} c^{2}+2 a^{2}+a c^{2}$$
$$2 a c \le a^{2}+c^{2}$$
$$ 0 \le 0 $$
The sum of all inequalities gives us a proof of the inequality.
0
def findvalues(formula,values=None,variables=None,**kwargs):
	"""finds a candidate for parameter "values" in "prove" function
	blabla"""
	formula=S(formula)
	addsymbols(formula)
	num,den=fractioncancel(formula)
	if variables==None:
		variables=sorted(formula.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,**kwargs))
	return tuple([x*x for x in tup])
findvalues(Sm('(x^2-4x+4)/y'))
Optimization terminated successfully.
         Current function value: 1.000000
         Iterations: 51
         Function evaluations: 109
(1.9999999989778676, 0.6689493176765962)
findvalues?