from ..core import Add, Dummy, Pow, expand_log
from ..core.function import _mexpand
from ..functions import LambertW, exp, log, root
from ..polys.polytools import Poly, factor
from ..simplify import collect, separatevars
from ..utilities import ordered
from .solvers import _invert, solve


def _filtered_gens(poly, symbol):
    """Process the generators of ``poly``, returning the set of generators that
    have ``symbol``.  If there are two generators that are inverses of each other,
    prefer the one that has no denominator.

    Examples
    ========

    >>> _filtered_gens((x + 1/x + exp(x)).as_poly(), x)
    {E**x, x}

    """
    gens = {g for g in poly.gens if symbol in g.free_symbols}
    for g in list(gens):
        ag = 1/g
        if g in gens and ag in gens:
            if ag.as_numer_denom()[1] != 1:
                g = ag
            gens.remove(g)
    return gens


def _mostfunc(lhs, func, X=None):
    """Returns the term in lhs which contains the most of the
    func-type things e.g. log(log(x)) wins over log(x) if both terms appear.

    ``func`` can be a function (exp, log, etc...) or any other Diofant object,
    like Pow.

    If ``X`` is not ``None``, then the function returns the term composed with the
    most ``func`` having the specified variable.

    Examples
    ========

    >>> _mostfunc(exp(x) + exp(exp(x) + 2), Pow)
    E**(E**x + 2)
    >>> _mostfunc(exp(x) + exp(exp(y) + 2), Pow, x)
    E**x
    >>> _mostfunc(exp(x) + exp(exp(y) + 2), Pow)
    E**(E**y + 2)
    >>> _mostfunc(x, exp, x) is None
    True
    >>> _mostfunc(exp(x) + exp(x*y), Pow, x)
    E**x

    """
    fterms = [tmp for tmp in lhs.atoms(func) if (not X or
                                                 X.is_Symbol and X in tmp.free_symbols or
                                                 not X.is_Symbol and tmp.has(X))]
    if len(fterms) == 1:
        return fterms[0]
    if fterms:
        return max(list(ordered(fterms)), key=lambda x: x.count(func))


def _linab(arg, symbol):
    """Return ``a, b, X`` assuming ``arg`` can be written as ``a*X + b``
    where ``X`` is a symbol-dependent factor and ``a`` and ``b`` are
    independent of ``symbol``.

    Examples
    ========

    >>> _linab(Integer(2), x)
    (2, 0, 1)
    >>> _linab(2*x, x)
    (2, 0, x)
    >>> _linab(y + y*x + 2*x, x)
    (y + 2, y, x)
    >>> _linab(3 + 2*exp(x), x)
    (2, 3, E**x)

    """
    arg = arg.expand()
    ind, dep = arg.as_independent(symbol)
    if not arg.is_Add:
        b = 0
        a, x = ind, dep
    else:
        b = ind
        a, x = separatevars(dep).as_independent(symbol, as_Add=False)
    if x.could_extract_minus_sign():
        a = -a
        x = -x
    return a, b, x


def _lambert(eq, x):
    """
    Given an expression assumed to be in the form
        ``F(X, a..f) = a*log(b*X + c) + d*X + f = 0``
    where X = g(x) and x = g^-1(X), return the Lambert solution if possible:
        ``x = g^-1(-c/b + (a/d)*W(d/(a*b)*exp(c*d/a/b)*exp(-f/a)))``.

    """
    eq = _mexpand(expand_log(eq))
    mainlog = _mostfunc(eq, log, x)
    if not mainlog:
        return []  # violated assumptions
    other = eq.subs({mainlog: 0})
    if isinstance(-other, log):
        eq = (eq - other).subs({mainlog: mainlog.args[0]})
        mainlog = mainlog.args[0]
        if not isinstance(mainlog, log):
            return []  # violated assumptions
        other = -(-other).args[0]
        eq += other
    if x not in other.free_symbols:
        return []  # violated assumptions
    d, f, X2 = _linab(other, x)
    logterm = collect(eq - other, mainlog)
    a = logterm.as_coefficient(mainlog)
    if a is None or x in a.free_symbols:
        return []  # violated assumptions
    logarg = mainlog.args[0]
    b, c, X1 = _linab(logarg, x)
    if X1 != X2:
        return []  # violated assumptions

    u = Dummy('rhs')
    sol = []
    # check only real solutions:
    for k in [-1, 0]:
        l = LambertW(d/(a*b)*exp(c*d/a/b)*exp(-f/a), k)
        # if W's arg is between -1/e and 0 there is
        # a -1 branch real solution, too.
        if k and not l.is_extended_real:
            continue
        rhs = -c/b + (a/d)*l

        solns = solve(X1 - u, x)
        for i, tmp in enumerate(solns):
            solns[i] = tmp[x].subs({u: rhs})
            sol.append(solns[i])
    return sol


def _solve_lambert(f, symbol, gens):
    """Return solution to ``f`` if it is a Lambert-type expression
    else raise NotImplementedError.

    The equality, ``f(x, a..f) = a*log(b*X + c) + d*X - f = 0`` has the
    solution,  `X = -c/b + (a/d)*W(d/(a*b)*exp(c*d/a/b)*exp(f/a))`. There
    are a variety of forms for `f(X, a..f)` as enumerated below:

    1a1)
      if B**B = R for R not [0, 1] then
      log(B) + log(log(B)) = log(log(R))
      X = log(B), a = 1, b = 1, c = 0, d = 1, f = log(log(R))
    1a2)
      if B*(b*log(B) + c)**a = R then
      log(B) + a*log(b*log(B) + c) = log(R)
      X = log(B); d=1, f=log(R)
    1b)
      if a*log(b*B + c) + d*B = R then
      X = B, f = R
    2)
      if d*p**(a*B + g) - b*B = c then
      log(d) + (a*B + g)*log(p) - log(c + b*B) = 0
      a = -1, d = a*log(p), f = -log(d) - g*log(p)

    """
    nrhs, lhs = f.as_independent(symbol, as_Add=True)
    rhs = -nrhs

    lamcheck = [tmp for tmp in gens
                if (tmp.func in [exp, log] or
                    (tmp.is_Pow and symbol in tmp.exp.free_symbols))]
    if not lamcheck:
        raise NotImplementedError()

    if lhs.is_Mul:
        lhs = expand_log(log(lhs))
        rhs = log(rhs)

    lhs = factor(lhs, deep=True)
    # make sure we have inverted as completely as possible
    r = Dummy()
    i, lhs = _invert(lhs - r, symbol)
    rhs = i.xreplace({r: rhs})

    # For the first ones:
    # 1a1) B**B = R != 0 (when 0, there is only a solution if the base is 0,
    #                     but if it is, the exp is 0 and 0**0=1
    #                     comes back as B*log(B) = log(R)
    # 1a2) B*(a + b*log(B))**p = R or with monomial expanded or with whole
    #                              thing expanded comes back unchanged
    #     log(B) + p*log(a + b*log(B)) = log(R)
    #     lhs is Mul:
    #         expand log of both sides to give:
    #         log(B) + log(log(B)) = log(log(R))
    # 1b) d*log(a*B + b) + c*B = R
    #     lhs is Add:
    #         isolate c*B and expand log of both sides:
    #         log(c) + log(B) = log(R - d*log(a*B + b))

    soln = []
    mainlog = _mostfunc(lhs, log, symbol)
    if mainlog:
        if lhs.is_Mul and rhs != 0:
            soln = _lambert(log(lhs) - log(rhs), symbol)
        elif lhs.is_Add:
            other = lhs.subs({mainlog: 0})
            if other and not other.is_Add and [
                    tmp for tmp in other.atoms(Pow)
                    if symbol in tmp.free_symbols]:
                if not rhs:
                    diff = log(other) - log(other - lhs)
                else:
                    diff = log(lhs - other) - log(rhs - other)
                soln = _lambert(expand_log(diff), symbol)
            else:
                # it's ready to go
                soln = _lambert(lhs - rhs, symbol)

    # 2) d*p**(a*B + b) + c*B = R
    #     collect on main pow
    #     log(R - c*B) - a*B*log(p) = log(d) + b*log(p)

    if not soln:
        mainpow = _mostfunc(lhs, Pow, symbol)
        if mainpow and symbol in mainpow.exp.free_symbols:
            lhs = collect(lhs, mainpow)
            if lhs.is_Mul and rhs != 0:
                soln = _lambert(expand_log(log(lhs) - log(rhs)), symbol)
            elif lhs.is_Add:
                # move all but mainpow-containing term to rhs
                other = lhs.subs({mainpow: 0})
                mainterm = lhs - other
                rhs = rhs - other
                diff = log(mainterm) - log(rhs)
                soln = _lambert(expand_log(diff), symbol)

    if not soln:
        raise NotImplementedError(f'{f} does not appear to have a solution in '
                                  'terms of LambertW')

    return list(ordered(soln))


