Repository: sympy/sympy
Pull request: Improve printing for `PrimeIdeal`
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->

Fixes #22974 


#### Brief description of what is fixed or changed

* It is now possible for a `Submodule` representing the maximal ideal of an `AlgebraicField` to have a reference back to that field.
* When the maximal ideal is constructed via `AlgebraicField.maximal_order()` or `AlgebraicField.integral_basis()` it _does_ have such a reference.
* As a consequence, we're able to improve the printing methods for `PrimeIdeal`, so that they use the alias of the associated field's primitive element, if any. In particular, we:
  - Support latex printing
  - Rename `_pretty()` --> `repr()` since this is not 2D printing.
  - Provide a `__str__()` method, which prints less info than the `__repr__()` method.


#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* polys
  * Improve printing for `PrimeIdeal`
<!-- END RELEASE NOTES -->


Problem statement:
polys/numberfields: `PrimeIdeal` printing methods need work.
Two methods in the `PrimeIdeal` class, `pretty()` (for printing) and `reduce_poly()`, involve generating polynomial expressions. Ideally, the generator symbol in these expressions should be able to be automatically the same as the `alias` symbol for the generator of the surrounding `AlgebraicField`.

To do this, instances of `PrimeIdeal` will need a reference to that field, which they currently lack.

Hints from issue discussion:
So, in #22977 these two methods were made private, so that they could still be changed after v1.10.
Here are some thoughts on improving them for v1.11:

## `PrimeIdeal.pretty`

First, this is misnamed, because it's not "pretty printing" in the sense of 2D displays. It's really just a `repr`-type string.

Apart from this, we should have better `__str__` and LaTeX printing for the `PrimeIdeal` class, along the lines discussed in the OP (i.e. where the alias of the primitive element of the associated `AlgebraicField` (if any) can be automatically used).

Before that can happen though, it would be best if #23338 could be resolved.


## `PrimeIdeal.reduce_poly`

I'd like to improve this, and keep it around because it is useful, but I'd like to keep it private (leading underscore) in case we want to change the argument or return types as the need becomes more clear.

