code: def _solve_lambert(f, symbol, gens):
      """Return solution to ``f`` if it is a Lambert-type expression
      else raise NotImplementedError.
  
      For ``f(X, a..f) = a*log(b*X + c) + d*X - f = 0`` the solution
      for ``X`` is ``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 in [0, 1] (since those cases would already
        be solved before getting here) then log of both sides gives
        log(B) + log(log(B)) = log(log(R)) and
        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 of both sides gives
        log(B) + a*log(b*log(B) + c) = log(R) and
        X = log(B), d=1, f=log(R)
      1b)
        if a*log(b*B + c) + d*B = R and
        X = B, f = R
      2a)
        if (b*B + c)*exp(d*B + g) = R then log of both sides gives
        log(b*B + c) + d*B + g = log(R) and
        X = B, a = 1, f = log(R) - g
      2b)
        if g*exp(d*B + h) - b*B = c then the log form is
        log(g) + d*B + h - log(b*B + c) = 0 and
        X = B, a = -1, f = -h - log(g)
      3)
        if d*p**(a*B + g) - b*B = c then the log form is
        log(d) + (a*B + g)*log(p) - log(b*B + c) = 0 and
        X = B, a = -1, d = a*log(p), f = -log(d) - g*log(p)
      """
  
      def _solve_even_degree_expr(expr, t, symbol):
          """Return the unique solutions of equations derived from
          ``expr`` by replacing ``t`` with ``+/- symbol``.
  
          Parameters
          ==========
  
          expr : Expr
              The expression which includes a dummy variable t to be
              replaced with +symbol and -symbol.
  
          symbol : Symbol
              The symbol for which a solution is being sought.
  
          Returns
          =======
  
          List of unique solution of the two equations generated by
          replacing ``t`` with positive and negative ``symbol``.
  
          Notes
          =====
  
          If ``expr = 2*log(t) + x/2` then solutions for
          ``2*log(x) + x/2 = 0`` and ``2*log(-x) + x/2 = 0`` are
          returned by this function. Though this may seem
          counter-intuitive, one must note that the ``expr`` being
          solved here has been derived from a different expression. For
          an expression like ``eq = x**2*g(x) = 1``, if we take the
          log of both sides we obtain ``log(x**2) + log(g(x)) = 0``. If
          x is positive then this simplifies to
          ``2*log(x) + log(g(x)) = 0``; the Lambert-solving routines will
          return solutions for this, but we must also consider the
          solutions for  ``2*log(-x) + log(g(x))`` since those must also
          be a solution of ``eq`` which has the same value when the ``x``
          in ``x**2`` is negated. If `g(x)` does not have even powers of
          symbol then we do not want to replace the ``x`` there with
          ``-x``. So the role of the ``t`` in the expression received by
          this function is to mark where ``+/-x`` should be inserted
          before obtaining the Lambert solutions.
  
          """
          nlhs, plhs = [
              expr.xreplace({t: sgn * symbol}) for sgn in (-1, 1)
          ]
          sols = _solve_lambert(nlhs, symbol, gens)
          if plhs != nlhs:
              sols.extend(_solve_lambert(plhs, symbol, gens))
          # uniq is needed for a case like
          # 2*log(t) - log(-z**2) + log(z + log(x) + log(z))
          # where substituting t with +/-x gives all the same solution;
          # uniq, rather than list(set()), is used to maintain canonical
          # order
          return list(uniq(sols))
  
      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_Add or lhs.is_Mul:
          # replacing all even_degrees of symbol with dummy variable t
          # since these will need special handling; non-Add/Mul do not
          # need this handling
          t = Dummy("t", **symbol.assumptions0)
          lhs = lhs.replace(
              lambda i: i.is_Pow  # find symbol**even
              and i.base == symbol
              and i.exp.is_even,
              lambda i: t**i.exp,  # replace t**even
          )
  
          if lhs.is_Add and lhs.has(t):
              t_indep = lhs.subs(t, 0)
              t_term = lhs - t_indep
              _rhs = rhs - t_indep
              if (
                  not t_term.is_Add
                  and _rhs
                  and not (t_term.has(S.ComplexInfinity, S.NaN))
              ):
                  eq = expand_log(log(t_term) - log(_rhs))
                  return _solve_even_degree_expr(eq, t, symbol)
          elif lhs.is_Mul and rhs:
              # this needs to happen whether t is present or not
              lhs = expand_log(log(lhs), force=True)
              rhs = log(rhs)
              if lhs.has(t) and lhs.is_Add:
                  # it expanded from Mul to Add
                  eq = lhs - rhs
                  return _solve_even_degree_expr(eq, t, symbol)
  
          # restore symbol in lhs
          lhs = lhs.xreplace({t: symbol})
  
      lhs = powsimp(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 forms:
      #
      # 1a1) B**B = R will arrive here as B*log(B) = log(R)
      #      lhs is Mul so take log of both sides:
      #        log(B) + log(log(B)) = log(log(R))
      # 1a2) B*(b*log(B) + c)**a = R will arrive unchanged so
      #      lhs is Mul, so take log of both sides:
      #        log(B) + a*log(b*log(B) + c) = log(R)
      # 1b) d*log(a*B + b) + c*B = R will arrive unchanged so
      #      lhs is Add, so isolate c*B and expand log of both sides:
      #        log(c) + log(B) = log(R - d*log(a*B + b))
  
      soln = []
      if not 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)
  
      # For the next forms,
      #
      #     collect on main exp
      #     2a) (b*B + c)*exp(d*B + g) = R
      #         lhs is mul, so take log of both sides:
      #           log(b*B + c) + d*B = log(R) - g
      #     2b) g*exp(d*B + h) - b*B = R
      #         lhs is add, so add b*B to both sides,
      #         take the log of both sides and rearrange to give
      #           log(R + b*B) - d*B = log(g) + h
  
      if not soln:
          mainexp = _mostfunc(lhs, exp, symbol)
          if mainexp:
              lhs = collect(lhs, mainexp)
              if lhs.is_Mul and rhs != 0:
                  soln = _lambert(
                      expand_log(log(lhs) - log(rhs)), symbol
                  )
              elif lhs.is_Add:
                  # move all but mainexp-containing term to rhs
                  other = lhs.subs(mainexp, 0)
                  mainterm = lhs - other
                  rhs = rhs - other
                  if (
                      mainterm.could_extract_minus_sign()
                      and rhs.could_extract_minus_sign()
                  ):
                      mainterm *= -1
                      rhs *= -1
                  diff = log(mainterm) - log(rhs)
                  soln = _lambert(expand_log(diff), symbol)
  
      # For the last form:
      #
      #  3) d*p**(a*B + g) - b*B = c
      #     collect on main pow, add b*B to both sides,
      #     take log of both sides and rearrange to give
      #       a*B*log(p) - log(b*B + c) = -log(d) - g*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:
                  # b*B = 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(
              "%s does not appear to have a solution in "
              "terms of LambertW" % f
          )
  
      return list(ordered(soln))
  
imports: from sympy.core.add import Add
  from sympy.core.exprtools import factor_terms
  from sympy.core.function import expand_log, _mexpand
  from sympy.core.power import Pow
  from sympy.core.singleton import S
  from sympy.core.sorting import ordered
  from sympy.core.symbol import Dummy
  from sympy.functions.elementary.exponential import LambertW, exp, log
  from sympy.functions.elementary.miscellaneous import root
  from sympy.polys.polyroots import roots
  from sympy.polys.polytools import Poly, factor
  from sympy.simplify.simplify import separatevars
  from sympy.simplify.radsimp import collect
  from sympy.simplify.simplify import powsimp
  from sympy.solvers.solvers import solve, _invert
  from sympy.utilities.iterables import uniq
label_desc: Use image node_id 5 for calling a global function with example usage: _solve_lambert(f, symbol, gens) and returns: list, list, _solve_even_degree_expr, _solve_even_degree_expr