def test_diop_general_sum_of_squares_quick():
    for i in range(3, 10):
        assert check_solutions(sum(i**2 for i in symbols(':%i' % i)) - i)
    raises(ValueError, lambda: _diop_general_sum_of_squares((x, y), 2))
    assert _diop_general_sum_of_squares((x, y, z), -2) == set()
    eq = x**2 + y**2 + z**2 - (1 + 4 + 9)
    assert diop_general_sum_of_squares(eq) == \
        set([(1, 2, 3)])
    eq = u**2 + v**2 + x**2 + y**2 + z**2 - 1313
    assert len(diop_general_sum_of_squares(eq, 3)) == 3
    # issue 11016
    var = symbols(':5') + (symbols('6', negative=True),)
    eq = Add(*[i**2 for i in var]) - 112
    assert diophantine(eq) == set(
        [(0, 1, 1, 5, 6, -7), (1, 1, 1, 3, 6, -8), (2, 3, 3, 4,
        5, -7), (0, 1, 1, 1, 3, -10), (0, 0, 4, 4, 4, -8), (1, 2, 3,
        3, 5, -8), (0, 1, 2, 3, 7, -7), (2, 2, 4, 4, 6, -6), (1, 1,
        3, 4, 6, -7), (0, 2, 3, 3, 3, -9), (0, 0, 2, 2, 2, -10), (1,
        1, 2, 3, 4, -9), (0, 1, 1, 2, 5, -9), (0, 0, 2, 6, 6, -6),
        (1, 3, 4, 5, 5, -6), (0, 2, 2, 2, 6, -8), (0, 3, 3, 3, 6,
        -7), (0, 2, 3, 5, 5, -7), (0, 1, 5, 5, 5, -6)])
    # handle negated squares with signsimp
    assert diophantine(12 - x**2 - y**2 - z**2) == set([(2, 2, 2)])
    # diophantine handles simplification, so classify_diop should
    # not have to look for additional patterns that are removed
    # by diophantine
    eq = a**2 + b**2 + c**2 + d**2 - 4
    raises(NotImplementedError, lambda: classify_diop(-eq))

    def _intersect(self, other):
        from sympy import Dummy
        from sympy.solvers.diophantine import diophantine
        from sympy.sets.sets import imageset
        if self.base_set is S.Integers:
            if isinstance(other, ImageSet) and other.base_set is S.Integers:
                f, g = self.lamda.expr, other.lamda.expr
                n, m = self.lamda.variables[0], other.lamda.variables[0]

                # Diophantine sorts the solutions according to the alphabetic
                # order of the variable names, since the result should not depend
                # on the variable name, they are replaced by the dummy variables
                # below
                a, b = Dummy('a'), Dummy('b')
                f, g = f.subs(n, a), g.subs(m, b)
                solns_set = diophantine(f - g)
                if solns_set == set():
                    return EmptySet()
                solns = list(diophantine(f - g))
                if len(solns) == 1:
                    t = list(solns[0][0].free_symbols)[0]
                else:
                    return None

                # since 'a' < 'b'
                return imageset(Lambda(t, f.subs(a, solns[0][0])), S.Integers)

def test_assumptions():
    """
    Test whether diophantine respects the assumptions.
    """
    #Test case taken from the below so question regarding assumptions in diophantine module
    #http://stackoverflow.com/questions/23301941/how-can-i-declare-natural-symbols-with-sympy
    m, n = symbols('m n', integer=True, positive=True)
    diof = diophantine(n ** 2 + m * n - 500)
    assert diof == set([(5, 20), (40, 10), (95, 5), (121, 4), (248, 2), (499, 1)])

    a, b = symbols('a b', integer=True, positive=False)
    diof = diophantine(a*b + 2*a + 3*b - 6)
    assert diof == set([(-15, -3), (-9, -4), (-7, -5), (-6, -6), (-5, -8), (-4, -14)])

