def test_squarefree():
    assert Poly(x-1, x).is_squarefree == True
    assert Poly((x-1)**2, x).is_squarefree == False

    assert Poly(3*x**2, x).as_squarefree() == Poly(3*x, x)
    assert Poly(x**2+2*x+1, x).as_squarefree() == Poly(x+1, x)
    assert Poly(x**5-x**4-x+1, x).as_squarefree() == Poly(x**4-1, x)

    assert poly_sqf(1, x) == [Poly(1, x)]
    assert poly_sqf(x, x) == [Poly(x, x)]

    assert poly_sqf(3*x**2, x) == [Poly(3, x), Poly(x, x)]
    assert poly_sqf(x**2+2*x+1, x) == [Poly(1, x), Poly(x+1, x)]

    assert poly_sqf(x**5-x**4-x+1, x) == \
        [Poly(x**3 + x**2 + x + 1, x), Poly(x-1, x)]
    assert poly_sqf(x**8+6*x**6+12*x**4+8*x**2, x) == \
        [Poly(1, x), Poly(x, x), Poly(x**2+2, x)]

    # Bronstein, Symbolic Integration, pp. 52

    A = Poly(x**4 - 3*x**2 + 6, x)
    D = Poly(x**6 - 5*x**4 + 5*x**2 + 4, x)

    f, g = D, A - D.diff(x).mul_term(t)

    res, R = poly_subresultants(f, g)
    S = poly_sqf(Poly(res, t))

    assert S == [Poly(45796, t), Poly(1, t), Poly(4*t**2 + 1, t)]

def gauss_lobatto_points(p):
    """
    Returns the list of Gauss-Lobatto points of the order 'p'.
    """
    x = Symbol("x")
    print "creating"
    e = legendre(p, x).diff(x)
    e = Poly(e, x)
    print "polydone"
    if e == 0:
        return []
    print "roots"
    if e == 1:
        r = []
    else:
        with workdps(40):
            r, err = polyroots(e.all_coeffs(), error=True)
        #if err > 1e-40:
        #    raise Exception("Internal Error: Root is not precise")
    print "done"
    p = []
    p.append("-1.0")
    for x in r:
        if abs(x) < 1e-40: x = 0
        p.append(str(x))
    p.append("1.0")
    return p

def main():
    x=Symbol("x")
    s = Poly(A(x), x)
    num = list(reversed(s.coeffs()))[:11]

    print s.as_basic()
    print num

    def _eval_expand_func(self, **hints):
        from sympy import exp, I, floor, Add, Poly, Dummy, exp_polar, unpolarify
        z, s, a = self.args
        if z == 1:
            return zeta(s, a)
        if s.is_Integer and s <= 0:
            t = Dummy('t')
            p = Poly((t + a)**(-s), t)
            start = 1/(1 - t)
            res = S(0)
            for c in reversed(p.all_coeffs()):
                res += c*start
                start = t*start.diff(t)
            return res.subs(t, z)

        if a.is_Rational:
            # See section 18 of
            #   Kelly B. Roach.  Hypergeometric Function Representations.
            #   In: Proceedings of the 1997 International Symposium on Symbolic and
            #   Algebraic Computation, pages 205-211, New York, 1997. ACM.
            # TODO should something be polarified here?
            add = S(0)
            mul = S(1)
            # First reduce a to the interaval (0, 1]
            if a > 1:
                n = floor(a)
                if n == a:
                    n -= 1
                a -= n
                mul = z**(-n)
                add = Add(*[-z**(k - n)/(a + k)**s for k in xrange(n)])
            elif a <= 0:
                n = floor(-a) + 1
                a += n
                mul = z**n
                add = Add(*[z**(n - 1 - k)/(a - k - 1)**s for k in xrange(n)])

            m, n = S([a.p, a.q])
            zet = exp_polar(2*pi*I/n)
            root = z**(1/n)
            return add + mul*n**(s - 1)*Add(
                *[polylog(s, zet**k*root)._eval_expand_func(**hints)
                  / (unpolarify(zet)**k*root)**m for k in xrange(n)])

        # TODO use minpoly instead of ad-hoc methods when issue 2789 is fixed
        if z.func is exp and (z.args[0]/(pi*I)).is_Rational or z in [-1, I, -I]:
            # TODO reference?
            if z == -1:
                p, q = S([1, 2])
            elif z == I:
                p, q = S([1, 4])
            elif z == -I:
                p, q = S([-1, 4])
            else:
                arg = z.args[0]/(2*pi*I)
                p, q = S([arg.p, arg.q])
            return Add(*[exp(2*pi*I*k*p/q)/q**s*zeta(s, (k + a)/q)
                         for k in xrange(q)])

        return lerchphi(z, s, a)