def bivariate_type(f, x, y, **kwargs):
    """Given an expression, f, 3 tests will be done to see what type
    of composite bivariate it might be, options for u(x, y) are::

        x*y
        x+y
        x*y+x
        x*y+y

    If it matches one of these types, ``u(x, y)``, ``P(u)`` and dummy
    variable ``u`` will be returned. Solving ``P(u)`` for ``u`` and
    equating the solutions to ``u(x, y)`` and then solving for ``x`` or
    ``y`` is equivalent to solving the original expression for ``x`` or
    ``y``. If ``x`` and ``y`` represent two functions in the same
    variable, e.g. ``x = g(t)`` and ``y = h(t)``, then if ``u(x, y) - p``
    can be solved for ``t`` then these represent the solutions to
    ``P(u) = 0`` when ``p`` are the solutions of ``P(u) = 0``.

    Only positive values of ``u`` are considered.

    Examples
    ========

    >>> eq = (x**2 - 3).subs({x: x + y})
    >>> bivariate_type(eq, x, y)
    (x + y, _u**2 - 3, _u)
    >>> uxy, pu, u = _
    >>> usol = solve(pu, u)
    >>> usol
    [{_u: sqrt(3)}]
    >>> [solve(uxy - s[u]) for s in solve(pu, u)]
    [[{x: -y + sqrt(3)}]]
    >>> all(eq.subs(s).equals(0) for sol in _ for s in sol)
    True

    """
    u = Dummy('u', positive=True)

    if kwargs.pop('first', True):
        p = Poly(f, x, y)
        f = p.as_expr()
        _x = Dummy()
        _y = Dummy()
        rv = bivariate_type(Poly(f.subs({x: _x, y: _y}), _x, _y), _x, _y, first=False)
        if rv:
            reps = {_x: x, _y: y}
            return rv[0].xreplace(reps), rv[1].xreplace(reps), rv[2]
        return

    p = f
    f = p.as_expr()

    # f(x*y)
    args = Add.make_args(p.as_expr())
    new = []
    for a in args:
        a = _mexpand(a.subs({x: u/y}))
        free = a.free_symbols
        if x in free or y in free:
            break
        new.append(a)
    else:
        return x*y, Add(*new), u

    def ok(f, v, c):
        new = _mexpand(f.subs({v: c}))
        free = new.free_symbols
        return None if (x in free or y in free) else new

    # f(a*x + b*y)
    new = []
    d = p.degree(x)
    if p.degree(y) == d:
        a = root(p.coeff_monomial(x**d), d)
        b = root(p.coeff_monomial(y**d), d)
        new = ok(f, x, (u - b*y)/a)
        if new is not None:
            return a*x + b*y, new, u

    # f(a*x*y + b*y)
    new = []
    d = p.degree(x)
    if p.degree(y) == d:
        for _ in range(2):
            a = root(p.coeff_monomial(x**d*y**d), d)
            b = root(p.coeff_monomial(y**d), d)
            new = ok(f, x, (u - b*y)/a/y)
            if new is not None:
                return a*x*y + b*y, new, u
            x, y = y, x

8c8a619cc5eefb811465a150250894b66cd61e91
<mask token>


def _linab(arg, symbol):
    """Return ``a, b, X`` assuming ``arg`` can be written as ``a*X + b``
    where ``X`` is a symbol-dependent factor and ``a`` and ``b`` are
    independent of ``symbol``.

    Examples
    ========

    >>> _linab(Integer(2), x)
    (2, 0, 1)
    >>> _linab(2*x, x)
    (2, 0, x)
    >>> _linab(y + y*x + 2*x, x)
    (y + 2, y, x)
    >>> _linab(3 + 2*exp(x), x)
    (2, 3, E**x)

    """
    arg = arg.expand()
    ind, dep = arg.as_independent(symbol)
    if not arg.is_Add:
        b = 0
        a, x = ind, dep
    else:
        b = ind
        a, x = separatevars(dep).as_independent(symbol, as_Add=False)
    if x.could_extract_minus_sign():
        a = -a
        x = -x
    return a, b, x


<mask token>


def _solve_lambert(f, symbol, gens):
    """Return solution to ``f`` if it is a Lambert-type expression
    else raise NotImplementedError.

    The equality, ``f(x, a..f) = a*log(b*X + c) + d*X - f = 0`` has the
    solution,  `X = -c/b + (a/d)*W(d/(a*b)*exp(c*d/a/b)*exp(f/a))`. There
    are a variety of forms for `f(X, a..f)` as enumerated below:

    1a1)
      if B**B = R for R not [0, 1] then
      log(B) + log(log(B)) = log(log(R))
      X = log(B), a = 1, b = 1, c = 0, d = 1, f = log(log(R))
    1a2)
      if B*(b*log(B) + c)**a = R then
      log(B) + a*log(b*log(B) + c) = log(R)
      X = log(B); d=1, f=log(R)
    1b)
      if a*log(b*B + c) + d*B = R then
      X = B, f = R
    2)
      if d*p**(a*B + g) - b*B = c then
      log(d) + (a*B + g)*log(p) - log(c + b*B) = 0
      a = -1, d = a*log(p), f = -log(d) - g*log(p)

    """
    nrhs, lhs = f.as_independent(symbol, as_Add=True)
    rhs = -nrhs
    lamcheck = [tmp for tmp in gens if tmp.func in [exp, log] or tmp.is_Pow and
        symbol in tmp.exp.free_symbols]
    if not lamcheck:
        raise NotImplementedError()
    if lhs.is_Mul:
        lhs = expand_log(log(lhs))
        rhs = log(rhs)
    lhs = factor(lhs, deep=True)
    r = Dummy()
    i, lhs = _invert(lhs - r, symbol)
    rhs = i.xreplace({r: rhs})
    soln = []
    mainlog = _mostfunc(lhs, log, symbol)
    if mainlog:
        if lhs.is_Mul and rhs != 0:
            soln = _lambert(log(lhs) - log(rhs), symbol)
        elif lhs.is_Add:
            other = lhs.subs({mainlog: 0})
            if other and not other.is_Add and [tmp for tmp in other.atoms(
                Pow) if symbol in tmp.free_symbols]:
                if not rhs:
                    diff = log(other) - log(other - lhs)
                else:
                    diff = log(lhs - other) - log(rhs - other)
                soln = _lambert(expand_log(diff), symbol)
            else:
                soln = _lambert(lhs - rhs, symbol)
    if not soln:
        mainpow = _mostfunc(lhs, Pow, symbol)
        if mainpow and symbol in mainpow.exp.free_symbols:
            lhs = collect(lhs, mainpow)
            if lhs.is_Mul and rhs != 0:
                soln = _lambert(expand_log(log(lhs) - log(rhs)), symbol)
            elif lhs.is_Add:
                other = lhs.subs({mainpow: 0})
                mainterm = lhs - other
                rhs = rhs - other
                diff = log(mainterm) - log(rhs)
                soln = _lambert(expand_log(diff), symbol)
    if not soln:
        raise NotImplementedError(
            f'{f} does not appear to have a solution in terms of LambertW')
    return list(ordered(soln))


