"""Prime ideals in number fields. """
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
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)
[docs]
class PrimeIdeal(IntegerPowerable):
r"""
A prime ideal in a ring of algebraic integers.
"""
[docs]
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 __str__(self):
if self.is_inert:
return f'({self.p})'
return f'({self.p}, {self.alpha.as_expr()})'
@property
def is_inert(self):
"""
Say whether the rational prime we divide is inert, i.e. stays prime in
our ring of integers.
"""
return self.f == self.ZK.n
[docs]
def repr(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].repr())
[ (11, x**3 + 5*x**2 + 4*x - 1) e=1, f=3 ]
>>> print(P[0].repr(field_gen=zeta))
[ (11, zeta**3 + 5*zeta**2 + 4*zeta - 1) e=1, f=3 ]
>>> print(P[0].repr(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.repr()
[docs]
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
[docs]
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__
[docs]
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
[docs]
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
[docs]
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)
[docs]
def reduce_element(self, elt):
"""
Reduce a :py:class:`~.PowerBasisElement` to a "small representative"
modulo this prime ideal.
Parameters
==========
elt : :py:class:`~.PowerBasisElement`
The element to be reduced.
Returns
=======
:py:class:`~.PowerBasisElement`
The reduced element.
See Also
========
reduce_ANP
reduce_alg_num
.Submodule.reduce_element
"""
return self.as_submodule().reduce_element(elt)
[docs]
def reduce_ANP(self, a):
"""
Reduce an :py:class:`~.ANP` to a "small representative" modulo this
prime ideal.
Parameters
==========
elt : :py:class:`~.ANP`
The element to be reduced.
Returns
=======
:py:class:`~.ANP`
The reduced element.
See Also
========
reduce_element
reduce_alg_num
.Submodule.reduce_element
"""
elt = self.ZK.parent.element_from_ANP(a)
red = self.reduce_element(elt)
return red.to_ANP()
[docs]
def reduce_alg_num(self, a):
"""
Reduce an :py:class:`~.AlgebraicNumber` to a "small representative"
modulo this prime ideal.
Parameters
==========
elt : :py:class:`~.AlgebraicNumber`
The element to be reduced.
Returns
=======
:py:class:`~.AlgebraicNumber`
The reduced element.
See Also
========
reduce_element
reduce_ANP
.Submodule.reduce_element
"""
elt = self.ZK.parent.element_from_alg_num(a)
red = self.reduce_element(elt)
return a.field_element(list(reversed(red.QQ_col.flat())))
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.convert_to(ZZ), denom=ZK.denom)
return beta
[docs]
@public
def prime_valuation(I, P):
r"""
Compute the *P*-adic valuation for an integral ideal *I*.
Examples
========
>>> from sympy import QQ
>>> from sympy.polys.numberfields import prime_valuation
>>> K = QQ.cyclotomic_field(5)
>>> 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()
if len(fl) == 1 and fl[0][1] == 1:
return [PrimeIdeal(ZK, p, ZK.parent.zero(), ZK.n, 1)]
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 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.to_list()))
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
[docs]
@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