def main():
    x=Symbol("x")
    s = Poly(A(x), x)
    num = [s.coeff(n) for n in range(11)]

    print s.as_basic()
    print num

def test_coeff():
    p = Poly(3*x**2*y + 4*x*y**2 + 1, x, y)

    assert p.coeff(0, 0) == p.coeff() == 1

    assert p.coeff(2, 1) == 3
    assert p.coeff(1, 2) == 4
    assert p.coeff(1, 1) == 0

 def linear_arg(arg):
     """ Test if arg is of form a*s+b, raise exception if not. """
     if not arg.is_polynomial(s):
         raise exception(fact)
     p = Poly(arg, s)
     if p.degree() != 1:
         raise exception(fact)
     return p.all_coeffs()

def test_has_any():
    x,y,z,t,u = symbols('x y z t u')
    f = Function("f")
    g = Function("g")
    p = Wild('p')

    assert sin(x).has(x)
    assert sin(x).has(sin)
    assert not sin(x).has(y)
    assert not sin(x).has(cos)
    assert f(x).has(x)
    assert f(x).has(f)
    assert not f(x).has(y)
    assert not f(x).has(g)

    assert f(x).diff(x).has(x)
    assert f(x).diff(x).has(f)
    assert f(x).diff(x).has(Derivative)
    assert not f(x).diff(x).has(y)
    assert not f(x).diff(x).has(g)
    assert not f(x).diff(x).has(sin)

    assert (x**2).has(Symbol)
    assert not (x**2).has(Wild)
    assert (2*p).has(Wild)

    i = Integer(4400)

    assert i.has(x) is False

    assert (i*x**i).has(x)
    assert (i*y**i).has(x) is False
    assert (i*y**i).has(x, y)

    expr = x**2*y + sin(2**t + log(z))

    assert expr.has(u) is False

    assert expr.has(x)
    assert expr.has(y)
    assert expr.has(z)
    assert expr.has(t)

    assert expr.has(x, y, z, t)
    assert expr.has(x, y, z, t, u)

    from sympy.physics.units import m, s

    assert (x*m/s).has(x)
    assert (x*m/s).has(y, z) is False

    poly = Poly(x**2 + x*y*sin(z), x, y, t)

    assert poly.has(x)
    assert poly.has(x, y, z)
    assert poly.has(x, y, z, t)

    assert FockState((x, y)).has(x)

def test_add_sub_term():
    f = Poly((), x, y, z)

    assert Poly(2, x) - Poly(1, x) == Poly(1, x)
    assert Poly(1, x) - Poly(1, x) == Poly(0, x)
    assert Poly(1, x) - Poly(2, x) == Poly(-1, x)

    assert f.add_term(12, (1,2,3)) == Poly(((12,), ((1,2,3),)), x,y,z)
    assert f.sub_term(12, (1,2,3)) == Poly(((-12,), ((1,2,3),)), x,y,z)

	def LM(self,poly):
		from sympy import S,Poly
		if poly==0:return poly
		poly = Poly(poly,self.variables)
		polyDict = poly.as_dict()
		exponents = polyDict.keys()
		largest = exponents[0]
		for i in xrange(len(polyDict)):
			if self.compare(exponents[i],largest): largest = exponents[i]
		return Poly({largest:S(1)},self.variables)