def test_diophantine_permute_sign():
    from sympy.abc import a, b, c, d, e
    eq = a**4 + b**4 - (2**4 + 3**4)
    base_sol = set([(2, 3)])
    assert diophantine(eq) == base_sol
    complete_soln = set(signed_permutations(base_sol.pop()))
    assert diophantine(eq, permute=True) == complete_soln

    eq = a**2 + b**2 + c**2 + d**2 + e**2 - 234
    assert len(diophantine(eq)) == 35
    assert len(diophantine(eq, permute=True)) == 62000
    soln = set([(-1, -1), (-1, 2), (1, -2), (1, 1)])
    assert diophantine(10*x**2 + 12*x*y + 12*y**2 - 34, permute=True) == soln

    def _intersect(self, other):
        from sympy.solvers.diophantine import diophantine
        if self.base_set is S.Integers:
            g = None
            if isinstance(other, ImageSet) and other.base_set is S.Integers:
                g = other.lamda.expr
                m = other.lamda.variables[0]
            elif other is S.Integers:
                m = g = Dummy('x')
            if g is not None:
                f = self.lamda.expr
                n = self.lamda.variables[0]
                # Diophantine sorts the solutions according to the alphabetic
                # order of the variable names, since the result should not depend
                # on the variable name, they are replaced by the dummy variables
                # below
                a, b = Dummy('a'), Dummy('b')
                f, g = f.subs(n, a), g.subs(m, b)
                solns_set = diophantine(f - g)
                if solns_set == set():
                    return EmptySet()
                solns = list(diophantine(f - g))

                if len(solns) != 1:
                    return

                # since 'a' < 'b', select soln for n
                nsol = solns[0][0]
                t = nsol.free_symbols.pop()
                return imageset(Lambda(n, f.subs(a, nsol.subs(t, n))), S.Integers)

        if other == S.Reals:
            from sympy.solvers.solveset import solveset_real
            from sympy.core.function import expand_complex
            if len(self.lamda.variables) > 1:
                return None

            f = self.lamda.expr
            n = self.lamda.variables[0]

            n_ = Dummy(n.name, real=True)
            f_ = f.subs(n, n_)

            re, im = f_.as_real_imag()
            im = expand_complex(im)

            return imageset(Lambda(n_, re),
                            self.base_set.intersect(
                                solveset_real(im, n_)))

def check_solutions(eq):
    """
    Determines whether solutions returned by diophantine() satisfy the original
    equation. Hope to generalize this so we can remove functions like check_ternay_quadratic,
    check_solutions_normal, check_solutions()
    """
    s = diophantine(eq)

    terms = factor_list(eq)[1]

    var = list(eq.free_symbols)
    var.sort(key=default_sort_key)

    okay = True

    while len(s) and okay:
        solution = s.pop()

        okay = False

        for term in terms:
            subeq = term[0]

            if simplify(_mexpand(Subs(subeq, var, solution).doit())) == 0:
                okay = True
                break

    return okay

    def _intersect(self, other):
        from sympy import Dummy
        from sympy.solvers.diophantine import diophantine
        from sympy.sets.sets import imageset

        if self.base_set is S.Integers:
            if isinstance(other, ImageSet) and other.base_set is S.Integers:
                f, g = self.lamda.expr, other.lamda.expr
                n, m = self.lamda.variables[0], other.lamda.variables[0]

                # Diophantine sorts the solutions according to the alphabetic
                # order of the variable names, since the result should not depend
                # on the variable name, they are replaced by the dummy variables
                # below
                a, b = Dummy("a"), Dummy("b")
                f, g = f.subs(n, a), g.subs(m, b)
                solns_set = diophantine(f - g)
                if solns_set == set():
                    return EmptySet()
                solns = list(diophantine(f - g))
                if len(solns) == 1:
                    t = list(solns[0][0].free_symbols)[0]
                else:
                    return None

                # since 'a' < 'b'
                return imageset(Lambda(t, f.subs(a, solns[0][0])), S.Integers)

        if other == S.Reals:
            from sympy.solvers.solveset import solveset_real
            from sympy.core.function import expand_complex

            if len(self.lamda.variables) > 1:
                return None

            f = self.lamda.expr
            n = self.lamda.variables[0]

            n_ = Dummy(n.name, real=True)
            f_ = f.subs(n, n_)

            re, im = f_.as_real_imag()
            im = expand_complex(im)

            return imageset(Lambda(n_, re), self.base_set.intersect(solveset_real(im, n_)))

