code: def unrad(eq, *syms, **flags):
      """
      Remove radicals with symbolic arguments and return (eq, cov),
      None, or raise an error.
  
      Explanation
      ===========
  
      None is returned if there are no radicals to remove.
  
      NotImplementedError is raised if there are radicals and they cannot be
      removed or if the relationship between the original symbols and the
      change of variable needed to rewrite the system as a polynomial cannot
      be solved.
  
      Otherwise the tuple, ``(eq, cov)``, is returned where:
  
      *eq*, ``cov``
          *eq* is an equation without radicals (in the symbol(s) of
          interest) whose solutions are a superset of the solutions to the
          original expression. *eq* might be rewritten in terms of a new
          variable; the relationship to the original variables is given by
          ``cov`` which is a list containing ``v`` and ``v**p - b`` where
          ``p`` is the power needed to clear the radical and ``b`` is the
          radical now expressed as a polynomial in the symbols of interest.
          For example, for sqrt(2 - x) the tuple would be
          ``(c, c**2 - 2 + x)``. The solutions of *eq* will contain
          solutions to the original equation (if there are any).
  
      *syms*
          An iterable of symbols which, if provided, will limit the focus of
          radical removal: only radicals with one or more of the symbols of
          interest will be cleared. All free symbols are used if *syms* is not
          set.
  
      *flags* are used internally for communication during recursive calls.
      Two options are also recognized:
  
          ``take``, when defined, is interpreted as a single-argument function
          that returns True if a given Pow should be handled.
  
      Radicals can be removed from an expression if:
  
          *   All bases of the radicals are the same; a change of variables is
              done in this case.
          *   If all radicals appear in one term of the expression.
          *   There are only four terms with sqrt() factors or there are less than
              four terms having sqrt() factors.
          *   There are only two terms with radicals.
  
      Examples
      ========
  
      >>> from sympy.solvers.solvers import unrad
      >>> from sympy.abc import x
      >>> from sympy import sqrt, Rational, root
  
      >>> unrad(sqrt(x)*x**Rational(1, 3) + 2)
      (x**5 - 64, [])
      >>> unrad(sqrt(x) + root(x + 1, 3))
      (-x**3 + x**2 + 2*x + 1, [])
      >>> eq = sqrt(x) + root(x, 3) - 2
      >>> unrad(eq)
      (_p**3 + _p**2 - 2, [_p, _p**6 - x])
  
      """
  
      uflags = {"check": False, "simplify": False}
  
      def _cov(p, e):
          if cov:
              # XXX - uncovered
              oldp, olde = cov
              if Poly(e, p).degree(p) in (1, 2):
                  cov[:] = [
                      p,
                      olde.subs(oldp, _vsolve(e, p, **uflags)[0]),
                  ]
              else:
                  raise NotImplementedError
          else:
              cov[:] = [p, e]
  
      def _canonical(eq, cov):
          if cov:
              # change symbol to vanilla so no solutions are eliminated
              p, e = cov
              rep = {p: Dummy(p.name)}
              eq = eq.xreplace(rep)
              cov = [p.xreplace(rep), e.xreplace(rep)]
  
          # remove constants and powers of factors since these don't change
          # the location of the root; XXX should factor or factor_terms be used?
          eq = factor_terms(
              _mexpand(eq.as_numer_denom()[0], recursive=True),
              clear=True,
          )
          if eq.is_Mul:
              args = []
              for f in eq.args:
                  if f.is_number:
                      continue
                  if f.is_Pow:
                      args.append(f.base)
                  else:
                      args.append(f)
              eq = Mul(*args)  # leave as Mul for more efficient solving
  
          # make the sign canonical
          margs = list(Mul.make_args(eq))
          changed = False
          for i, m in enumerate(margs):
              if m.could_extract_minus_sign():
                  margs[i] = -m
                  changed = True
          if changed:
              eq = Mul(*margs, evaluate=False)
  
          return eq, cov
  
      def _Q(pow):
          # return leading Rational of denominator of Pow's exponent
          c = pow.as_base_exp()[1].as_coeff_Mul()[0]
          if not c.is_Rational:
              return S.One
          return c.q
  
      # define the _take method that will determine whether a term is of interest
      def _take(d):
          # return True if coefficient of any factor's exponent's den is not 1
          for pow in Mul.make_args(d):
              if not pow.is_Pow:
                  continue
              if _Q(pow) == 1:
                  continue
              if pow.free_symbols & syms:
                  return True
          return False
  
      _take = flags.setdefault("_take", _take)
  
      if isinstance(eq, Eq):
          eq = eq.lhs - eq.rhs  # XXX legacy Eq as Eqn support
      elif not isinstance(eq, Expr):
          return
  
      cov, nwas, rpt = [
          flags.setdefault(k, v)
          for k, v in sorted({"cov": [], "n": None, "rpt": 0}.items())
      ]
  
      # preconditioning
      eq = powdenest(factor_terms(eq, radical=True, clear=True))
      eq = eq.as_numer_denom()[0]
      eq = _mexpand(eq, recursive=True)
      if eq.is_number:
          return
  
      # see if there are radicals in symbols of interest
      syms = set(syms) or eq.free_symbols  # _take uses this
      poly = eq.as_poly()
      gens = [g for g in poly.gens if _take(g)]
      if not gens:
          return
  
      # recast poly in terms of eigen-gens
      poly = eq.as_poly(*gens)
  
      # not a polynomial e.g. 1 + sqrt(x)*exp(sqrt(x)) with gen sqrt(x)
      if poly is None:
          return
  
      # - an exponent has a symbol of interest (don't handle)
      if any(g.exp.has(*syms) for g in gens):
          return
  
      def _rads_bases_lcm(poly):
          # if all the bases are the same or all the radicals are in one
          # term, `lcm` will be the lcm of the denominators of the
          # exponents of the radicals
          lcm = 1
          rads = set()
          bases = set()
          for g in poly.gens:
              q = _Q(g)
              if q != 1:
                  rads.add(g)
                  lcm = ilcm(lcm, q)
                  bases.add(g.base)
          return rads, bases, lcm
  
      rads, bases, lcm = _rads_bases_lcm(poly)
  
      covsym = Dummy("p", nonnegative=True)
  
      # only keep in syms symbols that actually appear in radicals;
      # and update gens
      newsyms = set()
      for r in rads:
          newsyms.update(syms & r.free_symbols)
      if newsyms != syms:
          syms = newsyms
      # get terms together that have common generators
      drad = dict(zip(rads, range(len(rads))))
      rterms = {(): []}
      args = Add.make_args(poly.as_expr())
      for t in args:
          if _take(t):
              common = set(t.as_poly().gens).intersection(rads)
              key = tuple(sorted([drad[i] for i in common]))
          else:
              key = ()
          rterms.setdefault(key, []).append(t)
      others = Add(*rterms.pop(()))
      rterms = [Add(*rterms[k]) for k in rterms.keys()]
  
      # the output will depend on the order terms are processed, so
      # make it canonical quickly
      rterms = list(reversed(list(ordered(rterms))))
  
      ok = False  # we don't have a solution yet
      depth = sqrt_depth(eq)
  
      if len(rterms) == 1 and not (rterms[0].is_Add and lcm > 2):
          eq = rterms[0] ** lcm - ((-others) ** lcm)
          ok = True
      else:
          if len(rterms) == 1 and rterms[0].is_Add:
              rterms = list(rterms[0].args)
          if len(bases) == 1:
              b = bases.pop()
              if len(syms) > 1:
                  x = b.free_symbols
              else:
                  x = syms
              x = list(ordered(x))[0]
              try:
                  inv = _vsolve(covsym**lcm - b, x, **uflags)
                  if not inv:
                      raise NotImplementedError
                  eq = (
                      poly.as_expr()
                      .subs(b, covsym**lcm)
                      .subs(x, inv[0])
                  )
                  _cov(covsym, covsym**lcm - b)
                  return _canonical(eq, cov)
              except NotImplementedError:
                  pass
  
          if len(rterms) == 2:
              if not others:
                  eq = rterms[0] ** lcm - (-rterms[1]) ** lcm
                  ok = True
              elif not log(lcm, 2).is_Integer:
                  # the lcm-is-power-of-two case is handled below
                  r0, r1 = rterms
                  if flags.get("_reverse", False):
                      r1, r0 = r0, r1
                  i0 = _rads0, _bases0, lcm0 = _rads_bases_lcm(
                      r0.as_poly()
                  )
                  i1 = _rads1, _bases1, lcm1 = _rads_bases_lcm(
                      r1.as_poly()
                  )
                  for reverse in range(2):
                      if reverse:
                          i0, i1 = i1, i0
                          r0, r1 = r1, r0
                      _rads1, _, lcm1 = i1
                      _rads1 = Mul(*_rads1)
                      t1 = _rads1**lcm1
                      c = covsym**lcm1 - t1
                      for x in syms:
                          try:
                              sol = _vsolve(c, x, **uflags)
                              if not sol:
                                  raise NotImplementedError
                              neweq = (
                                  r0.subs(x, sol[0])
                                  + covsym * r1 / _rads1
                                  + others
                              )
                              tmp = unrad(neweq, covsym)
                              if tmp:
                                  eq, newcov = tmp
                                  if newcov:
                                      newp, newc = newcov
                                      _cov(
                                          newp,
                                          c.subs(
                                              covsym,
                                              _vsolve(
                                                  newc, covsym, **uflags
                                              )[0],
                                          ),
                                      )
                                  else:
                                      _cov(covsym, c)
                              else:
                                  eq = neweq
                                  _cov(covsym, c)
                              ok = True
                              break
                          except NotImplementedError:
                              if reverse:
                                  raise NotImplementedError(
                                      "no successful change of variable found"
                                  )
                              else:
                                  pass
                      if ok:
                          break
          elif len(rterms) == 3:
              # two cube roots and another with order less than 5
              # (so an analytical solution can be found) or a base
              # that matches one of the cube root bases
              info = [_rads_bases_lcm(i.as_poly()) for i in rterms]
              RAD = 0
              BASES = 1
              LCM = 2
              if info[0][LCM] != 3:
                  info.append(info.pop(0))
                  rterms.append(rterms.pop(0))
              elif info[1][LCM] != 3:
                  info.append(info.pop(1))
                  rterms.append(rterms.pop(1))
              if info[0][LCM] == info[1][LCM] == 3:
                  if info[1][BASES] != info[2][BASES]:
                      info[0], info[1] = info[1], info[0]
                      rterms[0], rterms[1] = rterms[1], rterms[0]
                  if info[1][BASES] == info[2][BASES]:
                      eq = (
                          rterms[0] ** 3
                          + (rterms[1] + rterms[2] + others) ** 3
                      )
                      ok = True
                  elif info[2][LCM] < 5:
                      # a*root(A, 3) + b*root(B, 3) + others = c
                      a, b, c, d, A, B = [Dummy(i) for i in "abcdAB"]
                      # zz represents the unraded expression into which the
                      # specifics for this case are substituted
                      zz = (c - d) * (
                          A**3 * a**9
                          + 3 * A**2 * B * a**6 * b**3
                          - 3 * A**2 * a**6 * c**3
                          + 9 * A**2 * a**6 * c**2 * d
                          - 9 * A**2 * a**6 * c * d**2
                          + 3 * A**2 * a**6 * d**3
                          + 3 * A * B**2 * a**3 * b**6
                          + 21 * A * B * a**3 * b**3 * c**3
                          - 63 * A * B * a**3 * b**3 * c**2 * d
                          + 63 * A * B * a**3 * b**3 * c * d**2
                          - 21 * A * B * a**3 * b**3 * d**3
                          + 3 * A * a**3 * c**6
                          - 18 * A * a**3 * c**5 * d
                          + 45 * A * a**3 * c**4 * d**2
                          - 60 * A * a**3 * c**3 * d**3
                          + 45 * A * a**3 * c**2 * d**4
                          - 18 * A * a**3 * c * d**5
                          + 3 * A * a**3 * d**6
                          + B**3 * b**9
                          - 3 * B**2 * b**6 * c**3
                          + 9 * B**2 * b**6 * c**2 * d
                          - 9 * B**2 * b**6 * c * d**2
                          + 3 * B**2 * b**6 * d**3
                          + 3 * B * b**3 * c**6
                          - 18 * B * b**3 * c**5 * d
                          + 45 * B * b**3 * c**4 * d**2
                          - 60 * B * b**3 * c**3 * d**3
                          + 45 * B * b**3 * c**2 * d**4
                          - 18 * B * b**3 * c * d**5
                          + 3 * B * b**3 * d**6
                          - c**9
                          + 9 * c**8 * d
                          - 36 * c**7 * d**2
                          + 84 * c**6 * d**3
                          - 126 * c**5 * d**4
                          + 126 * c**4 * d**5
                          - 84 * c**3 * d**6
                          + 36 * c**2 * d**7
                          - 9 * c * d**8
                          + d**9
                      )
  
                      def _t(i):
                          b = Mul(*info[i][RAD])
                          return cancel(rterms[i] / b), Mul(
                              *info[i][BASES]
                          )
  
                      aa, AA = _t(0)
                      bb, BB = _t(1)
                      cc = -rterms[2]
                      dd = others
                      eq = zz.xreplace(
                          dict(
                              zip(
                                  (a, A, b, B, c, d),
                                  (aa, AA, bb, BB, cc, dd),
                              )
                          )
                      )
                      ok = True
          # handle power-of-2 cases
          if not ok:
              if log(lcm, 2).is_Integer and (
                  not others and len(rterms) == 4 or len(rterms) < 4
              ):
  
                  def _norm2(a, b):
                      return a**2 + b**2 + 2 * a * b
  
                  if len(rterms) == 4:
                      # (r0+r1)**2 - (r2+r3)**2
                      r0, r1, r2, r3 = rterms
                      eq = _norm2(r0, r1) - _norm2(r2, r3)
                      ok = True
                  elif len(rterms) == 3:
                      # (r1+r2)**2 - (r0+others)**2
                      r0, r1, r2 = rterms
                      eq = _norm2(r1, r2) - _norm2(r0, others)
                      ok = True
                  elif len(rterms) == 2:
                      # r0**2 - (r1+others)**2
                      r0, r1 = rterms
                      eq = r0**2 - _norm2(r1, others)
                      ok = True
  
      new_depth = sqrt_depth(eq) if ok else depth
      rpt += 1  # XXX how many repeats with others unchanging is enough?
      if not ok or (
          nwas is not None
          and len(rterms) == nwas
          and new_depth is not None
          and new_depth == depth
          and rpt > 3
      ):
          raise NotImplementedError("Cannot remove all radicals")
  
      flags.update({"cov": cov, "n": len(rterms), "rpt": rpt})
      neq = unrad(eq, *syms, **flags)
      if neq:
          eq, cov = neq
      eq, cov = _canonical(eq, cov)
      return eq, cov
  