The basic operation is an important one: we're mapping elements of a ring of integers `R` into the finite field `R/P`, where `P` is a prime ideal. So it's like the [`nfmodpr`](http://pari.math.u-bordeaux.fr/dochtml/html-stable/General_number_fields.html#nfmodpr) function of PARI.

The existing implementation is funny though, in that it accepts a `Poly`. It makes more sense to accept a class like `ANP` or `AlgebraicNumber` (or both), which represents an algebraic number more directly. (And the method should be called something simpler like `_reduce()`).

As for the return value, we're using a `Poly` mod `p`, which does make sense as a representation of an element of a finite field of order `p^n`, but again it's not clear if this should change in the future. We might just want a `dup` instead of a `Poly`, or at some point SymPy might provide `GF(p^n)` as a domain type (for `n > 1`), and then we might want it to be one of those (in particular if `GF(p^n)` lets you choose the irreducible poly by which it's defined).

Commit message: Improve printing for `PrimeIdeal`

* Support latex printing
* Rename `_pretty()` --> `repr()` since this is not 2D printing.
* Provide a `__str__()` method, which prints less info than the `__repr__()` method.

Patch to review:
diff --git a/sympy/polys/domains/algebraicfield.py b/sympy/polys/domains/algebraicfield.py
index 830d820ba24e..241ce819a2a8 100644
--- a/sympy/polys/domains/algebraicfield.py
+++ b/sympy/polys/domains/algebraicfield.py
@@ -213,8 +213,7 @@ class AlgebraicField(Field, CharacteristicZero, SimpleDomain):
     >>> K
     QQ<exp(2*I*pi/7)>
     >>> K.primes_above(11)
-    [[ (11, _x**3 + 5*_x**2 + 4*_x - 1) e=1, f=3 ],
-     [ (11, _x**3 - 4*_x**2 - 5*_x - 1) e=1, f=3 ]]
+    [(11, _x**3 + 5*_x**2 + 4*_x - 1), (11, _x**3 - 4*_x**2 - 5*_x - 1)]
 
     Notes
     =====
@@ -433,7 +432,7 @@ def from_GaussianRationalField(K1, a, K0):
 
     def _do_round_two(self):
         from sympy.polys.numberfields.basis import round_two
-        ZK, dK = round_two(self.ext.minpoly, radicals=self._nilradicals_mod_p)
+        ZK, dK = round_two(self, radicals=self._nilradicals_mod_p)
         self._maximal_order = ZK
         self._discriminant = dK
 
@@ -512,7 +511,6 @@ def integral_basis(self, fmt=None):
             return [self.to_alg_num(b) for b in B]
         return B
 
-
     def discriminant(self):
         """Get the discriminant of the field."""
         if self._discriminant is None:
diff --git a/sympy/polys/numberfields/basis.py b/sympy/polys/numberfields/basis.py
index 9cb4b0e0aab4..a4ff2df0b286 100644
--- a/sympy/polys/numberfields/basis.py
+++ b/sympy/polys/numberfields/basis.py
@@ -1,6 +1,7 @@
 """Computing integral bases for number fields. """
 
 from sympy.polys.polytools import Poly
+from sympy.polys.domains.algebraicfield import AlgebraicField
 from sympy.polys.domains.integerring import ZZ
 from sympy.polys.domains.rationalfield import QQ
 from sympy.polys.polyerrors import CoercionFailed
@@ -103,6 +104,10 @@ def round_two(T, radicals=None):
     polynomial *T* over :ref:`ZZ`. This computes an integral basis and the
     discriminant for the field $K = \mathbb{Q}[x]/(T(x))$.
 
+    Alternatively, you may pass an :py:class:`~.AlgebraicField` instance, in
+    place of the polynomial *T*, in which case the algorithm is applied to the
+    minimal polynomial for the field's primitive element.
+
     Ordinarily this function need not be called directly, as one can instead
     access the :py:meth:`~.AlgebraicField.maximal_order`,
     :py:meth:`~.AlgebraicField.integral_basis`, and
@@ -147,9 +152,10 @@ def round_two(T, radicals=None):
     Parameters
     ==========
 
-    T : :py:class:`~.Poly`
-        The irreducible monic polynomial over :ref:`ZZ` defining the number
-        field.
+    T : :py:class:`~.Poly`, :py:class:`~.AlgebraicField`
+        Either (1) the irreducible monic polynomial over :ref:`ZZ` defining the
+        number field, or (2) an :py:class:`~.AlgebraicField` representing the
+        number field itself.
 
     radicals : dict, optional
         This is a way for any $p$-radicals (if computed) to be returned by
@@ -182,6 +188,9 @@ def round_two(T, radicals=None):
     .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.*
 
     """
+    K = None
+    if isinstance(T, AlgebraicField):
+        K, T = T, T.ext.minpoly_of_element()
     if T.domain == QQ:
         try:
             T = Poly(T, domain=ZZ)
@@ -198,7 +207,7 @@ def round_two(T, radicals=None):
     # D must be 0 or 1 mod 4 (see Cohen Sec 4.4), which ensures we can write
     # it in the form D = D_0 * F**2, where D_0 is 1 or a fundamental discriminant.
     _, F = extract_fundamental_discriminant(D)
-    Ztheta = PowerBasis(T)
+    Ztheta = PowerBasis(K or T)
     H = Ztheta.whole_submodule()
     nilrad = None
     while F:
diff --git a/sympy/polys/numberfields/modules.py b/sympy/polys/numberfields/modules.py
index 22813c8065f3..bd852553d56c 100644
--- a/sympy/polys/numberfields/modules.py
+++ b/sympy/polys/numberfields/modules.py
@@ -10,8 +10,9 @@
 
 * For a :py:class:`~.PowerBasis`, the generators are the first $n$ powers
   (starting with the zeroth) of an algebraic integer $\theta$ of degree $n$.
-  The :py:class:`~.PowerBasis` is constructed by passing the minimal
-  polynomial of $\theta$.
+  The :py:class:`~.PowerBasis` is constructed by passing either the minimal
+  polynomial of $\theta$, or an :py:class:`~.AlgebraicField` having $\theta$
+  as its primitive element.
 
 * For a :py:class:`~.Submodule`, the generators are a set of
   $\mathbb{Q}$-linear combinations of the generators of another module. That
@@ -181,6 +182,7 @@
 from sympy.core.symbol import Dummy
 from sympy.polys.polytools import Poly
 from sympy.polys.densetools import dup_clear_denoms
+from sympy.polys.domains.algebraicfield import AlgebraicField
 from sympy.polys.domains.finitefield import FF
 from sympy.polys.domains.rationalfield import QQ
 from sympy.polys.domains.integerring import ZZ
@@ -436,6 +438,28 @@ def nearest_common_ancestor(self, other):
                 break
         return nca
 
+    @property
+    def number_field(self):
+        r"""
+        Return the associated :py:class:`~.AlgebraicField`, if any.
+
+        Explanation
+        ===========
+
+        A :py:class:`~.PowerBasis` can be constructed on a :py:class:`~.Poly`
+        $f$ or on an :py:class:`~.AlgebraicField` $K$. In the latter case, the
+        :py:class:`~.PowerBasis` and all its descendant modules will return $K$
+        as their ``.number_field`` property, while in the former case they will
+        all return ``None``.
+
+        Returns
+        =======
+
+        :py:class:`~.AlgebraicField`, ``None``
+
+        """
+        return self.power_basis_ancestor().number_field
+
     def is_compat_col(self, col):
         """Say whether *col* is a suitable column vector for this module."""
         return isinstance(col, DomainMatrix) and col.shape == (self.n, 1) and col.domain.is_ZZ
@@ -678,15 +702,30 @@ def __init__(self, T):
         Parameters
         ==========
 
-        T : :py:class:`~.Poly`
-            The monic, irreducible, univariate polynomial over :ref:`ZZ`, a
-            root of which is the generator of the power basis.
+        T : :py:class:`~.Poly`, :py:class:`~.AlgebraicField`
+            Either (1) the monic, irreducible, univariate polynomial over
+            :ref:`ZZ`, a root of which is the generator of the power basis,
+            or (2) an :py:class:`~.AlgebraicField` whose primitive element
+            is the generator of the power basis.
 
         """
+        K = None
+        if isinstance(T, AlgebraicField):
+            K, T = T, T.ext.minpoly_of_element()
+        if T.domain == QQ:
+            try:
+                T = Poly(T, domain=ZZ)
+            except CoercionFailed:
+                raise ValueError('A polynomial over ZZ is required')
+        self.K = K
         self.T = T
         self._n = T.degree()
         self._mult_tab = None
 
+    @property
+    def number_field(self):
+        return self.K
+
     def __repr__(self):
         return f'PowerBasis({self.T.as_expr()})'
 
@@ -1542,6 +1581,28 @@ def poly(self, x=None):
         """Obtain the number as a polynomial over :ref:`QQ`."""
         return self.numerator(x=x) // self.denom
 
+    @property
+    def is_rational(self):
+        """Say whether this element represents a rational number."""
+        return self.col[1:, :].is_zero_matrix
+
+    @property
+    def generator(self):
+        """
+        Return a :py:class:`~.Symbol` to be used when expressing this element
+        as a polynomial.
+
+        If we have an associated :py:class:`~.AlgebraicField` whose primitive
+        element has an alias symbol, we use that. Otherwise we use the variable
+        of the minimal polynomial defining the power basis to which we belong.
+        """
+        K = self.module.number_field
+        return K.ext.alias if K and K.ext.is_aliased else self.T.gen
+
+    def as_expr(self, x=None):
+        """Create a Basic expression from ``self``. """
+        return self.poly(x or self.generator).as_expr()
+
     def norm(self, T=None):
         """Compute the norm of this number."""
         T = T or self.T
diff --git a/sympy/polys/numberfields/primes.py b/sympy/polys/numberfields/primes.py
index dd3f6145cbfc..7fed6672cb3e 100644
--- a/sympy/polys/numberfields/primes.py
+++ b/sympy/polys/numberfields/primes.py
@@ -70,7 +70,12 @@ def __init__(self, ZK, p, alpha, f, e=None):
         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):
+    def __str__(self):
+        if self.alpha.is_rational:
+            return f'({self.p})'
+        return f'({self.p}, {self.alpha.as_expr()})'
+
+    def repr(self, field_gen=None, just_gens=False):
         """
         Print a representation of this prime ideal.
 
@@ -82,11 +87,11 @@ def pretty(self, field_gen=None, just_gens=False):
         >>> T = cyclotomic_poly(7, x)
         >>> K = QQ.algebraic_field((T, zeta))
         >>> P = K.primes_above(11)
-        >>> print(P[0].pretty())
+        >>> print(P[0].repr())
         [ (11, x**3 + 5*x**2 + 4*x - 1) e=1, f=3 ]
-        >>> print(P[0].pretty(field_gen=zeta))
+        >>> print(P[0].repr(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))
+        >>> print(P[0].repr(field_gen=zeta, just_gens=True))
         (11, zeta**3 + 5*zeta**2 + 4*zeta - 1)
 
         Parameters
@@ -114,7 +119,7 @@ def pretty(self, field_gen=None, just_gens=False):
         return f'[ {gens} e={e}, f={f} ]'
 
     def __repr__(self):
-        return self.pretty()
+        return self.repr()
 
     def as_submodule(self):
         r"""
diff --git a/sympy/polys/numberfields/tests/test_primes.py b/sympy/polys/numberfields/tests/test_primes.py
index 09e31e012918..93ec01375564 100644
--- a/sympy/polys/numberfields/tests/test_primes.py
+++ b/sympy/polys/numberfields/tests/test_primes.py
@@ -193,17 +193,16 @@ def test_decomp_8():
         x ** 3 + 9 * x ** 2 + 6 * x - 8,
         x ** 3 + 15 * x ** 2 - 9 * x + 13,
     )
