def test_ilcm():
    assert ilcm(0, 0) == 0
    assert ilcm(1, 0) == 0
    assert ilcm(0, 1) == 0
    assert ilcm(1, 1) == 1
    assert ilcm(2, 1) == 2
    assert ilcm(8, 2) == 8
    assert ilcm(8, 6) == 24
    assert ilcm(8, 7) == 56
    raises(ValueError, lambda: ilcm(8.1, 7))
    raises(ValueError, lambda: ilcm(8, 7.1))

def test_ilcm():
    assert ilcm(0, 0) == 0
    assert ilcm(1, 0) == 0
    assert ilcm(0, 1) == 0
    assert ilcm(1, 1) == 1
    assert ilcm(2, 1) == 2
    assert ilcm(8, 2) == 8
    assert ilcm(8, 6) == 24
    assert ilcm(8, 7) == 56

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

        symbol = symbols[0]

        strategy = guess_solve_strategy(f, symbol)

        if strategy == GS_POLY:
            poly = f.as_poly( symbol )
            assert poly is not None
            result = roots(poly, cubics=True, quartics=True).keys()

        elif strategy == GS_RATIONAL:
            P, Q = f.as_numer_denom()
            #TODO: check for Q != 0
            return solve(P, symbol, **flags)

        elif strategy == GS_POLY_CV_1:
            args = list(f.args)
            if isinstance(f, Add):
                # we must search for a suitable change of variable
                # collect exponents
                exponents_denom = list()
                for arg in args:
                    if isinstance(arg, Pow):
                        exponents_denom.append(arg.exp.q)
                    elif isinstance(arg, Mul):
                        for mul_arg in arg.args:
                            if isinstance(mul_arg, Pow):
                                exponents_denom.append(mul_arg.exp.q)
                assert len(exponents_denom) > 0
                if len(exponents_denom) == 1:
                    m = exponents_denom[0]
                else:
                    # get the GCD of the denominators
                    m = ilcm(*exponents_denom)
                # x -> y**m.
                # we assume positive for simplification purposes
                t = Symbol('t', positive=True, dummy=True)
                f_ = f.subs(symbol, t**m)
                if guess_solve_strategy(f_, t) != GS_POLY:
                    raise TypeError("Could not convert to a polynomial equation: %s" % f_)
                cv_sols = solve(f_, t)
                for sol in cv_sols:
                    result.append(sol**m)

            elif isinstance(f, Mul):
                for mul_arg in args:
                    result.extend(solve(mul_arg, symbol))

        elif strategy == GS_POLY_CV_2:
            m = 0
            args = list(f.args)
            if isinstance(f, Add):
                for arg in args:
                    if isinstance(arg, Pow):
                        m = min(m, arg.exp)
                    elif isinstance(arg, Mul):
                        for mul_arg in arg.args:
                            if isinstance(mul_arg, Pow):
                                m = min(m, mul_arg.exp)
            elif isinstance(f, Mul):
                for mul_arg in args:
                    if isinstance(mul_arg, Pow):
                        m = min(m, mul_arg.exp)
            f1 = simplify(f*symbol**(-m))
            result = solve(f1, symbol)
            # TODO: we might have introduced unwanted solutions

    def _intersect(self, other):
        from sympy.functions.elementary.integers import ceiling, floor
        from sympy.functions.elementary.complexes import sign

        if other is S.Naturals:
            return self._intersect(Interval(1, S.Infinity))

        if other is S.Integers:
            return self

        if other.is_Interval:
            if not all(i.is_number for i in other.args[:2]):
                return

            # In case of null Range, return an EmptySet.
            if self.size == 0:
                return S.EmptySet

            # trim down to self's size, and represent
            # as a Range with step 1.
            start = ceiling(max(other.inf, self.inf))
            if start not in other:
                start += 1
            end = floor(min(other.sup, self.sup))
            if end not in other:
                end -= 1
            return self.intersect(Range(start, end + 1))

        if isinstance(other, Range):
            from sympy.solvers.diophantine import diop_linear
            from sympy.core.numbers import ilcm

            # non-overlap quick exits
            if not other:
                return S.EmptySet
            if not self:
                return S.EmptySet
            if other.sup < self.inf:
                return S.EmptySet
            if other.inf > self.sup:
                return S.EmptySet

            # work with finite end at the start
            r1 = self
            if r1.start.is_infinite:
                r1 = r1.reversed
            r2 = other
            if r2.start.is_infinite:
                r2 = r2.reversed

            # this equation represents the values of the Range;
            # it's a linear equation
            eq = lambda r, i: r.start + i*r.step

            # we want to know when the two equations might
            # have integer solutions so we use the diophantine
            # solver
            a, b = diop_linear(eq(r1, Dummy()) - eq(r2, Dummy()))

            # check for no solution
            no_solution = a is None and b is None
            if no_solution:
                return S.EmptySet

            # there is a solution
            # -------------------

            # find the coincident point, c
            a0 = a.as_coeff_Add()[0]
            c = eq(r1, a0)

            # find the first point, if possible, in each range
            # since c may not be that point
            def _first_finite_point(r1, c):
                if c == r1.start:
                    return c
                # st is the signed step we need to take to
                # get from c to r1.start
                st = sign(r1.start - c)*step
                # use Range to calculate the first point:
                # we want to get as close as possible to
                # r1.start; the Range will not be null since
                # it will at least contain c
                s1 = Range(c, r1.start + st, st)[-1]
                if s1 == r1.start:
                    pass
                else:
                    # if we didn't hit r1.start then, if the
                    # sign of st didn't match the sign of r1.step
                    # we are off by one and s1 is not in r1
                    if sign(r1.step) != sign(st):
                        s1 -= st
                if s1 not in r1:
                    return
                return s1

            # calculate the step size of the new Range
            step = abs(ilcm(r1.step, r2.step))
            s1 = _first_finite_point(r1, c)
            if s1 is None:
#.........这里部分代码省略.........

def solve(f, *symbols, **flags):
    """Solves equations and systems of equations.

       Currently supported are univariate polynomial and transcendental
       equations and systems of linear and polynomial equations.  Input
       is formed as a single expression or an equation,  or an iterable
       container in case of an equation system.  The type of output may
       vary and depends heavily on the input. For more details refer to
       more problem specific functions.

       By default all solutions are simplified to make the output more
       readable. If this is not the expected behavior,  eg. because of
       speed issues, set simplified=False in function arguments.

       To solve equations and systems of equations of other kind, eg.
       recurrence relations of differential equations use rsolve() or
       dsolve() functions respectively.

       >>> from sympy import *
       >>> x,y = symbols('xy')

       Solve a polynomial equation:

       >>> solve(x**4-1, x)
       [1, -1, -I, I]

       Solve a linear system:

       >>> solve((x+5*y-2, -3*x+6*y-15), x, y)
       {x: -3, y: 1}

    """
    if not symbols:
        raise ValueError('no symbols were given')

    if len(symbols) == 1:
        if isinstance(symbols[0], (list, tuple, set)):
            symbols = symbols[0]

    symbols = map(sympify, symbols)

    if any(not s.is_Symbol for s in symbols):
        raise TypeError('not a Symbol')

    if not isinstance(f, (tuple, list, set)):
        f = sympify(f)

        if isinstance(f, Equality):
            f = f.lhs - f.rhs

        if len(symbols) != 1:
            raise NotImplementedError('multivariate equation')

        symbol = symbols[0]

        strategy = guess_solve_strategy(f, symbol)

        if strategy == GS_POLY:
            poly = f.as_poly( symbol )
            assert poly is not None
            result = roots(poly, cubics=True, quartics=True).keys()

        elif strategy == GS_RATIONAL:
            P, Q = f.as_numer_denom()
            #TODO: check for Q != 0
            return solve(P, symbol, **flags)

        elif strategy == GS_POLY_CV_1:
            args = list(f.args)
            if isinstance(f, Add):
                # we must search for a suitable change of variable
                # collect exponents
                exponents_denom = list()
                for arg in args:
                    if isinstance(arg, Pow):
                        exponents_denom.append(arg.exp.q)
                    elif isinstance(arg, Mul):
                        for mul_arg in arg.args:
                            if isinstance(mul_arg, Pow):
                                exponents_denom.append(mul_arg.exp.q)
                assert len(exponents_denom) > 0
                if len(exponents_denom) == 1:
                    m = exponents_denom[0]
                else:
                    # get the GCD of the denominators
                    m = ilcm(*exponents_denom)
                # x -> y**m.
                # we assume positive for simplification purposes
                t = Symbol('t', positive=True, dummy=True)
                f_ = f.subs(symbol, t**m)
                if guess_solve_strategy(f_, t) != GS_POLY:
                    raise TypeError("Could not convert to a polynomial equation: %s" % f_)
                cv_sols = solve(f_, t)
                result = list()
                for sol in cv_sols:
                    result.append(sol**m)

            elif isinstance(f, Mul):
                result = []
                for mul_arg in args:
#.........这里部分代码省略.........