imports: from __future__ import annotations
  from sympy.core import S, Add, Symbol, Dummy, Expr, Mul
  from sympy.core.assumptions import check_assumptions
  from sympy.core.exprtools import factor_terms
  from sympy.core.function import expand_mul, expand_log, Derivative, AppliedUndef, UndefinedFunction, nfloat, Function, expand_power_exp, _mexpand, expand, expand_func
  from sympy.core.logic import fuzzy_not
  from sympy.core.numbers import Float, Rational, _illegal
  from sympy.core.intfunc import integer_log, ilcm
  from sympy.core.power import Pow
  from sympy.core.relational import Eq, Ne
  from sympy.core.sorting import ordered, default_sort_key
  from sympy.core.sympify import sympify, _sympify
  from sympy.core.traversal import preorder_traversal
  from sympy.logic.boolalg import And, BooleanAtom
  from sympy.functions import log, exp, LambertW, cos, sin, tan, acos, asin, atan, Abs, re, im, arg, sqrt, atan2
  from sympy.functions.combinatorial.factorials import binomial
  from sympy.functions.elementary.hyperbolic import HyperbolicFunction
  from sympy.functions.elementary.piecewise import piecewise_fold, Piecewise
  from sympy.functions.elementary.trigonometric import TrigonometricFunction
  from sympy.integrals.integrals import Integral
  from sympy.ntheory.factor_ import divisors
  from sympy.simplify import simplify, collect, powsimp, posify, powdenest, nsimplify, denom, logcombine, sqrtdenest, fraction, separatevars
  from sympy.simplify.sqrtdenest import sqrt_depth
  from sympy.simplify.fu import TR1, TR2i, TR10, TR11
  from sympy.strategies.rl import rebuild
  from sympy.matrices.exceptions import NonInvertibleMatrixError
  from sympy.matrices import Matrix, zeros
  from sympy.polys import roots, cancel, factor, Poly
  from sympy.polys.solvers import sympy_eqs_to_ring, solve_lin_sys
  from sympy.polys.polyerrors import GeneratorsNeeded, PolynomialError
  from sympy.polys.polytools import gcd
  from sympy.utilities.lambdify import lambdify
  from sympy.utilities.misc import filldedent, debugf
  from sympy.utilities.iterables import connected_components, generate_bell, uniq, iterable, is_sequence, subsets, flatten, sift
  from sympy.utilities.decorator import conserve_mpmath_dps
  from mpmath import findroot
  from sympy.solvers.polysys import solve_poly_system
  from types import GeneratorType
  from collections import defaultdict
  from itertools import combinations, product
  import warnings
  from sympy.solvers.bivariate import bivariate_type, _solve_lambert, _filtered_gens
label_desc: Use image node_id 24 for calling a global function with example usage: unrad(eq) and returns: eq,  cov, eq,  cov, c, False, rads,  bases,  lcm, S, True, _canonical, unknown, cancel,  Mul