-    '''
     def display(T, p, radical, P, I, J):
         """Useful for inspection, when running test manually."""
         print('=' * 20)
         print(T, p, radical)
         for Pi in P:
-            print(f'  ({Pi.pretty()})')
+            print(f'  ({Pi!r})')
         print("I: ", I)
         print("J: ", J)
         print(f'Equal: {I == J}')
-    '''
+    inspect = False
     for g in cases:
         T = Poly(g)
         rad = {}
@@ -216,7 +215,8 @@ def display(T, p, radical, P, I, J):
             P = prime_decomp(p, T, dK=dK, ZK=ZK, radical=radical)
             I = prod(Pi**Pi.e for Pi in P)
             J = p * ZK
-            #display(T, p, radical, P, I, J)
+            if inspect:
+                display(T, p, radical, P, I, J)
             assert I == J
 
 
@@ -238,16 +238,31 @@ def test_PrimeIdeal_add():
     assert P0 + 7 * P0.ZK == P0.as_submodule()
 
 
-def test_pretty_printing():
-    d = -7
-    T = Poly(x ** 2 - d)
-    rad = {}
-    ZK, dK = round_two(T, radicals=rad)
-    p = 2
-    P = prime_decomp(p, T, dK=dK, ZK=ZK, radical=rad.get(p))
+def test_str():
+    # Without alias:
+    k = QQ.alg_field_from_poly(Poly(x**2 + 7))
+    frp = k.primes_above(2)[0]
+    assert str(frp) == '(2, 3*_x/2 + 1/2)'
+
+    frp = k.primes_above(3)[0]
+    assert str(frp) == '(3)'
+
+    # With alias:
+    k = QQ.alg_field_from_poly(Poly(x ** 2 + 7), alias='alpha')
+    frp = k.primes_above(2)[0]
+    assert str(frp) == '(2, 3*alpha/2 + 1/2)'
+
+    frp = k.primes_above(3)[0]
+    assert str(frp) == '(3)'
+
+
+def test_repr():
+    T = Poly(x**2 + 7)
+    ZK, dK = round_two(T)
+    P = prime_decomp(2, T, dK=dK, ZK=ZK)
     assert repr(P[0]) == '[ (2, (3*x + 1)/2) e=1, f=1 ]'