def test_diopcoverage():
    eq = (2*x + y + 1)**2
    assert diop_solve(eq) == set([(t_0, -2*t_0 - 1)])
    eq = 2*x**2 + 6*x*y + 12*x + 4*y**2 + 18*y + 18
    assert diop_solve(eq) == set([(t_0, -t_0 - 3), (2*t_0 - 3, -t_0)])
    assert diop_quadratic(x + y**2 - 3) == set([(-t**2 + 3, -t)])

    assert diop_linear(x + y - 3) == (t_0, 3 - t_0)

    assert base_solution_linear(0, 1, 2, t=None) == (0, 0)
    ans = (3*t - 1, -2*t + 1)
    assert base_solution_linear(4, 8, 12, t) == ans
    assert base_solution_linear(4, 8, 12, t=None) == tuple(_.subs(t, 0) for _ in ans)

    assert cornacchia(1, 1, 20) is None
    assert cornacchia(1, 1, 5) == set([(2, 1)])
    assert cornacchia(1, 2, 17) == set([(3, 2)])

    raises(ValueError, lambda: reconstruct(4, 20, 1))

    assert gaussian_reduce(4, 1, 3) == (1, 1)
    eq = -w**2 - x**2 - y**2 + z**2

    assert diop_general_pythagorean(eq) == \
        diop_general_pythagorean(-eq) == \
            (m1**2 + m2**2 - m3**2, 2*m1*m3,
            2*m2*m3, m1**2 + m2**2 + m3**2)

    assert check_param(S(3) + x/3, S(4) + x/2, S(2), x) == (None, None)
    assert check_param(S(3)/2, S(4) + x, S(2), x) == (None, None)
    assert check_param(S(4) + x, S(3)/2, S(2), x) == (None, None)

    assert _nint_or_floor(16, 10) == 2
    assert _odd(1) == (not _even(1)) == True
    assert _odd(0) == (not _even(0)) == False
    assert _remove_gcd(2, 4, 6) == (1, 2, 3)
    raises(TypeError, lambda: _remove_gcd((2, 4, 6)))
    assert sqf_normal(2 * 3**2 * 5, 2 * 5 * 11, 2 * 7**2 * 11)  == \
        (11, 1, 5)

    # it's ok if these pass some day when the solvers are implemented
    raises(NotImplementedError, lambda: diophantine(x**2 + y**2 + x*y + 2*y*z - 12))
    raises(NotImplementedError, lambda: diophantine(x**3 + y**2))
    assert diop_quadratic(x**2 + y**2 - 1**2 - 3**4) == \
        set([(-9, -1), (-9, 1), (-1, -9), (-1, 9), (1, -9), (1, 9), (9, -1), (9, 1)])

def _has_conflict(rd1, rd2, offset):
    """
    :param rd1: rectangular domain from Stencil1 (1D) <- (low, high, stride) (written domain)
    :param rd2: rectangular domain from Stencil2 (1D) <- (low, high, stride) (read domain)
    :param offset_vector: Offset from the center (relative to rd2), 1D projections
    :return: Whether OV from rectangular domain 2 reads from rectangular domain 1
    """

    n1 = sympy.Symbol("n1")
    n2 = sympy.Symbol("n2")

    # Diophantine equations:
    # offset(iterationspace1) + n1 * stride(iteration_space1) <- write vectors
    diophantine_eq1 = rd1[0] + n1 * rd1[2]
    # offset(iterationspace2) + n2 * stride(iteration_space2) + offset_vector  <- Read vectors
    diophantine_eq2 = rd2[0] + n2 * rd2[2] + offset
    # Since sympy solves eq = 0, we want dio_eq1 - dio_eq2 = 0.
    eqn = diophantine_eq1 - diophantine_eq2
    parameter = sympy.Symbol("t", integer=True)
    sat_param = sympy.Symbol("t_0", integer=True)
    if not eqn.free_symbols:
        return eqn == 0
    solutions = dio.diophantine(eqn, parameter)  # default parameter is "t"
    # print("Sols:", solutions)
    for sol in solutions:
        if len(eqn.free_symbols) != 2:
            if n1 in eqn.free_symbols:
                parametric_n1 = sol[0]
                parametric_n2 = 0
            else:
                parametric_n1 = 0
                parametric_n2 = sol[0]
        else:
            (parametric_n1, parametric_n2) = sol
        # Solutions is a set of tuples, each containing either a number or parametric expression
        # which give conditions on satisfiability.

        # are these satisfiable on the original bounds?
        # substitute the parametric forms in
        substituted_1 = diophantine_eq1.subs({n1: parametric_n1})
        substituted_2 = (
            rd2[0] + parametric_n2 * rd2[2]
        )  # we ditch the offset because it's irrelevant since the bounds are based on the center of the stencil

        # print(substituted_1, "\t", substituted_2)
        # print(rd1[0], rd1[1], rd2[0], rd2[1])
        # now do they satisfy the bounds?
        satisfactory_interval = ineq.reduce_rational_inequalities(
            [[rd1[0] <= substituted_1, rd1[1] > substituted_1, rd2[0] <= substituted_2, rd2[1] > substituted_2]],
            sat_param,
            relational=False,
        )
        # print(satisfactory_interval)

        if _contains_integer(satisfactory_interval):
            return True
    return False