<mask token>

<mask token>


def _mostfunc(lhs, func, X=None):
    """Returns the term in lhs which contains the most of the
    func-type things e.g. log(log(x)) wins over log(x) if both terms appear.

    ``func`` can be a function (exp, log, etc...) or any other Diofant object,
    like Pow.

    If ``X`` is not ``None``, then the function returns the term composed with the
    most ``func`` having the specified variable.

    Examples
    ========

    >>> _mostfunc(exp(x) + exp(exp(x) + 2), Pow)
    E**(E**x + 2)
    >>> _mostfunc(exp(x) + exp(exp(y) + 2), Pow, x)
    E**x
    >>> _mostfunc(exp(x) + exp(exp(y) + 2), Pow)
    E**(E**y + 2)
    >>> _mostfunc(x, exp, x) is None
    True
    >>> _mostfunc(exp(x) + exp(x*y), Pow, x)
    E**x

    """
    fterms = [tmp for tmp in lhs.atoms(func) if not X or X.is_Symbol and X in
        tmp.free_symbols or not X.is_Symbol and tmp.has(X)]
    if len(fterms) == 1:
        return fterms[0]
    if fterms:
        return max(list(ordered(fterms)), key=lambda x: x.count(func))


def _linab(arg, symbol):
    """Return ``a, b, X`` assuming ``arg`` can be written as ``a*X + b``
    where ``X`` is a symbol-dependent factor and ``a`` and ``b`` are
    independent of ``symbol``.

    Examples
    ========

    >>> _linab(Integer(2), x)
    (2, 0, 1)
    >>> _linab(2*x, x)
    (2, 0, x)
    >>> _linab(y + y*x + 2*x, x)
    (y + 2, y, x)
    >>> _linab(3 + 2*exp(x), x)
    (2, 3, E**x)

    """
    arg = arg.expand()
    ind, dep = arg.as_independent(symbol)
    if not arg.is_Add:
        b = 0
        a, x = ind, dep
    else:
        b = ind
        a, x = separatevars(dep).as_independent(symbol, as_Add=False)
    if x.could_extract_minus_sign():
        a = -a
        x = -x
    return a, b, x


def _lambert(eq, x):
    """
    Given an expression assumed to be in the form
        ``F(X, a..f) = a*log(b*X + c) + d*X + f = 0``
    where X = g(x) and x = g^-1(X), return the Lambert solution if possible:
        ``x = g^-1(-c/b + (a/d)*W(d/(a*b)*exp(c*d/a/b)*exp(-f/a)))``.

    """
    eq = _mexpand(expand_log(eq))
    mainlog = _mostfunc(eq, log, x)
    if not mainlog:
        return []
    other = eq.subs({mainlog: 0})
    if isinstance(-other, log):
        eq = (eq - other).subs({mainlog: mainlog.args[0]})
        mainlog = mainlog.args[0]
        if not isinstance(mainlog, log):
            return []
        other = -(-other).args[0]
        eq += other
    if x not in other.free_symbols:
        return []
    d, f, X2 = _linab(other, x)
    logterm = collect(eq - other, mainlog)
    a = logterm.as_coefficient(mainlog)
    if a is None or x in a.free_symbols:
        return []
    logarg = mainlog.args[0]
    b, c, X1 = _linab(logarg, x)
    if X1 != X2:
        return []
    u = Dummy('rhs')
    sol = []
    for k in [-1, 0]:
        l = LambertW(d / (a * b) * exp(c * d / a / b) * exp(-f / a), k)
        if k and not l.is_extended_real:
            continue
        rhs = -c / b + a / d * l
        solns = solve(X1 - u, x)
        for i, tmp in enumerate(solns):
            solns[i] = tmp[x].subs({u: rhs})
            sol.append(solns[i])
    return sol


def _solve_lambert(f, symbol, gens):
    """Return solution to ``f`` if it is a Lambert-type expression
    else raise NotImplementedError.

    The equality, ``f(x, a..f) = a*log(b*X + c) + d*X - f = 0`` has the
    solution,  `X = -c/b + (a/d)*W(d/(a*b)*exp(c*d/a/b)*exp(f/a))`. There
    are a variety of forms for `f(X, a..f)` as enumerated below:

    1a1)
      if B**B = R for R not [0, 1] then
      log(B) + log(log(B)) = log(log(R))
      X = log(B), a = 1, b = 1, c = 0, d = 1, f = log(log(R))
    1a2)
      if B*(b*log(B) + c)**a = R then
      log(B) + a*log(b*log(B) + c) = log(R)
      X = log(B); d=1, f=log(R)
    1b)
      if a*log(b*B + c) + d*B = R then
      X = B, f = R
    2)
      if d*p**(a*B + g) - b*B = c then
      log(d) + (a*B + g)*log(p) - log(c + b*B) = 0
      a = -1, d = a*log(p), f = -log(d) - g*log(p)

    """
    nrhs, lhs = f.as_independent(symbol, as_Add=True)
    rhs = -nrhs
    lamcheck = [tmp for tmp in gens if tmp.func in [exp, log] or tmp.is_Pow and
        symbol in tmp.exp.free_symbols]
    if not lamcheck:
        raise NotImplementedError()
    if lhs.is_Mul:
        lhs = expand_log(log(lhs))
        rhs = log(rhs)
    lhs = factor(lhs, deep=True)
    r = Dummy()
    i, lhs = _invert(lhs - r, symbol)
    rhs = i.xreplace({r: rhs})
    soln = []
    mainlog = _mostfunc(lhs, log, symbol)
    if mainlog:
        if lhs.is_Mul and rhs != 0:
            soln = _lambert(log(lhs) - log(rhs), symbol)
        elif lhs.is_Add:
            other = lhs.subs({mainlog: 0})
            if other and not other.is_Add and [tmp for tmp in other.atoms(
                Pow) if symbol in tmp.free_symbols]:
                if not rhs:
                    diff = log(other) - log(other - lhs)
                else:
                    diff = log(lhs - other) - log(rhs - other)
                soln = _lambert(expand_log(diff), symbol)
            else:
                soln = _lambert(lhs - rhs, symbol)
    if not soln:
        mainpow = _mostfunc(lhs, Pow, symbol)
        if mainpow and symbol in mainpow.exp.free_symbols:
            lhs = collect(lhs, mainpow)
            if lhs.is_Mul and rhs != 0:
                soln = _lambert(expand_log(log(lhs) - log(rhs)), symbol)
            elif lhs.is_Add:
                other = lhs.subs({mainpow: 0})
                mainterm = lhs - other
                rhs = rhs - other
                diff = log(mainterm) - log(rhs)
                soln = _lambert(expand_log(diff), symbol)
    if not soln:
        raise NotImplementedError(
            f'{f} does not appear to have a solution in terms of LambertW')
    return list(ordered(soln))