-    assert P[0].pretty(field_gen=theta) == '[ (2, (3*theta + 1)/2) e=1, f=1 ]'
-    assert P[0].pretty(field_gen=theta, just_gens=True) == '(2, (3*theta + 1)/2)'
+    assert P[0].repr(field_gen=theta) == '[ (2, (3*theta + 1)/2) e=1, f=1 ]'
+    assert P[0].repr(field_gen=theta, just_gens=True) == '(2, (3*theta + 1)/2)'
 
 
 def test_PrimeIdeal_reduce_poly():
diff --git a/sympy/printing/latex.py b/sympy/printing/latex.py
index d8c4ebbf67dd..99f8579b24e7 100644
--- a/sympy/printing/latex.py
+++ b/sympy/printing/latex.py
@@ -630,6 +630,13 @@ def _print_AlgebraicNumber(self, expr):
         else:
             return self._print(expr.as_expr())
 
+    def _print_PrimeIdeal(self, expr):
+        p = self._print(expr.p)
+        if expr.alpha.is_rational:
+            return rf'\left({p}\right)'
+        alpha = self._print(expr.alpha.as_expr())
+        return rf'\left({p}, {alpha}\right)'
+
     def _print_Pow(self, expr):
         # Treat x**Rational(1,n) as special case
         if expr.exp.is_Rational and abs(expr.exp.p) == 1 and expr.exp.q != 1 \