def poly_factorize(poly):
    """Factorize multivariate polynomials into a sum of products of monomials.

    This function can be used to decompose polynomials into a form which
    minimizes the number of additions and multiplications, and which thus
    can be evaluated efficently."""
    max_deg = {}

    if 'horner' in dir(sympy):
        return sympy.horner(poly)

    if not isinstance(poly, Poly):
        poly = Poly(sympy.expand(poly), *poly.atoms(Symbol))

    denom, poly = poly.as_integer()

    # Determine the order of factorization.  We proceed through the
    # symbols, starting with the one present in the highest order
    # in the polynomial.
    for i, sym in enumerate(poly.symbols):
        max_deg[i] = 0

    for monom in poly.monoms:
        for i, symvar in enumerate(monom):
            max_deg[i] = max(max_deg[i], symvar)

    ret_poly = 0
    monoms = list(poly.monoms)

    for isym, maxdeg in sorted(max_deg.items(), key=itemgetter(1), reverse=True):
        drop_idx = []
        new_poly = []

        for i, monom in enumerate(monoms):
            if monom[isym] > 0:
                drop_idx.append(i)
                new_poly.append((poly.coeff(*monom), monom))

        if not new_poly:
            continue

        ret_poly += sympy.factor(Poly(new_poly, *poly.symbols))

        for idx in reversed(drop_idx):
            del monoms[idx]

    # Add any remaining O(1) terms.
    new_poly = []
    for i, monom in enumerate(monoms):
        new_poly.append((poly.coeff(*monom), monom))

    if new_poly:
        ret_poly += Poly(new_poly, *poly.symbols)

    return ret_poly / denom

def test_as_poly_as_expr():
    f = x ** 2 + 2 * x * y

    assert f.as_poly().as_expr() == f
    assert f.as_poly(x, y).as_expr() == f

    assert (f + sin(x)).as_poly(x, y) is None

    p = Poly(f, x, y)

    assert p.as_poly() == p

def test_iterators():
    f = Poly(x**5 + x, x)

    assert list(f.iter_all_coeffs()) == \
        [1, 0, 0, 0, 1, 0]

    assert list(f.iter_all_monoms()) == \
        [(5,), (4,), (3,), (2,), (1,), (0,)]

    assert list(f.iter_all_terms()) == \
        [(1, (5,)), (0, (4,)), (0, (3,)), (0, (2,)), (1, (1,)), (0, (0,))]

def test_map_coeffs():
    p = Poly(x**2 + 2*x*y, x, y)
    q = p.map_coeffs(lambda c: 2*c)

    assert q.as_basic() == 2*x**2 + 4*x*y

    p = Poly(u*x**2 + v*x*y, x, y)
    q = p.map_coeffs(expand, complex=True)

    assert q.as_basic() == x**2*(I*im(u) + re(u)) + x*y*(I*im(v) + re(v))

    raises(PolynomialError, "p.map_coeffs(lambda c: x*c)")