def bivariate_type(f, x, y, **kwargs):
    """Given an expression, f, 3 tests will be done to see what type
    of composite bivariate it might be, options for u(x, y) are::

        x*y
        x+y
        x*y+x
        x*y+y

    If it matches one of these types, ``u(x, y)``, ``P(u)`` and dummy
    variable ``u`` will be returned. Solving ``P(u)`` for ``u`` and
    equating the solutions to ``u(x, y)`` and then solving for ``x`` or
    ``y`` is equivalent to solving the original expression for ``x`` or
    ``y``. If ``x`` and ``y`` represent two functions in the same
    variable, e.g. ``x = g(t)`` and ``y = h(t)``, then if ``u(x, y) - p``
    can be solved for ``t`` then these represent the solutions to
    ``P(u) = 0`` when ``p`` are the solutions of ``P(u) = 0``.

    Only positive values of ``u`` are considered.

    Examples
    ========

    >>> eq = (x**2 - 3).subs({x: x + y})
    >>> bivariate_type(eq, x, y)
    (x + y, _u**2 - 3, _u)
    >>> uxy, pu, u = _
    >>> usol = solve(pu, u)
    >>> usol
    [{_u: sqrt(3)}]
    >>> [solve(uxy - s[u]) for s in solve(pu, u)]
    [[{x: -y + sqrt(3)}]]
    >>> all(eq.subs(s).equals(0) for sol in _ for s in sol)
    True

    """
    u = Dummy('u', positive=True)
    if kwargs.pop('first', True):
        p = Poly(f, x, y)
        f = p.as_expr()
        _x = Dummy()
        _y = Dummy()
        rv = bivariate_type(Poly(f.subs({x: _x, y: _y}), _x, _y), _x, _y,
            first=False)
        if rv:
            reps = {_x: x, _y: y}
            return rv[0].xreplace(reps), rv[1].xreplace(reps), rv[2]
        return
    p = f
    f = p.as_expr()
    args = Add.make_args(p.as_expr())
    new = []
    for a in args:
        a = _mexpand(a.subs({x: u / y}))
        free = a.free_symbols
        if x in free or y in free:
            break
        new.append(a)
    else:
        return x * y, Add(*new), u

    def ok(f, v, c):
        new = _mexpand(f.subs({v: c}))
        free = new.free_symbols
        return None if x in free or y in free else new
    new = []
    d = p.degree(x)
    if p.degree(y) == d:
        a = root(p.coeff_monomial(x ** d), d)
        b = root(p.coeff_monomial(y ** d), d)
        new = ok(f, x, (u - b * y) / a)
        if new is not None:
            return a * x + b * y, new, u
    new = []
    d = p.degree(x)
    if p.degree(y) == d:
        for _ in range(2):
            a = root(p.coeff_monomial(x ** d * y ** d), d)
            b = root(p.coeff_monomial(y ** d), d)
            new = ok(f, x, (u - b * y) / a / y)
            if new is not None:
                return a * x * y + b * y, new, u
            x, y = y, x

<mask token>


def _filtered_gens(poly, symbol):
    """Process the generators of ``poly``, returning the set of generators that
    have ``symbol``.  If there are two generators that are inverses of each other,
    prefer the one that has no denominator.

    Examples
    ========

    >>> _filtered_gens((x + 1/x + exp(x)).as_poly(), x)
    {E**x, x}

    """
    gens = {g for g in poly.gens if symbol in g.free_symbols}
    for g in list(gens):
        ag = 1 / g
        if g in gens and ag in gens:
            if ag.as_numer_denom()[1] != 1:
                g = ag
            gens.remove(g)
    return gens


def _mostfunc(lhs, func, X=None):
    """Returns the term in lhs which contains the most of the
    func-type things e.g. log(log(x)) wins over log(x) if both terms appear.

    ``func`` can be a function (exp, log, etc...) or any other Diofant object,
    like Pow.

    If ``X`` is not ``None``, then the function returns the term composed with the
    most ``func`` having the specified variable.

    Examples
    ========

    >>> _mostfunc(exp(x) + exp(exp(x) + 2), Pow)
    E**(E**x + 2)
    >>> _mostfunc(exp(x) + exp(exp(y) + 2), Pow, x)
    E**x
    >>> _mostfunc(exp(x) + exp(exp(y) + 2), Pow)
    E**(E**y + 2)
    >>> _mostfunc(x, exp, x) is None
    True
    >>> _mostfunc(exp(x) + exp(x*y), Pow, x)
    E**x

    """
    fterms = [tmp for tmp in lhs.atoms(func) if not X or X.is_Symbol and X in
        tmp.free_symbols or not X.is_Symbol and tmp.has(X)]
    if len(fterms) == 1:
        return fterms[0]
    if fterms:
        return max(list(ordered(fterms)), key=lambda x: x.count(func))


def _linab(arg, symbol):
    """Return ``a, b, X`` assuming ``arg`` can be written as ``a*X + b``
    where ``X`` is a symbol-dependent factor and ``a`` and ``b`` are
    independent of ``symbol``.

    Examples
    ========

    >>> _linab(Integer(2), x)
    (2, 0, 1)
    >>> _linab(2*x, x)
    (2, 0, x)
    >>> _linab(y + y*x + 2*x, x)
    (y + 2, y, x)
    >>> _linab(3 + 2*exp(x), x)
    (2, 3, E**x)

    """
    arg = arg.expand()
    ind, dep = arg.as_independent(symbol)
    if not arg.is_Add:
        b = 0
        a, x = ind, dep
    else:
        b = ind
        a, x = separatevars(dep).as_independent(symbol, as_Add=False)
    if x.could_extract_minus_sign():
        a = -a
        x = -x
    return a, b, x


def _lambert(eq, x):
    """
    Given an expression assumed to be in the form
        ``F(X, a..f) = a*log(b*X + c) + d*X + f = 0``
    where X = g(x) and x = g^-1(X), return the Lambert solution if possible:
        ``x = g^-1(-c/b + (a/d)*W(d/(a*b)*exp(c*d/a/b)*exp(-f/a)))``.

    """
    eq = _mexpand(expand_log(eq))
    mainlog = _mostfunc(eq, log, x)
    if not mainlog:
        return []
    other = eq.subs({mainlog: 0})
    if isinstance(-other, log):
        eq = (eq - other).subs({mainlog: mainlog.args[0]})
        mainlog = mainlog.args[0]
        if not isinstance(mainlog, log):
            return []
        other = -(-other).args[0]
        eq += other
    if x not in other.free_symbols:
        return []
    d, f, X2 = _linab(other, x)
    logterm = collect(eq - other, mainlog)
    a = logterm.as_coefficient(mainlog)
    if a is None or x in a.free_symbols:
        return []
    logarg = mainlog.args[0]
    b, c, X1 = _linab(logarg, x)
    if X1 != X2:
        return []
    u = Dummy('rhs')
    sol = []
    for k in [-1, 0]:
        l = LambertW(d / (a * b) * exp(c * d / a / b) * exp(-f / a), k)
        if k and not l.is_extended_real:
            continue
        rhs = -c / b + a / d * l
        solns = solve(X1 - u, x)
        for i, tmp in enumerate(solns):
            solns[i] = tmp[x].subs({u: rhs})
            sol.append(solns[i])
    return sol