diff --git a/sympy/printing/tests/test_latex.py b/sympy/printing/tests/test_latex.py
index a158cdabcad4..de8e402c5802 100644
--- a/sympy/printing/tests/test_latex.py
+++ b/sympy/printing/tests/test_latex.py
@@ -206,6 +206,11 @@ def test_latex_basic():
         r"\zeta^{3} + 2 \zeta^{2} + 3 \zeta + 4"
     assert latex(k.ext.field_element([1, 2, 3, 4]), order='old') == \
         r"4 + 3 \zeta + 2 \zeta^{2} + \zeta^{3}"
+    assert latex(k.primes_above(19)[0]) == \
+        r"\left(19, \zeta^{2} + 5 \zeta + 1\right)"
+    assert latex(k.primes_above(19)[0], order='old') == \
+           r"\left(19, 1 + 5 \zeta + \zeta^{2}\right)"
+    assert latex(k.primes_above(7)[0]) == r"\left(7\right)"
 
     assert latex(1.5e20*x) == r"1.5 \cdot 10^{20} x"
     assert latex(1.5e20*x, mul_symbol='dot') == r"1.5 \cdot 10^{20} \cdot x"


Code review comments:
diff_hunk: @@ -678,15 +702,30 @@ def __init__(self, T):
         Parameters
         ==========
 
-        T : :py:class:`~.Poly`
-            The monic, irreducible, univariate polynomial over :ref:`ZZ`, a
-            root of which is the generator of the power basis.
+        T : :py:class:`~.Poly`, :py:class:`~.AlgebraicField`
+            Either (1) the monic, irreducible, univariate polynomial over
+            :ref:`ZZ`, a root of which is the generator of the power basis,
+            or (2) an :py:class:`~.AlgebraicField` whose primitive element
+            is the generator of the power basis.
 
         """
+        K = None
+        if isinstance(T, AlgebraicField):
+            K, T = T, T.ext.minpoly_of_element()
+        if T.domain == QQ:
+            try:
+                T = Poly(T, domain=ZZ)
path: sympy/polys/numberfields/modules.py
text: @user1:
The direct way to do this is `T.set_domain(ZZ)`.

@user1:
I guess it isn't necessary to check if the domain is `QQ` since it should also raise for e.g. `ZZ[x]` if it can't be coerced to `ZZ`.

Another possibility would be to use `clear_denoms` I guess.

@author:
`set_domain()` looks good for this. Thanks, was not aware of this method!