def test_as_poly_basic():
    x, y = symbols('xy')

    f = x**2 + 2*x*y

    assert f.as_poly().as_basic() == f
    assert f.as_poly(x, y).as_basic() == f

    assert (f + sin(x)).as_poly(x, y) is None

    p = Poly(f, x, y)

    assert p.as_poly() == p

	def reduce(self, debug=False, inputFilename=None, outputFilename=None):
		previousLeftBlockRegion = -1
		previousRightBlockRegion = -1
		equationCount = 0
		partIndex = 0
		equationInPartIndex = 0

		output = None
		if debug:
			input = open(inputFilename, 'r')
			output = open(outputFilename, 'w')
		else:
			input = sys.stdin

		for strEquation in input:
			equationCount += 1
			index, strEquation = strEquation.strip().split('\t')
			equation, symbols = decodeExpression(strEquation)
			equation = Poly(equation, symbols)
			equationCoefficients = [0] * len(symbols)
			for monom, coeff in equation.terms():
				try:
					i = list(monom).index(1)
					equationCoefficients[i] = 1 if coeff != 0 else 0
				except ValueError:
					continue

			equationLen = len(equationCoefficients)
			rightColumnIndex = equationLen - 1
			leftColumnIndex = 0

			while leftColumnIndex < equationLen and equationCoefficients[leftColumnIndex] != 1:
				leftColumnIndex += 1
			while rightColumnIndex > 0 and equationCoefficients[rightColumnIndex] != 1:
				rightColumnIndex -= 1

			if previousLeftBlockRegion != leftColumnIndex or previousRightBlockRegion != rightColumnIndex:
				if equationInPartIndex > 0:
					partIndex += 1
					equationInPartIndex = 0
				previousLeftBlockRegion = leftColumnIndex
				previousRightBlockRegion = rightColumnIndex
			if output is not None:
				output.write('{0}\t{1}\t{2}\n'.format(str(partIndex), index, strEquation))
			else:
				print('{0}\t{1}\t{2}'.format(str(partIndex), index, strEquation))
			equationInPartIndex += 1
		if debug:
			input.close()
		if output is not None:
			output.close()

def _simplify(expr):
    u"""Simplifie une expression.

    Alias de simplify (sympy 0.6.4).
    Mais simplify n'est pas garanti d'être stable dans le temps.
    (Cf. simplify.__doc__)."""
    return together(expand(Poly.cancel(powsimp(expr))))

def gauss_lobatto_points(p):
    """
    Returns the list of Gauss-Lobatto points of the order 'p'.
    """
    x = Symbol("x")
    e = (1-x**2)*legendre(p, x).diff(x)
    e = Poly(e, x)
    if e == 0:
        return []
    with workdps(40):
        r, err = polyroots(e.all_coeffs(), error=True)
    if err > 1e-40:
        raise Exception("Internal Error: Root is not precise")
    p = []
    for x in r:
        if abs(x) < 1e-40: x = 0
        p.append(str(x))
    return p

def test_mul_div_term():
    f = Poly(x*y**2 + 2*y, x, y)

    assert f.mul_term(0, (0, 1)) == Poly(0, x, y)

    assert f.mul_term(1) == Poly(x*y**2 + 2*y, x, y)
    assert f.mul_term(1, (0, 1)) == Poly(x*y**3 + 2*y**2, x, y)

    raises(ZeroDivisionError, "f.div_term(0)")

    assert f.div_term(1) == Poly(x*y**2 + 2*y, x, y)
    assert f.div_term(1, (0, 1)) == Poly(x*y + 2, x, y)

    assert f.mul_term(2) == Poly(2*x*y**2 + 4*y, x, y)
    assert f.mul_term(2, (0, 1)) == Poly(2*x*y**3 + 4*y**2, x, y)

    assert f.div_term(2) == Poly(x*y**2 / 2 + y, x, y)
    assert f.div_term(2, (0, 1)) == Poly(x*y / 2 + 1, x, y)

def new_poly(V_s, fields, n_coeffs):
    """Make a new polynomial function that has the same powers as V_s function
       but with coefficients C1, C2..."""
    
    from sympy import Poly
    from sympy import diff, Symbol, var, simplify, sympify, S
    from sympy.core.sympify import SympifyError
    
    P = Poly(V_s,*fields)
    d = P.as_dict()
    e = {}
    
    for key in d.iterkeys():
        #print d[key], str(d[key])
        for i in xrange(1, n_coeffs+1):
            if 'C' + str(i) in str(d[key]):
                e[key] = sympify('C'+str(i))

    P2 = Poly(e,*fields)
    return str(P2.as_basic())