def _solve_lambert(f, symbol, gens):
    """Return solution to ``f`` if it is a Lambert-type expression
    else raise NotImplementedError.

    The equality, ``f(x, a..f) = a*log(b*X + c) + d*X - f = 0`` has the
    solution,  `X = -c/b + (a/d)*W(d/(a*b)*exp(c*d/a/b)*exp(f/a))`. There
    are a variety of forms for `f(X, a..f)` as enumerated below:

    1a1)
      if B**B = R for R not [0, 1] then
      log(B) + log(log(B)) = log(log(R))
      X = log(B), a = 1, b = 1, c = 0, d = 1, f = log(log(R))
    1a2)
      if B*(b*log(B) + c)**a = R then
      log(B) + a*log(b*log(B) + c) = log(R)
      X = log(B); d=1, f=log(R)
    1b)
      if a*log(b*B + c) + d*B = R then
      X = B, f = R
    2)
      if d*p**(a*B + g) - b*B = c then
      log(d) + (a*B + g)*log(p) - log(c + b*B) = 0
      a = -1, d = a*log(p), f = -log(d) - g*log(p)

    """
    nrhs, lhs = f.as_independent(symbol, as_Add=True)
    rhs = -nrhs
    lamcheck = [tmp for tmp in gens if tmp.func in [exp, log] or tmp.is_Pow and
        symbol in tmp.exp.free_symbols]
    if not lamcheck:
        raise NotImplementedError()
    if lhs.is_Mul:
        lhs = expand_log(log(lhs))
        rhs = log(rhs)
    lhs = factor(lhs, deep=True)
    r = Dummy()
    i, lhs = _invert(lhs - r, symbol)
    rhs = i.xreplace({r: rhs})
    soln = []
    mainlog = _mostfunc(lhs, log, symbol)
    if mainlog:
        if lhs.is_Mul and rhs != 0:
            soln = _lambert(log(lhs) - log(rhs), symbol)
        elif lhs.is_Add:
            other = lhs.subs({mainlog: 0})
            if other and not other.is_Add and [tmp for tmp in other.atoms(
                Pow) if symbol in tmp.free_symbols]:
                if not rhs:
                    diff = log(other) - log(other - lhs)
                else:
                    diff = log(lhs - other) - log(rhs - other)
                soln = _lambert(expand_log(diff), symbol)
            else:
                soln = _lambert(lhs - rhs, symbol)
    if not soln:
        mainpow = _mostfunc(lhs, Pow, symbol)
        if mainpow and symbol in mainpow.exp.free_symbols:
            lhs = collect(lhs, mainpow)
            if lhs.is_Mul and rhs != 0:
                soln = _lambert(expand_log(log(lhs) - log(rhs)), symbol)
            elif lhs.is_Add:
                other = lhs.subs({mainpow: 0})
                mainterm = lhs - other
                rhs = rhs - other
                diff = log(mainterm) - log(rhs)
                soln = _lambert(expand_log(diff), symbol)
    if not soln:
        raise NotImplementedError(
            f'{f} does not appear to have a solution in terms of LambertW')
    return list(ordered(soln))


def bivariate_type(f, x, y, **kwargs):
    """Given an expression, f, 3 tests will be done to see what type
    of composite bivariate it might be, options for u(x, y) are::

        x*y
        x+y
        x*y+x
        x*y+y

    If it matches one of these types, ``u(x, y)``, ``P(u)`` and dummy
    variable ``u`` will be returned. Solving ``P(u)`` for ``u`` and
    equating the solutions to ``u(x, y)`` and then solving for ``x`` or
    ``y`` is equivalent to solving the original expression for ``x`` or
    ``y``. If ``x`` and ``y`` represent two functions in the same
    variable, e.g. ``x = g(t)`` and ``y = h(t)``, then if ``u(x, y) - p``
    can be solved for ``t`` then these represent the solutions to
    ``P(u) = 0`` when ``p`` are the solutions of ``P(u) = 0``.

    Only positive values of ``u`` are considered.

    Examples
    ========

    >>> eq = (x**2 - 3).subs({x: x + y})
    >>> bivariate_type(eq, x, y)
    (x + y, _u**2 - 3, _u)
    >>> uxy, pu, u = _
    >>> usol = solve(pu, u)
    >>> usol
    [{_u: sqrt(3)}]
    >>> [solve(uxy - s[u]) for s in solve(pu, u)]
    [[{x: -y + sqrt(3)}]]
    >>> all(eq.subs(s).equals(0) for sol in _ for s in sol)
    True

    """
    u = Dummy('u', positive=True)
    if kwargs.pop('first', True):
        p = Poly(f, x, y)
        f = p.as_expr()
        _x = Dummy()
        _y = Dummy()
        rv = bivariate_type(Poly(f.subs({x: _x, y: _y}), _x, _y), _x, _y,
            first=False)
        if rv:
            reps = {_x: x, _y: y}
            return rv[0].xreplace(reps), rv[1].xreplace(reps), rv[2]
        return
    p = f
    f = p.as_expr()
    args = Add.make_args(p.as_expr())
    new = []
    for a in args:
        a = _mexpand(a.subs({x: u / y}))
        free = a.free_symbols
        if x in free or y in free:
            break
        new.append(a)
    else:
        return x * y, Add(*new), u

    def ok(f, v, c):
        new = _mexpand(f.subs({v: c}))
        free = new.free_symbols
        return None if x in free or y in free else new
    new = []
    d = p.degree(x)
    if p.degree(y) == d:
        a = root(p.coeff_monomial(x ** d), d)
        b = root(p.coeff_monomial(y ** d), d)
        new = ok(f, x, (u - b * y) / a)
        if new is not None:
            return a * x + b * y, new, u
    new = []
    d = p.degree(x)
    if p.degree(y) == d:
        for _ in range(2):
            a = root(p.coeff_monomial(x ** d * y ** d), d)
            b = root(p.coeff_monomial(y ** d), d)
            new = ok(f, x, (u - b * y) / a / y)
            if new is not None:
                return a * x * y + b * y, new, u
            x, y = y, x

from ..core import Add, Dummy, Pow, expand_log
from ..core.function import _mexpand
from ..functions import LambertW, exp, log, root
from ..polys.polytools import Poly, factor
from ..simplify import collect, separatevars
from ..utilities import ordered
from .solvers import _invert, solve


