text: """Prime ideals in number fields. """
  
  from sympy.core.expr import Expr
  from sympy.polys.polytools import Poly
  from sympy.polys.domains.finitefield import FF
  from sympy.polys.domains.rationalfield import QQ
  from sympy.polys.domains.integerring import ZZ
  from sympy.polys.matrices.domainmatrix import DomainMatrix
  from sympy.polys.polyerrors import CoercionFailed, GeneratorsNeeded
  from sympy.polys.polyutils import IntegerPowerable
  from sympy.utilities.decorator import public
  from .basis import round_two, nilradical_mod_p
  from .exceptions import StructureError
  from .modules import ModuleEndomorphism, find_min_poly
  from .utilities import coeff_search, supplement_a_subspace
  
  
  def _check_formal_conditions_for_maximal_order(submodule):
      r"""
      Several functions in this module accept an argument which is to be a
      :py:class:`~.Submodule` representing the maximal order in a number field,
      such as returned by the :py:func:`~sympy.polys.numberfields.basis.round_two`
      algorithm.
  
      We do not attempt to check that the given ``Submodule`` actually represents
      a maximal order, but we do check a basic set of formal conditions that the
      ``Submodule`` must satisfy, at a minimum. The purpose is to catch an
      obviously ill-formed argument.
      """
      prefix = 'The submodule representing the maximal order should '
      cond = None
      if not submodule.is_power_basis_submodule():
          cond = 'be a direct submodule of a power basis.'
      elif not submodule.starts_with_unity():
          cond = 'have 1 as its first generator.'
      elif not submodule.is_sq_maxrank_HNF():
          cond = 'have square matrix, of maximal rank, in Hermite Normal Form.'
      if cond is not None:
          raise StructureError(prefix + cond)
  
  
  class PrimeIdeal(IntegerPowerable):
      r"""
      A prime ideal in a ring of algebraic integers.
      """
  
      def __init__(self, ZK, p, alpha, f, e=None):
          """
          Parameters
          ==========
  
          ZK : :py:class:`~.Submodule`
              The maximal order where this ideal lives.
          p : int
              The rational prime this ideal divides.
          alpha : :py:class:`~.PowerBasisElement`
              Such that the ideal is equal to ``p*ZK + alpha*ZK``.
          f : int
              The inertia degree.
          e : int, ``None``, optional
              The ramification index, if already known. If ``None``, we will
              compute it here.
  
          """
          _check_formal_conditions_for_maximal_order(ZK)
          self.ZK = ZK
          self.p = p
          self.alpha = alpha
          self.f = f
          self._test_factor = None
          self.e = e if e is not None else self.valuation(p * ZK)
  
      def _pretty(self, field_gen=None, just_gens=False):
          """
          Print a representation of this prime ideal.
  
          Examples
          ========
  
          >>> from sympy import cyclotomic_poly, QQ
          >>> from sympy.abc import x, zeta
          >>> T = cyclotomic_poly(7, x)
          >>> K = QQ.algebraic_field((T, zeta))
          >>> P = K.primes_above(11)
          >>> print(P[0]._pretty())
          [ (11, x**3 + 5*x**2 + 4*x - 1) e=1, f=3 ]
          >>> print(P[0]._pretty(field_gen=zeta))
          [ (11, zeta**3 + 5*zeta**2 + 4*zeta - 1) e=1, f=3 ]
          >>> print(P[0]._pretty(field_gen=zeta, just_gens=True))
          (11, zeta**3 + 5*zeta**2 + 4*zeta - 1)
  
          Parameters
          ==========
  
          field_gen : :py:class:`~.Symbol`, ``None``, optional (default=None)
              The symbol to use for the generator of the field. This will appear
              in our representation of ``self.alpha``. If ``None``, we use the
              variable of the defining polynomial of ``self.ZK``.
          just_gens : bool, optional (default=False)
              If ``True``, just print the "(p, alpha)" part, showing "just the
              generators" of the prime ideal. Otherwise, print a string of the
              form "[ (p, alpha) e=..., f=... ]", giving the ramification index
              and inertia degree, along with the generators.
  
          """
          field_gen = field_gen or self.ZK.parent.T.gen
          p, alpha, e, f = self.p, self.alpha, self.e, self.f
          alpha_rep = str(alpha.numerator(x=field_gen).as_expr())
          if alpha.denom > 1:
              alpha_rep = f'({alpha_rep})/{alpha.denom}'
          gens = f'({p}, {alpha_rep})'
          if just_gens:
              return gens
          return f'[ {gens} e={e}, f={f} ]'
  
      def __repr__(self):
          return self._pretty()
  
      def as_submodule(self):
          r"""
          Represent this prime ideal as a :py:class:`~.Submodule`.
  
          Explanation
          ===========
  
          The :py:class:`~.PrimeIdeal` class serves to bundle information about
          a prime ideal, such as its inertia degree, ramification index, and
          two-generator representation, as well as to offer helpful methods like
          :py:meth:`~.PrimeIdeal.valuation` and
          :py:meth:`~.PrimeIdeal.test_factor`.
  
          However, in order to be added and multiplied by other ideals or
          rational numbers, it must first be converted into a
          :py:class:`~.Submodule`, which is a class that supports these
          operations.
  
          In many cases, the user need not perform this conversion deliberately,
          since it is automatically performed by the arithmetic operator methods
          :py:meth:`~.PrimeIdeal.__add__` and :py:meth:`~.PrimeIdeal.__mul__`.
  
          Raising a :py:class:`~.PrimeIdeal` to a non-negative integer power is
          also supported.
  
          Examples
          ========
  
          >>> from sympy import Poly, cyclotomic_poly, prime_decomp
          >>> T = Poly(cyclotomic_poly(7))
          >>> P0 = prime_decomp(7, T)[0]
          >>> print(P0**6 == 7*P0.ZK)
          True
  
          Note that, on both sides of the equation above, we had a
          :py:class:`~.Submodule`. In the next equation we recall that adding
          ideals yields their GCD. This time, we need a deliberate conversion
          to :py:class:`~.Submodule` on the right:
  
          >>> print(P0 + 7*P0.ZK == P0.as_submodule())
          True
  
          Returns
          =======
  
          :py:class:`~.Submodule`
              Will be equal to ``self.p * self.ZK + self.alpha * self.ZK``.
  
          See Also
          ========
  
          __add__
          __mul__
  
          """
          M = self.p * self.ZK + self.alpha * self.ZK
          # Pre-set expensive boolean properties whose value we already know:
          M._starts_with_unity = False
          M._is_sq_maxrank_HNF = True
          return M
  
      def __eq__(self, other):
          if isinstance(other, PrimeIdeal):
              return self.as_submodule() == other.as_submodule()
          return NotImplemented
  
      def __add__(self, other):
          """
          Convert to a :py:class:`~.Submodule` and add to another
          :py:class:`~.Submodule`.
  
          See Also
          ========
  
          as_submodule
  
          """
          return self.as_submodule() + other
  
      __radd__ = __add__
  
      def __mul__(self, other):
          """
          Convert to a :py:class:`~.Submodule` and multiply by another
          :py:class:`~.Submodule` or a rational number.
  
          See Also
          ========
  
          as_submodule
  
          """
          return self.as_submodule() * other
  
      __rmul__ = __mul__
  
      def _zeroth_power(self):
          return self.ZK
  
      def _first_power(self):
          return self
  
      def test_factor(self):
          r"""
          Compute a test factor for this prime ideal.
  
          Explanation
          ===========
  
          Write $\mathfrak{p}$ for this prime ideal, $p$ for the rational prime
          it divides. Then, for computing $\mathfrak{p}$-adic valuations it is
          useful to have a number $\beta \in \mathbb{Z}_K$ such that
          $p/\mathfrak{p} = p \mathbb{Z}_K + \beta \mathbb{Z}_K$.
  
          Essentially, this is the same as the number $\Psi$ (or the "reagent")
          from Kummer's 1847 paper (*Ueber die Zerlegung...*, Crelle vol. 35) in
          which ideal divisors were invented.
          """
          if self._test_factor is None:
              self._test_factor = _compute_test_factor(self.p, [self.alpha], self.ZK)
          return self._test_factor
  
      def valuation(self, I):
          r"""
          Compute the $\mathfrak{p}$-adic valuation of integral ideal I at this
          prime ideal.
  
          Parameters
          ==========
  
          I : :py:class:`~.Submodule`
  
          See Also
          ========
  
          prime_valuation
  
          """
          return prime_valuation(I, self)
  
      def _reduce_poly(self, f, gen=None):
          r"""
          Reduce a univariate :py:class:`~.Poly` *f*, or an :py:class:`~.Expr`
          expressing the same, modulo this :py:class:`~.PrimeIdeal`.
  
          Explanation
          ===========
  
          If our second generator $\alpha$ is zero, then we simply reduce the
          coefficients of *f* mod the rational prime $p$ lying under this ideal.
  
          Otherwise we first reduce *f* mod $\alpha$ (as a polynomial in the same
          variable as *f*), and then mod $p$.
  
          Examples
          ========
  
          >>> from sympy import QQ, cyclotomic_poly, symbols
          >>> zeta = symbols('zeta')
          >>> Phi = cyclotomic_poly(7, zeta)
          >>> k = QQ.algebraic_field((Phi, zeta))
          >>> P = k.primes_above(11)
          >>> frp = P[0]
          >>> B = k.integral_basis(fmt='sympy')
          >>> print([frp._reduce_poly(b, zeta) for b in B])
          [1, zeta, zeta**2, -5*zeta**2 - 4*zeta + 1, -zeta**2 - zeta - 5,
           4*zeta**2 - zeta - 1]
  
          Parameters
          ==========
  
          f : :py:class:`~.Poly`, :py:class:`~.Expr`
              The univariate polynomial to be reduced.
  
          gen : :py:class:`~.Symbol`, None, optional (default=None)
              Symbol to use as the variable in the polynomials. If *f* is a
              :py:class:`~.Poly` or a non-constant :py:class:`~.Expr`, this
              replaces its variable. If *f* is a constant :py:class:`~.Expr`,
              then *gen* must be supplied.
  
          Returns
          =======
  
          :py:class:`~.Poly`, :py:class:`~.Expr`
              Type is same as that of given *f*. If returning a
              :py:class:`~.Poly`, its domain will be the finite field
              $\mathbb{F}_p$.
  
          Raises
          ======
  
          GeneratorsNeeded
              If *f* is a constant :py:class:`~.Expr` and *gen* is ``None``.
          NotImplementedError
              If *f* is other than :py:class:`~.Poly` or :py:class:`~.Expr`,
              or is not univariate.
  
          """
          if isinstance(f, Expr):
              try:
                  g = Poly(f)
              except GeneratorsNeeded as e:
                  if gen is None:
                      raise e from None
                  g = Poly(f, gen)
              return self._reduce_poly(g).as_expr()
          if isinstance(f, Poly) and f.is_univariate:
              a = self.alpha.poly(f.gen)
              if a != 0:
                  f = f.rem(a)
              return f.set_modulus(self.p)
          raise NotImplementedError
  
  
  def _compute_test_factor(p, gens, ZK):
      r"""
      Compute the test factor for a :py:class:`~.PrimeIdeal` $\mathfrak{p}$.
  
      Parameters
      ==========
  
      p : int
          The rational prime $\mathfrak{p}$ divides
  
      gens : list of :py:class:`PowerBasisElement`
          A complete set of generators for $\mathfrak{p}$ over *ZK*, EXCEPT that
          an element equivalent to rational *p* can and should be omitted (since
          it has no effect except to waste time).
  
      ZK : :py:class:`~.Submodule`
          The maximal order where the prime ideal $\mathfrak{p}$ lives.
  
      Returns
      =======
  
      :py:class:`~.PowerBasisElement`
  
      References
      ==========
  
      .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.*
      (See Proposition 4.8.15.)
  
      """
      _check_formal_conditions_for_maximal_order(ZK)
      E = ZK.endomorphism_ring()
      matrices = [E.inner_endomorphism(g).matrix(modulus=p) for g in gens]
      B = DomainMatrix.zeros((0, ZK.n), FF(p)).vstack(*matrices)
      # A nonzero element of the nullspace of B will represent a
      # lin comb over the omegas which (i) is not a multiple of p
      # (since it is nonzero over FF(p)), while (ii) is such that
      # its product with each g in gens _is_ a multiple of p (since
      # B represents multiplication by these generators). Theory
      # predicts that such an element must exist, so nullspace should
      # be non-trivial.
      x = B.nullspace()[0, :].transpose()
      beta = ZK.parent(ZK.matrix * x, denom=ZK.denom)
      return beta
  
  
  @public
  def prime_valuation(I, P):
      r"""
      Compute the *P*-adic valuation for an integral ideal *I*.
  
      Examples
      ========
  
      >>> from sympy import QQ
      >>> from sympy.abc import theta
      >>> from sympy.polys import cyclotomic_poly
      >>> from sympy.polys.numberfields import prime_valuation
      >>> T = cyclotomic_poly(5)
      >>> K = QQ.algebraic_field((T, theta))
      >>> P = K.primes_above(5)
      >>> ZK = K.maximal_order()
      >>> print(prime_valuation(25*ZK, P[0]))
      8
  
      Parameters
      ==========
  
      I : :py:class:`~.Submodule`
          An integral ideal whose valuation is desired.
  
      P : :py:class:`~.PrimeIdeal`
          The prime at which to compute the valuation.
  
      Returns
      =======
  
      int
  
      See Also
      ========
  
      .PrimeIdeal.valuation
  
      References
      ==========
  
      .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.*
         (See Algorithm 4.8.17.)
  
      """
      p, ZK = P.p, P.ZK
      n, W, d = ZK.n, ZK.matrix, ZK.denom
  
      A = W.convert_to(QQ).inv() * I.matrix * d / I.denom
      # Although A must have integer entries, given that I is an integral ideal,
      # as a DomainMatrix it will still be over QQ, so we convert back:
      A = A.convert_to(ZZ)
      D = A.det()
      if D % p != 0:
          return 0
  
      beta = P.test_factor()
  
      f = d ** n // W.det()
      need_complete_test = (f % p == 0)
      v = 0
      while True:
          # Entering the loop, the cols of A represent lin combs of omegas.
          # Turn them into lin combs of thetas:
          A = W * A
          # And then one column at a time...
          for j in range(n):
              c = ZK.parent(A[:, j], denom=d)
              c *= beta
              # ...turn back into lin combs of omegas, after multiplying by beta:
              c = ZK.represent(c).flat()
              for i in range(n):
                  A[i, j] = c[i]
          if A[n - 1, n - 1].element % p != 0:
              break
          A = A / p
          # As noted above, domain converts to QQ even when division goes evenly.
          # So must convert back, even when we don't "need_complete_test".
          if need_complete_test:
              # In this case, having a non-integer entry is actually just our
              # halting condition.
              try:
                  A = A.convert_to(ZZ)
              except CoercionFailed:
                  break
          else:
              # In this case theory says we should not have any non-integer entries.
              A = A.convert_to(ZZ)
          v += 1
      return v
  
  
  def _two_elt_rep(gens, ZK, p, f=None, Np=None):
      r"""
      Given a set of *ZK*-generators of a prime ideal, compute a set of just two
      *ZK*-generators for the same ideal, one of which is *p* itself.
  
      Parameters
      ==========
  
      gens : list of :py:class:`PowerBasisElement`
          Generators for the prime ideal over *ZK*, the ring of integers of the
          field $K$.
  
      ZK : :py:class:`~.Submodule`
          The maximal order in $K$.
  
      p : int
          The rational prime divided by the prime ideal.
  
      f : int, optional
          The inertia degree of the prime ideal, if known.
  
      Np : int, optional
          The norm $p^f$ of the prime ideal, if known.
          NOTE: There is no reason to supply both *f* and *Np*. Either one will
          save us from having to compute the norm *Np* ourselves. If both are known,
          *Np* is preferred since it saves one exponentiation.
  
      Returns
      =======
  
      :py:class:`~.PowerBasisElement` representing a single algebraic integer
      alpha such that the prime ideal is equal to ``p*ZK + alpha*ZK``.
  
      References
      ==========
  
      .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.*
      (See Algorithm 4.7.10.)
  
      """
      _check_formal_conditions_for_maximal_order(ZK)
      pb = ZK.parent
      T = pb.T
      # Detect the special cases in which either (a) all generators are multiples
      # of p, or (b) there are no generators (so `all` is vacuously true):
      if all((g % p).equiv(0) for g in gens):
          return pb.zero()
  
      if Np is None:
          if f is not None:
              Np = p**f
          else:
              Np = abs(pb.submodule_from_gens(gens).matrix.det())
  
      omega = ZK.basis_element_pullbacks()
      beta = [p*om for om in omega[1:]]  # note: we omit omega[0] == 1
      beta += gens
      search = coeff_search(len(beta), 1)
      for c in search:
          alpha = sum(ci*betai for ci, betai in zip(c, beta))
          # Note: It may be tempting to reduce alpha mod p here, to try to work
          # with smaller numbers, but must not do that, as it can result in an
          # infinite loop! E.g. try factoring 2 in Q(sqrt(-7)).
          n = alpha.norm(T) // Np
          if n % p != 0:
              # Now can reduce alpha mod p.
              return alpha % p
  
  
  def _prime_decomp_easy_case(p, ZK):
      r"""
      Compute the decomposition of rational prime *p* in the ring of integers
      *ZK* (given as a :py:class:`~.Submodule`), in the "easy case", i.e. the
      case where *p* does not divide the index of $\theta$ in *ZK*, where
      $\theta$ is the generator of the ``PowerBasis`` of which *ZK* is a
      ``Submodule``.
      """
      T = ZK.parent.T
      T_bar = Poly(T, modulus=p)
      lc, fl = T_bar.factor_list()
      return [PrimeIdeal(ZK, p,
                         ZK.parent.element_from_poly(Poly(t, domain=ZZ)),
                         t.degree(), e)
              for t, e in fl]
  
  
  def _prime_decomp_compute_kernel(I, p, ZK):
      r"""
      Parameters
      ==========
  
      I : :py:class:`~.Module`
          An ideal of ``ZK/pZK``.
      p : int
          The rational prime being factored.
      ZK : :py:class:`~.Submodule`
          The maximal order.
  
      Returns
      =======
  
      Pair ``(N, G)``, where:
  
          ``N`` is a :py:class:`~.Module` representing the kernel of the map
          ``a |--> a**p - a`` on ``(O/pO)/I``, guaranteed to be a module with
          unity.
  
          ``G`` is a :py:class:`~.Module` representing a basis for the separable
          algebra ``A = O/I`` (see Cohen).
  
      """
      W = I.matrix
      n, r = W.shape
      # Want to take the Fp-basis given by the columns of I, adjoin (1, 0, ..., 0)
      # (which we know is not already in there since I is a basis for a prime ideal)
      # and then supplement this with additional columns to make an invertible n x n
      # matrix. This will then represent a full basis for ZK, whose first r columns
      # are pullbacks of the basis for I.
      if r == 0:
          B = W.eye(n, ZZ)
      else:
          B = W.hstack(W.eye(n, ZZ)[:, 0])
      if B.shape[1] < n:
          B = supplement_a_subspace(B.convert_to(FF(p))).convert_to(ZZ)
  
      G = ZK.submodule_from_matrix(B)
      # Must compute G's multiplication table _before_ discarding the first r
      # columns. (See Step 9 in Alg 6.2.9 in Cohen, where the betas are actually
      # needed in order to represent each product of gammas. However, once we've
      # found the representations, then we can ignore the betas.)
      G.compute_mult_tab()
      G = G.discard_before(r)
  
      phi = ModuleEndomorphism(G, lambda x: x**p - x)
      N = phi.kernel(modulus=p)
      assert N.starts_with_unity()
      return N, G
  
  
  def _prime_decomp_maximal_ideal(I, p, ZK):
      r"""
      We have reached the case where we have a maximal (hence prime) ideal *I*,
      which we know because the quotient ``O/I`` is a field.
  
      Parameters
      ==========
  
      I : :py:class:`~.Module`
          An ideal of ``O/pO``.
      p : int
          The rational prime being factored.
      ZK : :py:class:`~.Submodule`
          The maximal order.
  
      Returns
      =======
  
      :py:class:`~.PrimeIdeal` instance representing this prime
  
      """
      m, n = I.matrix.shape
      f = m - n
      G = ZK.matrix * I.matrix
      gens = [ZK.parent(G[:, j], denom=ZK.denom) for j in range(G.shape[1])]
      alpha = _two_elt_rep(gens, ZK, p, f=f)
      return PrimeIdeal(ZK, p, alpha, f)
  
  
  def _prime_decomp_split_ideal(I, p, N, G, ZK):
      r"""
      Perform the step in the prime decomposition algorithm where we have determined
      the the quotient ``ZK/I`` is _not_ a field, and we want to perform a non-trivial
      factorization of *I* by locating an idempotent element of ``ZK/I``.
      """
      assert I.parent == ZK and G.parent is ZK and N.parent is G
      # Since ZK/I is not a field, the kernel computed in the previous step contains
      # more than just the prime field Fp, and our basis N for the nullspace therefore
      # contains at least a second column (which represents an element outside Fp).
      # Let alpha be such an element:
      alpha = N(1).to_parent()
      assert alpha.module is G
  
      alpha_powers = []
      m = find_min_poly(alpha, FF(p), powers=alpha_powers)
      # TODO (future work):
      #  We don't actually need full factorization, so might use a faster method
      #  to just break off a single non-constant factor m1?
      lc, fl = m.factor_list()
      m1 = fl[0][0]
      m2 = m.quo(m1)
      U, V, g = m1.gcdex(m2)
      # Sanity check: theory says m is squarefree, so m1, m2 should be coprime:
      assert g == 1
      E = list(reversed(Poly(U * m1, domain=ZZ).rep.rep))
      eps1 = sum(E[i]*alpha_powers[i] for i in range(len(E)))
      eps2 = 1 - eps1
      idemps = [eps1, eps2]
      factors = []
      for eps in idemps:
          e = eps.to_parent()
          assert e.module is ZK
          D = I.matrix.convert_to(FF(p)).hstack(*[
              (e * om).column(domain=FF(p)) for om in ZK.basis_elements()
          ])
          W = D.columnspace().convert_to(ZZ)
          H = ZK.submodule_from_matrix(W)
          factors.append(H)
      return factors
  
  
  @public
  def prime_decomp(p, T=None, ZK=None, dK=None, radical=None):
      r"""
      Compute the decomposition of rational prime *p* in a number field.
  
      Explanation
      ===========
  
      Ordinarily this should be accessed through the
      :py:meth:`~.AlgebraicField.primes_above` method of an
      :py:class:`~.AlgebraicField`.
  
      Examples
      ========
  
      >>> from sympy import Poly, QQ
      >>> from sympy.abc import x, theta
      >>> T = Poly(x ** 3 + x ** 2 - 2 * x + 8)
      >>> K = QQ.algebraic_field((T, theta))
      >>> print(K.primes_above(2))
      [[ (2, x**2 + 1) e=1, f=1 ], [ (2, (x**2 + 3*x + 2)/2) e=1, f=1 ],
       [ (2, (3*x**2 + 3*x)/2) e=1, f=1 ]]
  
      Parameters
      ==========
  
      p : int
          The rational prime whose decomposition is desired.
  
      T : :py:class:`~.Poly`, optional
          Monic irreducible polynomial defining the number field $K$ in which to
          factor. NOTE: at least one of *T* or *ZK* must be provided.
  
      ZK : :py:class:`~.Submodule`, optional
          The maximal order for $K$, if already known.
          NOTE: at least one of *T* or *ZK* must be provided.
  
      dK : int, optional
          The discriminant of the field $K$, if already known.
  
      radical : :py:class:`~.Submodule`, optional
          The nilradical mod *p* in the integers of $K$, if already known.
  
      Returns
      =======
  
      List of :py:class:`~.PrimeIdeal` instances.
  
      References
      ==========
  
      .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.*
         (See Algorithm 6.2.9.)
  
      """
      if T is None and ZK is None:
          raise ValueError('At least one of T or ZK must be provided.')
      if ZK is not None:
          _check_formal_conditions_for_maximal_order(ZK)
      if T is None:
          T = ZK.parent.T
      radicals = {}
      if dK is None or ZK is None:
          ZK, dK = round_two(T, radicals=radicals)
      dT = T.discriminant()
      f_squared = dT // dK
      if f_squared % p != 0:
          return _prime_decomp_easy_case(p, ZK)
      radical = radical or radicals.get(p) or nilradical_mod_p(ZK, p)
      stack = [radical]
      primes = []
      while stack:
          I = stack.pop()
          N, G = _prime_decomp_compute_kernel(I, p, ZK)
          if N.n == 1:
              P = _prime_decomp_maximal_ideal(I, p, ZK)
              primes.append(P)
          else:
              I1, I2 = _prime_decomp_split_ideal(I, p, N, G, ZK)
              stack.extend([I1, I2])
      return primes
  