def analyze_dependencies(vec, space1, space2):
    x, y = sympy.symbols("x y", integer=True)  # used for setting up diophantine equations
    # TODO: Make it work on superficially colliding domains (but in reality spatially separated)

    results = []
    for dim, (s1, o1, v_comp, s2, o2) in enumerate(zip(space1.stride, space1.lower, vec, space2.stride, space2.lower)):
        if isinstance(v_comp, sympy.Expr):
            v_comp = v_comp.subs({"index_{}".format(dim): s1 * x})
        results.append(dio.diophantine(s1 * x + o1 + v_comp - s2 * y - o2))
    return results

def test_quadratic_simple_hyperbolic_case():
    # Simple Hyperbolic case: A = C = 0 and B != 0
    assert diop_solve(3*x*y + 34*x - 12*y + 1) == \
        set([(-Integer(133), -Integer(11)), (Integer(5), -Integer(57))])
    assert diop_solve(6*x*y + 2*x + 3*y + 1) == set([])
    assert diop_solve(-13*x*y + 2*x - 4*y - 54) == set([(Integer(27), Integer(0))])
    assert diop_solve(-27*x*y - 30*x - 12*y - 54) == set([(-Integer(14), -Integer(1))])
    assert diop_solve(2*x*y + 5*x + 56*y + 7) == set([(-Integer(161), -Integer(3)),\
        (-Integer(47),-Integer(6)), (-Integer(35), -Integer(12)), (-Integer(29), -Integer(69)),\
        (-Integer(27), Integer(64)), (-Integer(21), Integer(7)),(-Integer(9), Integer(1)),\
        (Integer(105), -Integer(2))])
    assert diop_solve(6*x*y + 9*x + 2*y + 3) == set([])
    assert diop_solve(x*y + x + y + 1) == set([(-Integer(1), t), (t, -Integer(1))])
    assert diophantine(48*x*y)

def check_integrality(eq):
    """
    Check that the solutions returned by diophantine() are integers.
    This should be seldom needed except for general quadratic
    equations which are solved with rational transformations.
    """
    def _check_values(x):
        """ Check a number of values. """
        for i in range(-4, 4):
            if not isinstance(simplify(x.subs(t, i)), Integer):
                return False
        return True

    for soln in diophantine(eq, param=t):
        for x in soln:
            if not _check_values(x):
                return False

    return True

def check_solutions(eq):
    """
    Determines whether solutions returned by diophantine() satisfy the original
    equation. Hope to generalize this so we can remove functions like check_ternay_quadratic,
    check_solutions_normal, check_solutions()
    """
    s = diophantine(eq)

    factors = Mul.make_args(eq)

    var = list(eq.free_symbols)
    var.sort(key=default_sort_key)

    while s:
        solution = s.pop()
        for f in factors:
            if diop_simplify(f.subs(zip(var, solution))) == 0:
                break
        else:
            return False
    return True