def _filtered_gens(poly, symbol):
    """Process the generators of ``poly``, returning the set of generators that
    have ``symbol``.  If there are two generators that are inverses of each other,
    prefer the one that has no denominator.

    Examples
    ========

    >>> _filtered_gens((x + 1/x + exp(x)).as_poly(), x)
    {E**x, x}

    """
    gens = {g for g in poly.gens if symbol in g.free_symbols}
    for g in list(gens):
        ag = 1 / g
        if g in gens and ag in gens:
            if ag.as_numer_denom()[1] != 1:
                g = ag
            gens.remove(g)
    return gens


def _mostfunc(lhs, func, X=None):
    """Returns the term in lhs which contains the most of the
    func-type things e.g. log(log(x)) wins over log(x) if both terms appear.

    ``func`` can be a function (exp, log, etc...) or any other Diofant object,
    like Pow.

    If ``X`` is not ``None``, then the function returns the term composed with the
    most ``func`` having the specified variable.

    Examples
    ========

    >>> _mostfunc(exp(x) + exp(exp(x) + 2), Pow)
    E**(E**x + 2)
    >>> _mostfunc(exp(x) + exp(exp(y) + 2), Pow, x)
    E**x
    >>> _mostfunc(exp(x) + exp(exp(y) + 2), Pow)
    E**(E**y + 2)
    >>> _mostfunc(x, exp, x) is None
    True
    >>> _mostfunc(exp(x) + exp(x*y), Pow, x)
    E**x

    """
    fterms = [tmp for tmp in lhs.atoms(func) if not X or X.is_Symbol and X in
        tmp.free_symbols or not X.is_Symbol and tmp.has(X)]
    if len(fterms) == 1:
        return fterms[0]
    if fterms:
        return max(list(ordered(fterms)), key=lambda x: x.count(func))


def _linab(arg, symbol):
    """Return ``a, b, X`` assuming ``arg`` can be written as ``a*X + b``
    where ``X`` is a symbol-dependent factor and ``a`` and ``b`` are
    independent of ``symbol``.

    Examples
    ========

    >>> _linab(Integer(2), x)
    (2, 0, 1)
    >>> _linab(2*x, x)
    (2, 0, x)
    >>> _linab(y + y*x + 2*x, x)
    (y + 2, y, x)
    >>> _linab(3 + 2*exp(x), x)
    (2, 3, E**x)

    """
    arg = arg.expand()
    ind, dep = arg.as_independent(symbol)
    if not arg.is_Add:
        b = 0
        a, x = ind, dep
    else:
        b = ind
        a, x = separatevars(dep).as_independent(symbol, as_Add=False)
    if x.could_extract_minus_sign():
        a = -a
        x = -x
    return a, b, x


def _lambert(eq, x):
    """
    Given an expression assumed to be in the form
        ``F(X, a..f) = a*log(b*X + c) + d*X + f = 0``
    where X = g(x) and x = g^-1(X), return the Lambert solution if possible:
        ``x = g^-1(-c/b + (a/d)*W(d/(a*b)*exp(c*d/a/b)*exp(-f/a)))``.

    """
    eq = _mexpand(expand_log(eq))
    mainlog = _mostfunc(eq, log, x)
    if not mainlog:
        return []
    other = eq.subs({mainlog: 0})
    if isinstance(-other, log):
        eq = (eq - other).subs({mainlog: mainlog.args[0]})
        mainlog = mainlog.args[0]
        if not isinstance(mainlog, log):
            return []
        other = -(-other).args[0]
        eq += other
    if x not in other.free_symbols:
        return []
    d, f, X2 = _linab(other, x)
    logterm = collect(eq - other, mainlog)
    a = logterm.as_coefficient(mainlog)
    if a is None or x in a.free_symbols:
        return []
    logarg = mainlog.args[0]
    b, c, X1 = _linab(logarg, x)
    if X1 != X2:
        return []
    u = Dummy('rhs')
    sol = []
    for k in [-1, 0]:
        l = LambertW(d / (a * b) * exp(c * d / a / b) * exp(-f / a), k)
        if k and not l.is_extended_real:
            continue
        rhs = -c / b + a / d * l
        solns = solve(X1 - u, x)
        for i, tmp in enumerate(solns):
            solns[i] = tmp[x].subs({u: rhs})
            sol.append(solns[i])
    return sol


def _solve_lambert(f, symbol, gens):
    """Return solution to ``f`` if it is a Lambert-type expression
    else raise NotImplementedError.

    The equality, ``f(x, a..f) = a*log(b*X + c) + d*X - f = 0`` has the
    solution,  `X = -c/b + (a/d)*W(d/(a*b)*exp(c*d/a/b)*exp(f/a))`. There
    are a variety of forms for `f(X, a..f)` as enumerated below:

    1a1)
      if B**B = R for R not [0, 1] then
      log(B) + log(log(B)) = log(log(R))
      X = log(B), a = 1, b = 1, c = 0, d = 1, f = log(log(R))
    1a2)
      if B*(b*log(B) + c)**a = R then
      log(B) + a*log(b*log(B) + c) = log(R)
      X = log(B); d=1, f=log(R)
    1b)
      if a*log(b*B + c) + d*B = R then
      X = B, f = R
    2)
      if d*p**(a*B + g) - b*B = c then
      log(d) + (a*B + g)*log(p) - log(c + b*B) = 0
      a = -1, d = a*log(p), f = -log(d) - g*log(p)

    """
    nrhs, lhs = f.as_independent(symbol, as_Add=True)
    rhs = -nrhs
    lamcheck = [tmp for tmp in gens if tmp.func in [exp, log] or tmp.is_Pow and
        symbol in tmp.exp.free_symbols]
    if not lamcheck:
        raise NotImplementedError()
    if lhs.is_Mul:
        lhs = expand_log(log(lhs))
        rhs = log(rhs)
    lhs = factor(lhs, deep=True)
    r = Dummy()
    i, lhs = _invert(lhs - r, symbol)
    rhs = i.xreplace({r: rhs})
    soln = []
    mainlog = _mostfunc(lhs, log, symbol)
    if mainlog:
        if lhs.is_Mul and rhs != 0:
            soln = _lambert(log(lhs) - log(rhs), symbol)
        elif lhs.is_Add:
            other = lhs.subs({mainlog: 0})
            if other and not other.is_Add and [tmp for tmp in other.atoms(
                Pow) if symbol in tmp.free_symbols]:
                if not rhs:
                    diff = log(other) - log(other - lhs)
                else:
                    diff = log(lhs - other) - log(rhs - other)
                soln = _lambert(expand_log(diff), symbol)
            else:
                soln = _lambert(lhs - rhs, symbol)
    if not soln:
        mainpow = _mostfunc(lhs, Pow, symbol)
        if mainpow and symbol in mainpow.exp.free_symbols:
            lhs = collect(lhs, mainpow)
            if lhs.is_Mul and rhs != 0:
                soln = _lambert(expand_log(log(lhs) - log(rhs)), symbol)
            elif lhs.is_Add:
                other = lhs.subs({mainpow: 0})
                mainterm = lhs - other
                rhs = rhs - other
                diff = log(mainterm) - log(rhs)
                soln = _lambert(expand_log(diff), symbol)
    if not soln:
        raise NotImplementedError(
            f'{f} does not appear to have a solution in terms of LambertW')
    return list(ordered(soln))