def poly_lcm(f, g, *symbols):
    """Computes least common multiple of two polynomials.

       Given two univariate polynomials,  the LCM is computed  via
       f*g = gcd(f, g)*lcm(f, g) formula. In multivariate case, we
       compute the unique generator of the intersection of the two
       ideals, generated by f and g.  This is done by computing  a
       Groebner basis, with respect to any lexicographic ordering,
       of t*f and (1 - t)*g, where t is an unrelated symbol and
       filtering out solution that does not contain t.

       For more information on the implemented algorithm refer to:

       [1] D. Cox, J. Little, D. O'Shea, Ideals, Varieties and
           Algorithms, Springer, Second Edition, 1997, pp. 187

    """
    if not isinstance(f, Poly):
        f = Poly(f, *symbols)
    elif symbols:
        raise SymbolsError("Redundant symbols were given")

    f, g = f.unify_with(g)

    symbols, flags = f.symbols, f.flags

    if f.is_monomial and g.is_monomial:
        monom = monomial_lcm(f.LM, g.LM)

        fc, gc = f.LC, g.LC

        if fc.is_Rational and gc.is_Rational:
            coeff = Integer(ilcm(fc.p, gc.p))
        else:
            coeff = S.One

        return Poly((coeff, monom), *symbols, **flags)

    fc, f = f.as_primitive()
    gc, g = g.as_primitive()

    lcm = ilcm(int(fc), int(gc))

    if f.is_multivariate:
        t = Symbol('t', dummy=True)
        lex = { 'order' : 'lex' }

        f_monoms = [ (1,) + monom for monom in f.monoms ]

        F = Poly((f.coeffs, f_monoms), t, *symbols, **lex)

        g_monoms = [ (0,) + monom for monom in g.monoms ] + \
                   [ (1,) + monom for monom in g.monoms ]

        g_coeffs = list(g.coeffs) + [ -coeff for coeff in g.coeffs ]
        G = Poly(dict(zip(g_monoms, g_coeffs)), t, *symbols, **lex)

        def independent(h):
            return all(not monom[0] for monom in h.monoms)

        H = [ h for h in poly_groebner((F, G)) if independent(h) ]

        if lcm != 1:
            h_coeffs = [ coeff*lcm for coeff in H[0].coeffs ]
        else:
            h_coeffs = H[0].coeffs

        h_monoms = [ monom[1:] for monom in H[0].monoms ]

        return Poly(dict(zip(h_monoms, h_coeffs)), *symbols, **flags)
    else:
        h = poly_div(f * g, poly_gcd(f, g))[0]

        if lcm != 1:
            return h.mul_term(lcm / h.LC)
        else:
            return h.as_monic()

def intersection_sets(a, b):
    from sympy.solvers.diophantine import diop_linear
    from sympy.core.numbers import ilcm
    from sympy import sign

    # non-overlap quick exits
    if not b:
        return S.EmptySet
    if not a:
        return S.EmptySet
    if b.sup < a.inf:
        return S.EmptySet
    if b.inf > a.sup:
        return S.EmptySet

    # work with finite end at the start
    r1 = a
    if r1.start.is_infinite:
        r1 = r1.reversed
    r2 = b
    if r2.start.is_infinite:
        r2 = r2.reversed

    # this equation represents the values of the Range;
    # it's a linear equation
    eq = lambda r, i: r.start + i*r.step

    # we want to know when the two equations might
    # have integer solutions so we use the diophantine
    # solver
    va, vb = diop_linear(eq(r1, Dummy()) - eq(r2, Dummy()))

    # check for no solution
    no_solution = va is None and vb is None
    if no_solution:
        return S.EmptySet

    # there is a solution
    # -------------------

    # find the coincident point, c
    a0 = va.as_coeff_Add()[0]
    c = eq(r1, a0)

    # find the first point, if possible, in each range
    # since c may not be that point
    def _first_finite_point(r1, c):
        if c == r1.start:
            return c
        # st is the signed step we need to take to
        # get from c to r1.start
        st = sign(r1.start - c)*step
        # use Range to calculate the first point:
        # we want to get as close as possible to
        # r1.start; the Range will not be null since
        # it will at least contain c
        s1 = Range(c, r1.start + st, st)[-1]
        if s1 == r1.start:
            pass
        else:
            # if we didn't hit r1.start then, if the
            # sign of st didn't match the sign of r1.step
            # we are off by one and s1 is not in r1
            if sign(r1.step) != sign(st):
                s1 -= st
        if s1 not in r1:
            return
        return s1

    # calculate the step size of the new Range
    step = abs(ilcm(r1.step, r2.step))
    s1 = _first_finite_point(r1, c)
    if s1 is None:
        return S.EmptySet
    s2 = _first_finite_point(r2, c)
    if s2 is None:
        return S.EmptySet

    # replace the corresponding start or stop in
    # the original Ranges with these points; the
    # result must have at least one point since
    # we know that s1 and s2 are in the Ranges
    def _updated_range(r, first):
        st = sign(r.step)*step
        if r.start.is_finite:
            rv = Range(first, r.stop, st)
        else:
            rv = Range(r.start, first + st, st)
        return rv
    r1 = _updated_range(a, s1)
    r2 = _updated_range(b, s2)

    # work with them both in the increasing direction
    if sign(r1.step) < 0:
        r1 = r1.reversed
    if sign(r2.step) < 0:
        r2 = r2.reversed

    # return clipped Range with positive step; it
    # can't be empty at this point
#.........这里部分代码省略.........