def p100():
	eq = x**2 - y**2 - 2*x*y - x + y
	n = symbols('n',integer=True)
	s = diophantine(eq,n)
	x_1,y_1 = s.pop()
	x_2,y_2 = s.pop()
	x_3,y_3 = s.pop()
	x_4,y_4 = s.pop()
	xxs = (x_1, x_2, x_3, x_4)
	yys = (y_1,y_2,y_3,y_4)
	print 'n eqnum x y '
	for i in xrange(9):
		for e,(xx,yy) in enumerate(zip(xxs,yys)):
			try:
				x_val = simplify(xx.subs({n:i}))
				if x_val > 0:
					y_val = simplify(yy.subs({n:i}))
					if x_val+y_val>10**12:
						print 'DONE', i, e, x_val, y_val, x_val+y_val, x_val+y_val>10**12
						return x_val
					print i, e, x_val, y_val, x_val+y_val, x_val+y_val>10**12
			except KeyError:
				print 'strange sympy error', i, e

def test_not_implemented():
    eq = x**2 + y**4 - 1**2 - 3**4
    assert diophantine(eq, syms=[x, y]) == set([(9, 1), (1, 3)])

def test_issue_9106():
    eq = -48 - 2*x*(3*x - 1) + y*(3*y - 1)
    v = (x, y)
    for sol in diophantine(eq):
        assert not diop_simplify(eq.xreplace(dict(zip(v, sol))))

def test_input_format():
    raises(TypeError, lambda: diophantine(sin(x)))

    def _intersect(self, other):
        from sympy.solvers.diophantine import diophantine
        if self.base_set is S.Integers:
            g = None
            if isinstance(other, ImageSet) and other.base_set is S.Integers:
                g = other.lamda.expr
                m = other.lamda.variables[0]
            elif other is S.Integers:
                m = g = Dummy('x')
            if g is not None:
                f = self.lamda.expr
                n = self.lamda.variables[0]
                # Diophantine sorts the solutions according to the alphabetic
                # order of the variable names, since the result should not depend
                # on the variable name, they are replaced by the dummy variables
                # below
                a, b = Dummy('a'), Dummy('b')
                f, g = f.subs(n, a), g.subs(m, b)
                solns_set = diophantine(f - g)
                if solns_set == set():
                    return EmptySet()
                solns = list(diophantine(f - g))

                if len(solns) != 1:
                    return

                # since 'a' < 'b', select soln for n
                nsol = solns[0][0]
                t = nsol.free_symbols.pop()
                return imageset(Lambda(n, f.subs(a, nsol.subs(t, n))), S.Integers)

        if other == S.Reals:
            from sympy.solvers.solveset import solveset_real
            from sympy.core.function import expand_complex
            if len(self.lamda.variables) > 1:
                return None

            f = self.lamda.expr
            n = self.lamda.variables[0]

            n_ = Dummy(n.name, real=True)
            f_ = f.subs(n, n_)

            re, im = f_.as_real_imag()
            im = expand_complex(im)

            return imageset(Lambda(n_, re),
                            self.base_set.intersect(
                                solveset_real(im, n_)))

        elif isinstance(other, Interval):
            from sympy.solvers.solveset import (invert_real, invert_complex,
                                                solveset)

            f = self.lamda.expr
            n = self.lamda.variables[0]
            base_set = self.base_set
            new_inf, new_sup = None, None
            new_lopen, new_ropen = other.left_open, other.right_open

            if f.is_real:
                inverter = invert_real
            else:
                inverter = invert_complex

            g1, h1 = inverter(f, other.inf, n)
            g2, h2 = inverter(f, other.sup, n)

            if all(isinstance(i, FiniteSet) for i in (h1, h2)):
                if g1 == n:
                    if len(h1) == 1:
                        new_inf = h1.args[0]
                if g2 == n:
                    if len(h2) == 1:
                        new_sup = h2.args[0]
                # TODO: Design a technique to handle multiple-inverse
                # functions

                # Any of the new boundary values cannot be determined
                if any(i is None for i in (new_sup, new_inf)):
                    return

                range_set = S.EmptySet

                if all(i.is_real for i in (new_sup, new_inf)):
                    new_interval = Interval(new_inf, new_sup, new_lopen, new_ropen)
                    range_set = base_set._intersect(new_interval)
                else:
                    if other.is_subset(S.Reals):
                        solutions = solveset(f, n, S.Reals)
                        if not isinstance(range_set, (ImageSet, ConditionSet)):
                            range_set = solutions._intersect(other)
                        else:
                            return

                if range_set is S.EmptySet:
                    return S.EmptySet
                elif isinstance(range_set, Range) and range_set.size is not S.Infinity:
                    range_set = FiniteSet(*list(range_set))