def bivariate_type(f, x, y, **kwargs):
    """Given an expression, f, 3 tests will be done to see what type
    of composite bivariate it might be, options for u(x, y) are::

        x*y
        x+y
        x*y+x
        x*y+y

    If it matches one of these types, ``u(x, y)``, ``P(u)`` and dummy
    variable ``u`` will be returned. Solving ``P(u)`` for ``u`` and
    equating the solutions to ``u(x, y)`` and then solving for ``x`` or
    ``y`` is equivalent to solving the original expression for ``x`` or
    ``y``. If ``x`` and ``y`` represent two functions in the same
    variable, e.g. ``x = g(t)`` and ``y = h(t)``, then if ``u(x, y) - p``
    can be solved for ``t`` then these represent the solutions to
    ``P(u) = 0`` when ``p`` are the solutions of ``P(u) = 0``.

    Only positive values of ``u`` are considered.

    Examples
    ========

    >>> eq = (x**2 - 3).subs({x: x + y})
    >>> bivariate_type(eq, x, y)
    (x + y, _u**2 - 3, _u)
    >>> uxy, pu, u = _
    >>> usol = solve(pu, u)
    >>> usol
    [{_u: sqrt(3)}]
    >>> [solve(uxy - s[u]) for s in solve(pu, u)]
    [[{x: -y + sqrt(3)}]]
    >>> all(eq.subs(s).equals(0) for sol in _ for s in sol)
    True

    """
    u = Dummy('u', positive=True)
    if kwargs.pop('first', True):
        p = Poly(f, x, y)
        f = p.as_expr()
        _x = Dummy()
        _y = Dummy()
        rv = bivariate_type(Poly(f.subs({x: _x, y: _y}), _x, _y), _x, _y,
            first=False)
        if rv:
            reps = {_x: x, _y: y}
            return rv[0].xreplace(reps), rv[1].xreplace(reps), rv[2]
        return
    p = f
    f = p.as_expr()
    args = Add.make_args(p.as_expr())
    new = []
    for a in args:
        a = _mexpand(a.subs({x: u / y}))
        free = a.free_symbols
        if x in free or y in free:
            break
        new.append(a)
    else:
        return x * y, Add(*new), u

    def ok(f, v, c):
        new = _mexpand(f.subs({v: c}))
        free = new.free_symbols
        return None if x in free or y in free else new
    new = []
    d = p.degree(x)
    if p.degree(y) == d:
        a = root(p.coeff_monomial(x ** d), d)
        b = root(p.coeff_monomial(y ** d), d)
        new = ok(f, x, (u - b * y) / a)
        if new is not None:
            return a * x + b * y, new, u
    new = []
    d = p.degree(x)
    if p.degree(y) == d:
        for _ in range(2):
            a = root(p.coeff_monomial(x ** d * y ** d), d)
            b = root(p.coeff_monomial(y ** d), d)
            new = ok(f, x, (u - b * y) / a / y)
            if new is not None:
                return a * x * y + b * y, new, u
            x, y = y, x

from ..core import Add, Dummy, Pow, expand_log
from ..core.function import _mexpand
from ..functions import LambertW, exp, log, root
from ..polys.polytools import Poly, factor
from ..simplify import collect, separatevars
from ..utilities import ordered
from .solvers import _invert, solve


def _filtered_gens(poly, symbol):
    """Process the generators of ``poly``, returning the set of generators that
    have ``symbol``.  If there are two generators that are inverses of each other,
    prefer the one that has no denominator.

    Examples
    ========

    >>> _filtered_gens((x + 1/x + exp(x)).as_poly(), x)
    {E**x, x}

    """
    gens = {g for g in poly.gens if symbol in g.free_symbols}
    for g in list(gens):
        ag = 1/g
        if g in gens and ag in gens:
            if ag.as_numer_denom()[1] != 1:
                g = ag
            gens.remove(g)
    return gens


def _mostfunc(lhs, func, X=None):
    """Returns the term in lhs which contains the most of the
    func-type things e.g. log(log(x)) wins over log(x) if both terms appear.

    ``func`` can be a function (exp, log, etc...) or any other Diofant object,
    like Pow.

    If ``X`` is not ``None``, then the function returns the term composed with the
    most ``func`` having the specified variable.

    Examples
    ========

    >>> _mostfunc(exp(x) + exp(exp(x) + 2), Pow)
    E**(E**x + 2)
    >>> _mostfunc(exp(x) + exp(exp(y) + 2), Pow, x)
    E**x
    >>> _mostfunc(exp(x) + exp(exp(y) + 2), Pow)
    E**(E**y + 2)
    >>> _mostfunc(x, exp, x) is None
    True
    >>> _mostfunc(exp(x) + exp(x*y), Pow, x)
    E**x

    """
    fterms = [tmp for tmp in lhs.atoms(func) if (not X or
                                                 X.is_Symbol and X in tmp.free_symbols or
                                                 not X.is_Symbol and tmp.has(X))]
    if len(fterms) == 1:
        return fterms[0]
    if fterms:
        return max(list(ordered(fterms)), key=lambda x: x.count(func))


def _linab(arg, symbol):
    """Return ``a, b, X`` assuming ``arg`` can be written as ``a*X + b``
    where ``X`` is a symbol-dependent factor and ``a`` and ``b`` are
    independent of ``symbol``.

    Examples
    ========

    >>> _linab(Integer(2), x)
    (2, 0, 1)
    >>> _linab(2*x, x)
    (2, 0, x)
    >>> _linab(y + y*x + 2*x, x)
    (y + 2, y, x)
    >>> _linab(3 + 2*exp(x), x)
    (2, 3, E**x)

    """
    arg = arg.expand()
    ind, dep = arg.as_independent(symbol)
    if not arg.is_Add:
        b = 0
        a, x = ind, dep
    else:
        b = ind
        a, x = separatevars(dep).as_independent(symbol, as_Add=False)
    if x.could_extract_minus_sign():
        a = -a
        x = -x
    return a, b, x


def _lambert(eq, x):
    """
    Given an expression assumed to be in the form
        ``F(X, a..f) = a*log(b*X + c) + d*X + f = 0``
    where X = g(x) and x = g^-1(X), return the Lambert solution if possible:
        ``x = g^-1(-c/b + (a/d)*W(d/(a*b)*exp(c*d/a/b)*exp(-f/a)))``.

    """
    eq = _mexpand(expand_log(eq))
    mainlog = _mostfunc(eq, log, x)
    if not mainlog:
        return []  # violated assumptions
    other = eq.subs({mainlog: 0})
    if isinstance(-other, log):
        eq = (eq - other).subs({mainlog: mainlog.args[0]})
        mainlog = mainlog.args[0]
        if not isinstance(mainlog, log):
            return []  # violated assumptions
        other = -(-other).args[0]
        eq += other
    if x not in other.free_symbols:
        return []  # violated assumptions
    d, f, X2 = _linab(other, x)
    logterm = collect(eq - other, mainlog)
    a = logterm.as_coefficient(mainlog)
    if a is None or x in a.free_symbols:
        return []  # violated assumptions
    logarg = mainlog.args[0]
    b, c, X1 = _linab(logarg, x)
    if X1 != X2:
        return []  # violated assumptions

    u = Dummy('rhs')
    sol = []
    # check only real solutions:
    for k in [-1, 0]:
        l = LambertW(d/(a*b)*exp(c*d/a/b)*exp(-f/a), k)
        # if W's arg is between -1/e and 0 there is
        # a -1 branch real solution, too.
        if k and not l.is_extended_real:
            continue
        rhs = -c/b + (a/d)*l

        solns = solve(X1 - u, x)
        for i, tmp in enumerate(solns):
            solns[i] = tmp[x].subs({u: rhs})
            sol.append(solns[i])
    return sol


