Coverage for /usr/lib/python3/dist-packages/sympy/polys/numberfields/primes.py: 19%
225 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""Prime ideals in number fields. """
3from sympy.polys.polytools import Poly
4from sympy.polys.domains.finitefield import FF
5from sympy.polys.domains.rationalfield import QQ
6from sympy.polys.domains.integerring import ZZ
7from sympy.polys.matrices.domainmatrix import DomainMatrix
8from sympy.polys.polyerrors import CoercionFailed
9from sympy.polys.polyutils import IntegerPowerable
10from sympy.utilities.decorator import public
11from .basis import round_two, nilradical_mod_p
12from .exceptions import StructureError
13from .modules import ModuleEndomorphism, find_min_poly
14from .utilities import coeff_search, supplement_a_subspace
17def _check_formal_conditions_for_maximal_order(submodule):
18 r"""
19 Several functions in this module accept an argument which is to be a
20 :py:class:`~.Submodule` representing the maximal order in a number field,
21 such as returned by the :py:func:`~sympy.polys.numberfields.basis.round_two`
22 algorithm.
24 We do not attempt to check that the given ``Submodule`` actually represents
25 a maximal order, but we do check a basic set of formal conditions that the
26 ``Submodule`` must satisfy, at a minimum. The purpose is to catch an
27 obviously ill-formed argument.
28 """
29 prefix = 'The submodule representing the maximal order should '
30 cond = None
31 if not submodule.is_power_basis_submodule():
32 cond = 'be a direct submodule of a power basis.'
33 elif not submodule.starts_with_unity():
34 cond = 'have 1 as its first generator.'
35 elif not submodule.is_sq_maxrank_HNF():
36 cond = 'have square matrix, of maximal rank, in Hermite Normal Form.'
37 if cond is not None:
38 raise StructureError(prefix + cond)
41class PrimeIdeal(IntegerPowerable):
42 r"""
43 A prime ideal in a ring of algebraic integers.
44 """
46 def __init__(self, ZK, p, alpha, f, e=None):
47 """
48 Parameters
49 ==========
51 ZK : :py:class:`~.Submodule`
52 The maximal order where this ideal lives.
53 p : int
54 The rational prime this ideal divides.
55 alpha : :py:class:`~.PowerBasisElement`
56 Such that the ideal is equal to ``p*ZK + alpha*ZK``.
57 f : int
58 The inertia degree.
59 e : int, ``None``, optional
60 The ramification index, if already known. If ``None``, we will
61 compute it here.
63 """
64 _check_formal_conditions_for_maximal_order(ZK)
65 self.ZK = ZK
66 self.p = p
67 self.alpha = alpha
68 self.f = f
69 self._test_factor = None
70 self.e = e if e is not None else self.valuation(p * ZK)
72 def __str__(self):
73 if self.is_inert:
74 return f'({self.p})'
75 return f'({self.p}, {self.alpha.as_expr()})'
77 @property
78 def is_inert(self):
79 """
80 Say whether the rational prime we divide is inert, i.e. stays prime in
81 our ring of integers.
82 """
83 return self.f == self.ZK.n
85 def repr(self, field_gen=None, just_gens=False):
86 """
87 Print a representation of this prime ideal.
89 Examples
90 ========
92 >>> from sympy import cyclotomic_poly, QQ
93 >>> from sympy.abc import x, zeta
94 >>> T = cyclotomic_poly(7, x)
95 >>> K = QQ.algebraic_field((T, zeta))
96 >>> P = K.primes_above(11)
97 >>> print(P[0].repr())
98 [ (11, x**3 + 5*x**2 + 4*x - 1) e=1, f=3 ]
99 >>> print(P[0].repr(field_gen=zeta))
100 [ (11, zeta**3 + 5*zeta**2 + 4*zeta - 1) e=1, f=3 ]
101 >>> print(P[0].repr(field_gen=zeta, just_gens=True))
102 (11, zeta**3 + 5*zeta**2 + 4*zeta - 1)
104 Parameters
105 ==========
107 field_gen : :py:class:`~.Symbol`, ``None``, optional (default=None)
108 The symbol to use for the generator of the field. This will appear
109 in our representation of ``self.alpha``. If ``None``, we use the
110 variable of the defining polynomial of ``self.ZK``.
111 just_gens : bool, optional (default=False)
112 If ``True``, just print the "(p, alpha)" part, showing "just the
113 generators" of the prime ideal. Otherwise, print a string of the
114 form "[ (p, alpha) e=..., f=... ]", giving the ramification index
115 and inertia degree, along with the generators.
117 """
118 field_gen = field_gen or self.ZK.parent.T.gen
119 p, alpha, e, f = self.p, self.alpha, self.e, self.f
120 alpha_rep = str(alpha.numerator(x=field_gen).as_expr())
121 if alpha.denom > 1:
122 alpha_rep = f'({alpha_rep})/{alpha.denom}'
123 gens = f'({p}, {alpha_rep})'
124 if just_gens:
125 return gens
126 return f'[ {gens} e={e}, f={f} ]'
128 def __repr__(self):
129 return self.repr()
131 def as_submodule(self):
132 r"""
133 Represent this prime ideal as a :py:class:`~.Submodule`.
135 Explanation
136 ===========
138 The :py:class:`~.PrimeIdeal` class serves to bundle information about
139 a prime ideal, such as its inertia degree, ramification index, and
140 two-generator representation, as well as to offer helpful methods like
141 :py:meth:`~.PrimeIdeal.valuation` and
142 :py:meth:`~.PrimeIdeal.test_factor`.
144 However, in order to be added and multiplied by other ideals or
145 rational numbers, it must first be converted into a
146 :py:class:`~.Submodule`, which is a class that supports these
147 operations.
149 In many cases, the user need not perform this conversion deliberately,
150 since it is automatically performed by the arithmetic operator methods
151 :py:meth:`~.PrimeIdeal.__add__` and :py:meth:`~.PrimeIdeal.__mul__`.
153 Raising a :py:class:`~.PrimeIdeal` to a non-negative integer power is
154 also supported.
156 Examples
157 ========
159 >>> from sympy import Poly, cyclotomic_poly, prime_decomp
160 >>> T = Poly(cyclotomic_poly(7))
161 >>> P0 = prime_decomp(7, T)[0]
162 >>> print(P0**6 == 7*P0.ZK)
163 True
165 Note that, on both sides of the equation above, we had a
166 :py:class:`~.Submodule`. In the next equation we recall that adding
167 ideals yields their GCD. This time, we need a deliberate conversion
168 to :py:class:`~.Submodule` on the right:
170 >>> print(P0 + 7*P0.ZK == P0.as_submodule())
171 True
173 Returns
174 =======
176 :py:class:`~.Submodule`
177 Will be equal to ``self.p * self.ZK + self.alpha * self.ZK``.
179 See Also
180 ========
182 __add__
183 __mul__
185 """
186 M = self.p * self.ZK + self.alpha * self.ZK
187 # Pre-set expensive boolean properties whose value we already know:
188 M._starts_with_unity = False
189 M._is_sq_maxrank_HNF = True
190 return M
192 def __eq__(self, other):
193 if isinstance(other, PrimeIdeal):
194 return self.as_submodule() == other.as_submodule()
195 return NotImplemented
197 def __add__(self, other):
198 """
199 Convert to a :py:class:`~.Submodule` and add to another
200 :py:class:`~.Submodule`.
202 See Also
203 ========
205 as_submodule
207 """
208 return self.as_submodule() + other
210 __radd__ = __add__
212 def __mul__(self, other):
213 """
214 Convert to a :py:class:`~.Submodule` and multiply by another
215 :py:class:`~.Submodule` or a rational number.
217 See Also
218 ========
220 as_submodule
222 """
223 return self.as_submodule() * other
225 __rmul__ = __mul__
227 def _zeroth_power(self):
228 return self.ZK
230 def _first_power(self):
231 return self
233 def test_factor(self):
234 r"""
235 Compute a test factor for this prime ideal.
237 Explanation
238 ===========
240 Write $\mathfrak{p}$ for this prime ideal, $p$ for the rational prime
241 it divides. Then, for computing $\mathfrak{p}$-adic valuations it is
242 useful to have a number $\beta \in \mathbb{Z}_K$ such that
243 $p/\mathfrak{p} = p \mathbb{Z}_K + \beta \mathbb{Z}_K$.
245 Essentially, this is the same as the number $\Psi$ (or the "reagent")
246 from Kummer's 1847 paper (*Ueber die Zerlegung...*, Crelle vol. 35) in
247 which ideal divisors were invented.
248 """
249 if self._test_factor is None:
250 self._test_factor = _compute_test_factor(self.p, [self.alpha], self.ZK)
251 return self._test_factor
253 def valuation(self, I):
254 r"""
255 Compute the $\mathfrak{p}$-adic valuation of integral ideal I at this
256 prime ideal.
258 Parameters
259 ==========
261 I : :py:class:`~.Submodule`
263 See Also
264 ========
266 prime_valuation
268 """
269 return prime_valuation(I, self)
271 def reduce_element(self, elt):
272 """
273 Reduce a :py:class:`~.PowerBasisElement` to a "small representative"
274 modulo this prime ideal.
276 Parameters
277 ==========
279 elt : :py:class:`~.PowerBasisElement`
280 The element to be reduced.
282 Returns
283 =======
285 :py:class:`~.PowerBasisElement`
286 The reduced element.
288 See Also
289 ========
291 reduce_ANP
292 reduce_alg_num
293 .Submodule.reduce_element
295 """
296 return self.as_submodule().reduce_element(elt)
298 def reduce_ANP(self, a):
299 """
300 Reduce an :py:class:`~.ANP` to a "small representative" modulo this
301 prime ideal.
303 Parameters
304 ==========
306 elt : :py:class:`~.ANP`
307 The element to be reduced.
309 Returns
310 =======
312 :py:class:`~.ANP`
313 The reduced element.
315 See Also
316 ========
318 reduce_element
319 reduce_alg_num
320 .Submodule.reduce_element
322 """
323 elt = self.ZK.parent.element_from_ANP(a)
324 red = self.reduce_element(elt)
325 return red.to_ANP()
327 def reduce_alg_num(self, a):
328 """
329 Reduce an :py:class:`~.AlgebraicNumber` to a "small representative"
330 modulo this prime ideal.
332 Parameters
333 ==========
335 elt : :py:class:`~.AlgebraicNumber`
336 The element to be reduced.
338 Returns
339 =======
341 :py:class:`~.AlgebraicNumber`
342 The reduced element.
344 See Also
345 ========
347 reduce_element
348 reduce_ANP
349 .Submodule.reduce_element
351 """
352 elt = self.ZK.parent.element_from_alg_num(a)
353 red = self.reduce_element(elt)
354 return a.field_element(list(reversed(red.QQ_col.flat())))
357def _compute_test_factor(p, gens, ZK):
358 r"""
359 Compute the test factor for a :py:class:`~.PrimeIdeal` $\mathfrak{p}$.
361 Parameters
362 ==========
364 p : int
365 The rational prime $\mathfrak{p}$ divides
367 gens : list of :py:class:`PowerBasisElement`
368 A complete set of generators for $\mathfrak{p}$ over *ZK*, EXCEPT that
369 an element equivalent to rational *p* can and should be omitted (since
370 it has no effect except to waste time).
372 ZK : :py:class:`~.Submodule`
373 The maximal order where the prime ideal $\mathfrak{p}$ lives.
375 Returns
376 =======
378 :py:class:`~.PowerBasisElement`
380 References
381 ==========
383 .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.*
384 (See Proposition 4.8.15.)
386 """
387 _check_formal_conditions_for_maximal_order(ZK)
388 E = ZK.endomorphism_ring()
389 matrices = [E.inner_endomorphism(g).matrix(modulus=p) for g in gens]
390 B = DomainMatrix.zeros((0, ZK.n), FF(p)).vstack(*matrices)
391 # A nonzero element of the nullspace of B will represent a
392 # lin comb over the omegas which (i) is not a multiple of p
393 # (since it is nonzero over FF(p)), while (ii) is such that
394 # its product with each g in gens _is_ a multiple of p (since
395 # B represents multiplication by these generators). Theory
396 # predicts that such an element must exist, so nullspace should
397 # be non-trivial.
398 x = B.nullspace()[0, :].transpose()
399 beta = ZK.parent(ZK.matrix * x, denom=ZK.denom)
400 return beta
403@public
404def prime_valuation(I, P):
405 r"""
406 Compute the *P*-adic valuation for an integral ideal *I*.
408 Examples
409 ========
411 >>> from sympy import QQ
412 >>> from sympy.polys.numberfields import prime_valuation
413 >>> K = QQ.cyclotomic_field(5)
414 >>> P = K.primes_above(5)
415 >>> ZK = K.maximal_order()
416 >>> print(prime_valuation(25*ZK, P[0]))
417 8
419 Parameters
420 ==========
422 I : :py:class:`~.Submodule`
423 An integral ideal whose valuation is desired.
425 P : :py:class:`~.PrimeIdeal`
426 The prime at which to compute the valuation.
428 Returns
429 =======
431 int
433 See Also
434 ========
436 .PrimeIdeal.valuation
438 References
439 ==========
441 .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.*
442 (See Algorithm 4.8.17.)
444 """
445 p, ZK = P.p, P.ZK
446 n, W, d = ZK.n, ZK.matrix, ZK.denom
448 A = W.convert_to(QQ).inv() * I.matrix * d / I.denom
449 # Although A must have integer entries, given that I is an integral ideal,
450 # as a DomainMatrix it will still be over QQ, so we convert back:
451 A = A.convert_to(ZZ)
452 D = A.det()
453 if D % p != 0:
454 return 0
456 beta = P.test_factor()
458 f = d ** n // W.det()
459 need_complete_test = (f % p == 0)
460 v = 0
461 while True:
462 # Entering the loop, the cols of A represent lin combs of omegas.
463 # Turn them into lin combs of thetas:
464 A = W * A
465 # And then one column at a time...
466 for j in range(n):
467 c = ZK.parent(A[:, j], denom=d)
468 c *= beta
469 # ...turn back into lin combs of omegas, after multiplying by beta:
470 c = ZK.represent(c).flat()
471 for i in range(n):
472 A[i, j] = c[i]
473 if A[n - 1, n - 1].element % p != 0:
474 break
475 A = A / p
476 # As noted above, domain converts to QQ even when division goes evenly.
477 # So must convert back, even when we don't "need_complete_test".
478 if need_complete_test:
479 # In this case, having a non-integer entry is actually just our
480 # halting condition.
481 try:
482 A = A.convert_to(ZZ)
483 except CoercionFailed:
484 break
485 else:
486 # In this case theory says we should not have any non-integer entries.
487 A = A.convert_to(ZZ)
488 v += 1
489 return v
492def _two_elt_rep(gens, ZK, p, f=None, Np=None):
493 r"""
494 Given a set of *ZK*-generators of a prime ideal, compute a set of just two
495 *ZK*-generators for the same ideal, one of which is *p* itself.
497 Parameters
498 ==========
500 gens : list of :py:class:`PowerBasisElement`
501 Generators for the prime ideal over *ZK*, the ring of integers of the
502 field $K$.
504 ZK : :py:class:`~.Submodule`
505 The maximal order in $K$.
507 p : int
508 The rational prime divided by the prime ideal.
510 f : int, optional
511 The inertia degree of the prime ideal, if known.
513 Np : int, optional
514 The norm $p^f$ of the prime ideal, if known.
515 NOTE: There is no reason to supply both *f* and *Np*. Either one will
516 save us from having to compute the norm *Np* ourselves. If both are known,
517 *Np* is preferred since it saves one exponentiation.
519 Returns
520 =======
522 :py:class:`~.PowerBasisElement` representing a single algebraic integer
523 alpha such that the prime ideal is equal to ``p*ZK + alpha*ZK``.
525 References
526 ==========
528 .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.*
529 (See Algorithm 4.7.10.)
531 """
532 _check_formal_conditions_for_maximal_order(ZK)
533 pb = ZK.parent
534 T = pb.T
535 # Detect the special cases in which either (a) all generators are multiples
536 # of p, or (b) there are no generators (so `all` is vacuously true):
537 if all((g % p).equiv(0) for g in gens):
538 return pb.zero()
540 if Np is None:
541 if f is not None:
542 Np = p**f
543 else:
544 Np = abs(pb.submodule_from_gens(gens).matrix.det())
546 omega = ZK.basis_element_pullbacks()
547 beta = [p*om for om in omega[1:]] # note: we omit omega[0] == 1
548 beta += gens
549 search = coeff_search(len(beta), 1)
550 for c in search:
551 alpha = sum(ci*betai for ci, betai in zip(c, beta))
552 # Note: It may be tempting to reduce alpha mod p here, to try to work
553 # with smaller numbers, but must not do that, as it can result in an
554 # infinite loop! E.g. try factoring 2 in Q(sqrt(-7)).
555 n = alpha.norm(T) // Np
556 if n % p != 0:
557 # Now can reduce alpha mod p.
558 return alpha % p
561def _prime_decomp_easy_case(p, ZK):
562 r"""
563 Compute the decomposition of rational prime *p* in the ring of integers
564 *ZK* (given as a :py:class:`~.Submodule`), in the "easy case", i.e. the
565 case where *p* does not divide the index of $\theta$ in *ZK*, where
566 $\theta$ is the generator of the ``PowerBasis`` of which *ZK* is a
567 ``Submodule``.
568 """
569 T = ZK.parent.T
570 T_bar = Poly(T, modulus=p)
571 lc, fl = T_bar.factor_list()
572 if len(fl) == 1 and fl[0][1] == 1:
573 return [PrimeIdeal(ZK, p, ZK.parent.zero(), ZK.n, 1)]
574 return [PrimeIdeal(ZK, p,
575 ZK.parent.element_from_poly(Poly(t, domain=ZZ)),
576 t.degree(), e)
577 for t, e in fl]
580def _prime_decomp_compute_kernel(I, p, ZK):
581 r"""
582 Parameters
583 ==========
585 I : :py:class:`~.Module`
586 An ideal of ``ZK/pZK``.
587 p : int
588 The rational prime being factored.
589 ZK : :py:class:`~.Submodule`
590 The maximal order.
592 Returns
593 =======
595 Pair ``(N, G)``, where:
597 ``N`` is a :py:class:`~.Module` representing the kernel of the map
598 ``a |--> a**p - a`` on ``(O/pO)/I``, guaranteed to be a module with
599 unity.
601 ``G`` is a :py:class:`~.Module` representing a basis for the separable
602 algebra ``A = O/I`` (see Cohen).
604 """
605 W = I.matrix
606 n, r = W.shape
607 # Want to take the Fp-basis given by the columns of I, adjoin (1, 0, ..., 0)
608 # (which we know is not already in there since I is a basis for a prime ideal)
609 # and then supplement this with additional columns to make an invertible n x n
610 # matrix. This will then represent a full basis for ZK, whose first r columns
611 # are pullbacks of the basis for I.
612 if r == 0:
613 B = W.eye(n, ZZ)
614 else:
615 B = W.hstack(W.eye(n, ZZ)[:, 0])
616 if B.shape[1] < n:
617 B = supplement_a_subspace(B.convert_to(FF(p))).convert_to(ZZ)
619 G = ZK.submodule_from_matrix(B)
620 # Must compute G's multiplication table _before_ discarding the first r
621 # columns. (See Step 9 in Alg 6.2.9 in Cohen, where the betas are actually
622 # needed in order to represent each product of gammas. However, once we've
623 # found the representations, then we can ignore the betas.)
624 G.compute_mult_tab()
625 G = G.discard_before(r)
627 phi = ModuleEndomorphism(G, lambda x: x**p - x)
628 N = phi.kernel(modulus=p)
629 assert N.starts_with_unity()
630 return N, G
633def _prime_decomp_maximal_ideal(I, p, ZK):
634 r"""
635 We have reached the case where we have a maximal (hence prime) ideal *I*,
636 which we know because the quotient ``O/I`` is a field.
638 Parameters
639 ==========
641 I : :py:class:`~.Module`
642 An ideal of ``O/pO``.
643 p : int
644 The rational prime being factored.
645 ZK : :py:class:`~.Submodule`
646 The maximal order.
648 Returns
649 =======
651 :py:class:`~.PrimeIdeal` instance representing this prime
653 """
654 m, n = I.matrix.shape
655 f = m - n
656 G = ZK.matrix * I.matrix
657 gens = [ZK.parent(G[:, j], denom=ZK.denom) for j in range(G.shape[1])]
658 alpha = _two_elt_rep(gens, ZK, p, f=f)
659 return PrimeIdeal(ZK, p, alpha, f)
662def _prime_decomp_split_ideal(I, p, N, G, ZK):
663 r"""
664 Perform the step in the prime decomposition algorithm where we have determined
665 the the quotient ``ZK/I`` is _not_ a field, and we want to perform a non-trivial
666 factorization of *I* by locating an idempotent element of ``ZK/I``.
667 """
668 assert I.parent == ZK and G.parent is ZK and N.parent is G
669 # Since ZK/I is not a field, the kernel computed in the previous step contains
670 # more than just the prime field Fp, and our basis N for the nullspace therefore
671 # contains at least a second column (which represents an element outside Fp).
672 # Let alpha be such an element:
673 alpha = N(1).to_parent()
674 assert alpha.module is G
676 alpha_powers = []
677 m = find_min_poly(alpha, FF(p), powers=alpha_powers)
678 # TODO (future work):
679 # We don't actually need full factorization, so might use a faster method
680 # to just break off a single non-constant factor m1?
681 lc, fl = m.factor_list()
682 m1 = fl[0][0]
683 m2 = m.quo(m1)
684 U, V, g = m1.gcdex(m2)
685 # Sanity check: theory says m is squarefree, so m1, m2 should be coprime:
686 assert g == 1
687 E = list(reversed(Poly(U * m1, domain=ZZ).rep.rep))
688 eps1 = sum(E[i]*alpha_powers[i] for i in range(len(E)))
689 eps2 = 1 - eps1
690 idemps = [eps1, eps2]
691 factors = []
692 for eps in idemps:
693 e = eps.to_parent()
694 assert e.module is ZK
695 D = I.matrix.convert_to(FF(p)).hstack(*[
696 (e * om).column(domain=FF(p)) for om in ZK.basis_elements()
697 ])
698 W = D.columnspace().convert_to(ZZ)
699 H = ZK.submodule_from_matrix(W)
700 factors.append(H)
701 return factors
704@public
705def prime_decomp(p, T=None, ZK=None, dK=None, radical=None):
706 r"""
707 Compute the decomposition of rational prime *p* in a number field.
709 Explanation
710 ===========
712 Ordinarily this should be accessed through the
713 :py:meth:`~.AlgebraicField.primes_above` method of an
714 :py:class:`~.AlgebraicField`.
716 Examples
717 ========
719 >>> from sympy import Poly, QQ
720 >>> from sympy.abc import x, theta
721 >>> T = Poly(x ** 3 + x ** 2 - 2 * x + 8)
722 >>> K = QQ.algebraic_field((T, theta))
723 >>> print(K.primes_above(2))
724 [[ (2, x**2 + 1) e=1, f=1 ], [ (2, (x**2 + 3*x + 2)/2) e=1, f=1 ],
725 [ (2, (3*x**2 + 3*x)/2) e=1, f=1 ]]
727 Parameters
728 ==========
730 p : int
731 The rational prime whose decomposition is desired.
733 T : :py:class:`~.Poly`, optional
734 Monic irreducible polynomial defining the number field $K$ in which to
735 factor. NOTE: at least one of *T* or *ZK* must be provided.
737 ZK : :py:class:`~.Submodule`, optional
738 The maximal order for $K$, if already known.
739 NOTE: at least one of *T* or *ZK* must be provided.
741 dK : int, optional
742 The discriminant of the field $K$, if already known.
744 radical : :py:class:`~.Submodule`, optional
745 The nilradical mod *p* in the integers of $K$, if already known.
747 Returns
748 =======
750 List of :py:class:`~.PrimeIdeal` instances.
752 References
753 ==========
755 .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.*
756 (See Algorithm 6.2.9.)
758 """
759 if T is None and ZK is None:
760 raise ValueError('At least one of T or ZK must be provided.')
761 if ZK is not None:
762 _check_formal_conditions_for_maximal_order(ZK)
763 if T is None:
764 T = ZK.parent.T
765 radicals = {}
766 if dK is None or ZK is None:
767 ZK, dK = round_two(T, radicals=radicals)
768 dT = T.discriminant()
769 f_squared = dT // dK
770 if f_squared % p != 0:
771 return _prime_decomp_easy_case(p, ZK)
772 radical = radical or radicals.get(p) or nilradical_mod_p(ZK, p)
773 stack = [radical]
774 primes = []
775 while stack:
776 I = stack.pop()
777 N, G = _prime_decomp_compute_kernel(I, p, ZK)
778 if N.n == 1:
779 P = _prime_decomp_maximal_ideal(I, p, ZK)
780 primes.append(P)
781 else:
782 I1, I2 = _prime_decomp_split_ideal(I, p, N, G, ZK)
783 stack.extend([I1, I2])
784 return primes