#.........这里部分代码省略.........

def test_diophantine():
    assert check_solutions((x - y)*(y - z)*(z - x))
    assert check_solutions((x - y)*(x**2 + y**2 - z**2))
    assert check_solutions((x - 3*y + 7*z)*(x**2 + y**2 - z**2))
    assert check_solutions((x**2 - 3*y**2 - 1))
    assert check_solutions(y**2 + 7*x*y)
    assert check_solutions(x**2 - 3*x*y + y**2)
    assert check_solutions(z*(x**2 - y**2 - 15))
    assert check_solutions(x*(2*y - 2*z + 5))
    assert check_solutions((x**2 - 3*y**2 - 1)*(x**2 - y**2 - 15))
    assert check_solutions((x**2 - 3*y**2 - 1)*(y - 7*z))
    assert check_solutions((x**2 + y**2 - z**2)*(x - 7*y - 3*z + 4*w))
    # Following test case caused problems in parametric representation
    # But this can be solved by factroing out y.
    # No need to use methods for ternary quadratic equations.
    assert check_solutions(y**2 - 7*x*y + 4*y*z)
    assert check_solutions(x**2 - 2*x + 1)

    assert diophantine(x - y) == diophantine(Eq(x, y))
    assert diophantine(3*x*pi - 2*y*pi) == set([(2*t_0, 3*t_0)])
    assert diophantine(x**2 + y**2 + z**2 - 14) == set([(1, 2, 3)])
    assert diophantine(x**2 + 15*x/14 - 3) == set()
    # test issue 11049
    eq = 92*x**2 - 99*y**2 - z**2
    coeff = eq.as_coefficients_dict()
    assert _diop_ternary_quadratic_normal((x, y, z), coeff) == \
        (9, 7, 51)
    assert diophantine(eq) == set([(
        891*p**2 + 9*q**2, -693*p**2 - 102*p*q + 7*q**2,
        5049*p**2 - 1386*p*q - 51*q**2)])
    eq = 2*x**2 + 2*y**2 - z**2
    coeff = eq.as_coefficients_dict()
    assert _diop_ternary_quadratic_normal((x, y, z), coeff) == \
        (1, 1, 2)
    assert diophantine(eq) == set([(
        2*p**2 - q**2, -2*p**2 + 4*p*q - q**2,
        4*p**2 - 4*p*q + 2*q**2)])
    eq = 411*x**2+57*y**2-221*z**2
    coeff = eq.as_coefficients_dict()
    assert _diop_ternary_quadratic_normal((x, y, z), coeff) == \
        (2021, 2645, 3066)
    assert diophantine(eq) == \
        set([(115197*p**2 - 446641*q**2, -150765*p**2 + 1355172*p*q -
        584545*q**2, 174762*p**2 - 301530*p*q + 677586*q**2)])
    eq = 573*x**2+267*y**2-984*z**2
    coeff = eq.as_coefficients_dict()
    assert _diop_ternary_quadratic_normal((x, y, z), coeff) == \
        (49, 233, 127)
    assert diophantine(eq) == \
        set([(4361*p**2 - 16072*q**2, -20737*p**2 + 83312*p*q - 76424*q**2,
        11303*p**2 - 41474*p*q + 41656*q**2)])
    # this produces factors during reconstruction
    eq = x**2 + 3*y**2 - 12*z**2
    coeff = eq.as_coefficients_dict()
    assert _diop_ternary_quadratic_normal((x, y, z), coeff) == \
        (0, 2, 1)
    assert diophantine(eq) == \
        set([(24*p*q, 2*p**2 - 24*q**2, p**2 + 12*q**2)])
    # solvers have not been written for every type
    raises(NotImplementedError, lambda: diophantine(x*y**2 + 1))

    # rational expressions
    assert diophantine(1/x) == set()
    assert diophantine(1/x + 1/y - S.Half)
    set([(6, 3), (-2, 1), (4, 4), (1, -2), (3, 6)])

def test_issue_9539():
    assert diophantine(6*w + 9*y + 20*x - z) == \
        set([(t_0, t_1, t_1 + t_2, 6*t_0 + 29*t_1 + 9*t_2)])