def _solve_lambert(f, symbol, gens):
    """Return solution to ``f`` if it is a Lambert-type expression
    else raise NotImplementedError.

    The equality, ``f(x, a..f) = a*log(b*X + c) + d*X - f = 0`` has the
    solution,  `X = -c/b + (a/d)*W(d/(a*b)*exp(c*d/a/b)*exp(f/a))`. There
    are a variety of forms for `f(X, a..f)` as enumerated below:

    1a1)
      if B**B = R for R not [0, 1] then
      log(B) + log(log(B)) = log(log(R))
      X = log(B), a = 1, b = 1, c = 0, d = 1, f = log(log(R))
    1a2)
      if B*(b*log(B) + c)**a = R then
      log(B) + a*log(b*log(B) + c) = log(R)
      X = log(B); d=1, f=log(R)
    1b)
      if a*log(b*B + c) + d*B = R then
      X = B, f = R
    2)
      if d*p**(a*B + g) - b*B = c then
      log(d) + (a*B + g)*log(p) - log(c + b*B) = 0
      a = -1, d = a*log(p), f = -log(d) - g*log(p)

    """
    nrhs, lhs = f.as_independent(symbol, as_Add=True)
    rhs = -nrhs

    lamcheck = [tmp for tmp in gens
                if (tmp.func in [exp, log] or
                    (tmp.is_Pow and symbol in tmp.exp.free_symbols))]
    if not lamcheck:
        raise NotImplementedError()

    if lhs.is_Mul:
        lhs = expand_log(log(lhs))
        rhs = log(rhs)

    lhs = factor(lhs, deep=True)
    # make sure we have inverted as completely as possible
    r = Dummy()
    i, lhs = _invert(lhs - r, symbol)
    rhs = i.xreplace({r: rhs})

    # For the first ones:
    # 1a1) B**B = R != 0 (when 0, there is only a solution if the base is 0,
    #                     but if it is, the exp is 0 and 0**0=1
    #                     comes back as B*log(B) = log(R)
    # 1a2) B*(a + b*log(B))**p = R or with monomial expanded or with whole
    #                              thing expanded comes back unchanged
    #     log(B) + p*log(a + b*log(B)) = log(R)
    #     lhs is Mul:
    #         expand log of both sides to give:
    #         log(B) + log(log(B)) = log(log(R))
    # 1b) d*log(a*B + b) + c*B = R
    #     lhs is Add:
    #         isolate c*B and expand log of both sides:
    #         log(c) + log(B) = log(R - d*log(a*B + b))

    soln = []
    mainlog = _mostfunc(lhs, log, symbol)
    if mainlog:
        if lhs.is_Mul and rhs != 0:
            soln = _lambert(log(lhs) - log(rhs), symbol)
        elif lhs.is_Add:
            other = lhs.subs({mainlog: 0})
            if other and not other.is_Add and [
                    tmp for tmp in other.atoms(Pow)
                    if symbol in tmp.free_symbols]:
                if not rhs:
                    diff = log(other) - log(other - lhs)
                else:
                    diff = log(lhs - other) - log(rhs - other)
                soln = _lambert(expand_log(diff), symbol)
            else:
                # it's ready to go
                soln = _lambert(lhs - rhs, symbol)

    # 2) d*p**(a*B + b) + c*B = R
    #     collect on main pow
    #     log(R - c*B) - a*B*log(p) = log(d) + b*log(p)

    if not soln:
        mainpow = _mostfunc(lhs, Pow, symbol)
        if mainpow and symbol in mainpow.exp.free_symbols:
            lhs = collect(lhs, mainpow)
            if lhs.is_Mul and rhs != 0:
                soln = _lambert(expand_log(log(lhs) - log(rhs)), symbol)
            elif lhs.is_Add:
                # move all but mainpow-containing term to rhs
                other = lhs.subs({mainpow: 0})
                mainterm = lhs - other
                rhs = rhs - other
                diff = log(mainterm) - log(rhs)
                soln = _lambert(expand_log(diff), symbol)

    if not soln:
        raise NotImplementedError(f'{f} does not appear to have a solution in '
                                  'terms of LambertW')

    return list(ordered(soln))


def bivariate_type(f, x, y, **kwargs):
    """Given an expression, f, 3 tests will be done to see what type
    of composite bivariate it might be, options for u(x, y) are::

        x*y
        x+y
        x*y+x
        x*y+y

    If it matches one of these types, ``u(x, y)``, ``P(u)`` and dummy
    variable ``u`` will be returned. Solving ``P(u)`` for ``u`` and
    equating the solutions to ``u(x, y)`` and then solving for ``x`` or
    ``y`` is equivalent to solving the original expression for ``x`` or
    ``y``. If ``x`` and ``y`` represent two functions in the same
    variable, e.g. ``x = g(t)`` and ``y = h(t)``, then if ``u(x, y) - p``
    can be solved for ``t`` then these represent the solutions to
    ``P(u) = 0`` when ``p`` are the solutions of ``P(u) = 0``.

    Only positive values of ``u`` are considered.

    Examples
    ========

    >>> eq = (x**2 - 3).subs({x: x + y})
    >>> bivariate_type(eq, x, y)
    (x + y, _u**2 - 3, _u)
    >>> uxy, pu, u = _
    >>> usol = solve(pu, u)
    >>> usol
    [{_u: sqrt(3)}]
    >>> [solve(uxy - s[u]) for s in solve(pu, u)]
    [[{x: -y + sqrt(3)}]]
    >>> all(eq.subs(s).equals(0) for sol in _ for s in sol)
    True

    """
    u = Dummy('u', positive=True)

    if kwargs.pop('first', True):
        p = Poly(f, x, y)
        f = p.as_expr()
        _x = Dummy()
        _y = Dummy()
        rv = bivariate_type(Poly(f.subs({x: _x, y: _y}), _x, _y), _x, _y, first=False)
        if rv:
            reps = {_x: x, _y: y}
            return rv[0].xreplace(reps), rv[1].xreplace(reps), rv[2]
        return

    p = f
    f = p.as_expr()

    # f(x*y)
    args = Add.make_args(p.as_expr())
    new = []
    for a in args:
        a = _mexpand(a.subs({x: u/y}))
        free = a.free_symbols
        if x in free or y in free:
            break
        new.append(a)
    else:
        return x*y, Add(*new), u

    def ok(f, v, c):
        new = _mexpand(f.subs({v: c}))
        free = new.free_symbols
        return None if (x in free or y in free) else new

    # f(a*x + b*y)
    new = []
    d = p.degree(x)
    if p.degree(y) == d:
        a = root(p.coeff_monomial(x**d), d)
        b = root(p.coeff_monomial(y**d), d)
        new = ok(f, x, (u - b*y)/a)
        if new is not None:
            return a*x + b*y, new, u

    # f(a*x*y + b*y)
    new = []
    d = p.degree(x)
    if p.degree(y) == d:
        for _ in range(2):
            a = root(p.coeff_monomial(x**d*y**d), d)
            b = root(p.coeff_monomial(y**d), d)
            new = ok(f, x, (u - b*y)/a/y)
            if new is not None:
                return a*x*y + b*y, new, u
            x, y = y, x
