Repository: sympy/sympy

Issue:
Some functions misinterpret `AlgebraicNumber`s.
In certain cases, the functions `minimal_polynomial()` and `field_isomorphism()`
fail to operate on the numbers that are actually represented by instances of
`AlgebraicNumber`.


## Explanation

An `AlgebraicNumber` instance represents an element `alpha` of a number field `Q(theta)`. The primitive
element `theta` for the number field is represented by its minimal polynomial, and a particular root thereof.
The element `alpha` is represented as a polynomial in `theta`.

The `minimal_polynomial()` and `field_isomorphism()` functions are treating `AlgebraicNumber`
instances as if they represent not `alpha` but `theta`. In other words, these functions are operating
on the primitive element of the number field, not on the represented element itself.


## Example 1, `minimal_polynomial()`

```python
>>> from sympy import AlgebraicNumber, minimal_polynomial, sqrt, S
>>> from sympy.abc import x
>>> theta = sqrt(2) + sqrt(3)
>>> a = AlgebraicNumber(theta, [S(1) / 2, 0, S(-9) / 2, 0])
>>> a.as_expr()
sqrt(2)
```

The `AlgebraicNumber` `a` represents the number `sqrt(2)`, as an element of the number field `Q(theta)`.

(The coefficients `[1/2, 0, -9/2, 0]` in the definition of `a` correspond to the fact that `sqrt(2) = theta**3/2 - 9*theta/2`.)

But when we ask for the minimal polynomial of `a`, expecting `x**2 - 2`, we instead get the minimal polynomial for `theta`:


```python
>>> minimal_polynomial(a, x)
x**4 - 10*x**2 + 1
>>> minimal_polynomial(theta, x)
x**4 - 10*x**2 + 1
```


## Example 2, `field_isomorphism()`

```python
>>> from sympy import to_number_field, field_isomorphism, sqrt
>>> theta = sqrt(2) + sqrt(3)
>>> eta   = sqrt(2) + sqrt(5)
>>> a = to_number_field(sqrt(2), theta)
>>> b = to_number_field(sqrt(2), eta)
>>> a.as_expr()
sqrt(2)
>>> b.as_expr()
sqrt(2)
>>> print(field_isomorphism(a, b))
None
```

Here, `a` and `b` are `AlgebraicNumber`s representing the number `sqrt(2)` in the fields
`Q(theta)` and `Q(eta)`, respectively.

When we ask whether there is a field isomorphism from `Q(a)` to `Q(b)` we are asking
for an embedding from `Q(sqrt(2))` into `Q(sqrt(2))`. The answer should be `[1, 0]`,
meaning we can map `a` to `1*b + 0`. But instead we get `None`, because what `field_isomorphism()`
is actually telling us is that there is no embedding from `Q(theta)` into `Q(eta)` (which is true, but
it's not the question we asked).



Discussion:

Patch:
diff --git a/doc/src/modules/core.rst b/doc/src/modules/core.rst
index c6fbe3ec8f..0ade07ed6d 100644
--- a/doc/src/modules/core.rst
+++ b/doc/src/modules/core.rst
@@ -129,6 +129,8 @@ AlgebraicNumber
 .. autoclass:: AlgebraicNumber
    :members:
 
+   .. automethod:: AlgebraicNumber.__new__
+
 NumberSymbol
 ^^^^^^^^^^^^
 .. autoclass:: NumberSymbol
diff --git a/sympy/core/numbers.py b/sympy/core/numbers.py
index 3b0bca77f5..cd7bdb2c1e 100644
--- a/sympy/core/numbers.py
+++ b/sympy/core/numbers.py
@@ -2509,9 +2509,24 @@ def __invert__(self):
 
 
 class AlgebraicNumber(Expr):
-    """Class for representing algebraic numbers in SymPy. """
+    r"""
+    Class for representing algebraic numbers in SymPy.
+
+    Symbolically, an instance of this class represents an element
+    $\alpha \in \mathbb{Q}(\theta) \hookrightarrow \mathbb{C}$. That is, the
+    algebraic number $\alpha$ is represented as an element of a particular
+    number field $\mathbb{Q}(\theta)$, with a particular embedding of this
+    field into the complex numbers.
+
+    Formally, the primitive element $\theta$ is given by two data points: (1)
+    its minimal polynomial (which defines $\mathbb{Q}(\theta)$), and (2) a
+    particular complex number that is a root of this polynomial (which defines
+    the embedding $\mathbb{Q}(\theta) \hookrightarrow \mathbb{C}$). Finally,
+    the algebraic number $\alpha$ which we represent is then given by the
+    coefficients of a polynomial in $\theta$.
+    """
 
-    __slots__ = ('rep', 'root', 'alias', 'minpoly')
+    __slots__ = ('rep', 'root', 'alias', 'minpoly', '_own_minpoly')
 
     is_AlgebraicNumber = True
     is_algebraic = True
@@ -2526,11 +2541,180 @@ class AlgebraicNumber(Expr):
     free_symbols: tSet[Basic] = set()
 
     def __new__(cls, expr, coeffs=None, alias=None, **args):
-        """Construct a new algebraic number. """
+        r"""
+        Construct a new algebraic number $\alpha$ belonging to a number field
+        $k = \mathbb{Q}(\theta)$.
+
+        There are four instance attributes to be determined:
+
+        ===========  ============================================================================
+        Attribute    Type/Meaning
+        ===========  ============================================================================
+        ``root``     :py:class:`~.Expr` for $\theta$ as a complex number
+        ``minpoly``  :py:class:`~.Poly`, the minimal polynomial of $\theta$
+        ``rep``      :py:class:`~sympy.polys.polyclasses.DMP` giving $\alpha$ as poly in $\theta$
+        ``alias``    :py:class:`~.Symbol` for $\theta$, or ``None``
+        ===========  ============================================================================
+
+        See Parameters section for how they are determined.
+
+        Parameters
+        ==========
+
+        expr : :py:class:`~.Expr`, or pair $(m, r)$
+            There are three distinct modes of construction, depending on what
+            is passed as *expr*.
+
+            **(1)** *expr* is an :py:class:`~.AlgebraicNumber`:
+            In this case we begin by copying all four instance attributes from
+            *expr*. If *coeffs* were also given, we compose the two coeff
+            polynomials (see below). If an *alias* was given, it overrides.
+
+            **(2)** *expr* is any other type of :py:class:`~.Expr`:
+            Then ``root`` will equal *expr*. Therefore it
+            must express an algebraic quantity, and we will compute its
+            ``minpoly``.
+
+            **(3)** *expr* is an ordered pair $(m, r)$ giving the
+            ``minpoly`` $m$, and a ``root`` $r$ thereof, which together
+            define $\theta$. In this case $m$ may be either a univariate
+            :py:class:`~.Poly` or any :py:class:`~.Expr` which represents the
+            same, while $r$ must be some :py:class:`~.Expr` representing a
+            complex number that is a root of $m$, including both explicit
+            expressions in radicals, and instances of
+            :py:class:`~.ComplexRootOf` or :py:class:`~.AlgebraicNumber`.
+
+        coeffs : list, :py:class:`~.ANP`, None, optional (default=None)
+            This defines ``rep``, giving the algebraic number $\alpha$ as a
+            polynomial in $\theta$.
+
+            If a list, the elements should be integers or rational numbers.
+            If an :py:class:`~.ANP`, we take its coefficients (using its
+            :py:meth:`~.ANP.to_list()` method). If ``None``, then the list of
+            coefficients defaults to ``[1, 0]``, meaning that $\alpha = \theta$
+            is the primitive element of the field.
+
+            If *expr* was an :py:class:`~.AlgebraicNumber`, let $g(x)$ be its
+            ``rep`` polynomial, and let $f(x)$ be the polynomial defined by
+            *coeffs*. Then ``self.rep`` will represent the composition
+            $(f \circ g)(x)$.
+
+        alias : str, :py:class:`~.Symbol`, None, optional (default=None)
+            This is a way to provide a name for the primitive element. We
+            described several ways in which the *expr* argument can define the
+            value of the primitive element, but none of these methods gave it
+            a name. Here, for example, *alias* could be set as
+            ``Symbol('theta')``, in order to make this symbol appear when
+            $\alpha$ is printed, or rendered as a polynomial, using the
+            :py:meth:`~.as_poly()` method.
+
+        Examples
+        ========
+
+        Recall that we are constructing an algebraic number as a field element
+        $\alpha \in \mathbb{Q}(\theta)$.
+
+        >>> from sympy import AlgebraicNumber, sqrt, CRootOf, S
+        >>> from sympy.abc import x
+
+        Example (1): $\alpha = \theta = \sqrt{2}$
+
+        >>> a1 = AlgebraicNumber(sqrt(2))
+        >>> a1.minpoly_of_element().as_expr(x)
+        x**2 - 2
+        >>> a1.evalf(10)
+        1.414213562
+
+        Example (2): $\alpha = 3 \sqrt{2} - 5$, $\theta = \sqrt{2}$. We can
+        either build on the last example:
+
+        >>> a2 = AlgebraicNumber(a1, [3, -5])
+        >>> a2.as_expr()
+        -5 + 3*sqrt(2)
+
+        or start from scratch:
+
+        >>> a2 = AlgebraicNumber(sqrt(2), [3, -5])
+        >>> a2.as_expr()
+        -5 + 3*sqrt(2)
+
+        Example (3): $\alpha = 6 \sqrt{2} - 11$, $\theta = \sqrt{2}$. Again we
+        can build on the previous example, and we see that the coeff polys are
+        composed:
+
+        >>> a3 = AlgebraicNumber(a2, [2, -1])
+        >>> a3.as_expr()
+        -11 + 6*sqrt(2)
+
+        reflecting the fact that $(2x - 1) \circ (3x - 5) = 6x - 11$.
+
+        Example (4): $\alpha = \sqrt{2}$, $\theta = \sqrt{2} + \sqrt{3}$. The
+        easiest way is to use the :py:func:`~.to_number_field()` function:
+
+        >>> from sympy import to_number_field
+        >>> a4 = to_number_field(sqrt(2), sqrt(2) + sqrt(3))
+        >>> a4.minpoly_of_element().as_expr(x)
+        x**2 - 2
+        >>> a4.to_root()
+        sqrt(2)
+        >>> a4.primitive_element()
+        sqrt(2) + sqrt(3)
+        >>> a4.coeffs()
+        [1/2, 0, -9/2, 0]
+
+        but if you already knew the right coefficients, you could construct it
+        directly:
+
+        >>> a4 = AlgebraicNumber(sqrt(2) + sqrt(3), [S(1)/2, 0, S(-9)/2, 0])
+        >>> a4.to_root()
+        sqrt(2)
+        >>> a4.primitive_element()
+        sqrt(2) + sqrt(3)
+
+        Example (5): Construct the Golden Ratio as an element of the 5th
+        cyclotomic field, supposing we already know its coefficients. This time
+        we introduce the alias $\zeta$ for the primitive element of the field:
+
+        >>> from sympy import cyclotomic_poly
+        >>> from sympy.abc import zeta
+        >>> a5 = AlgebraicNumber(CRootOf(cyclotomic_poly(5), -1),
+        ...                  [-1, -1, 0, 0], alias=zeta)
+        >>> a5.as_poly().as_expr()
+        -zeta**3 - zeta**2
+        >>> a5.evalf()
+        1.61803398874989
+
+        (The index ``-1`` to ``CRootOf`` selects the complex root with the
+        largest real and imaginary parts, which in this case is
+        $\mathrm{e}^{2i\pi/5}$. See :py:class:`~.ComplexRootOf`.)
+
+        Example (6): Building on the last example, construct the number
+        $2 \phi \in \mathbb{Q}(\phi)$, where $\phi$ is the Golden Ratio:
+
+        >>> from sympy.abc import phi
+        >>> a6 = AlgebraicNumber(a5.to_root(), coeffs=[2, 0], alias=phi)
+        >>> a6.as_poly().as_expr()
+        2*phi
+        >>> a6.primitive_element().evalf()
+        1.61803398874989
+
+        Note that we needed to use ``a5.to_root()``, since passing ``a5`` as
+        the first argument would have constructed the number $2 \phi$ as an
+        element of the field $\mathbb{Q}(\zeta)$:
+
+        >>> a6_wrong = AlgebraicNumber(a5, coeffs=[2, 0])
+        >>> a6_wrong.as_poly().as_expr()
+        -2*zeta**3 - 2*zeta**2
+        >>> a6_wrong.primitive_element().evalf()
+        0.309016994374947 + 0.951056516295154*I
+
+        """
         from sympy.polys.polyclasses import ANP, DMP
         from sympy.polys.numberfields import minimal_polynomial
 
         expr = sympify(expr)
+        rep0 = None
+        alias0 = None
 
         if isinstance(expr, (tuple, Tuple)):
             minpoly, root = expr
@@ -2539,7 +2723,8 @@ def __new__(cls, expr, coeffs=None, alias=None, **args):
                 from sympy.polys.polytools import Poly
                 minpoly = Poly(minpoly)
         elif expr.is_AlgebraicNumber:
-            minpoly, root = expr.minpoly, expr.root
+            minpoly, root, rep0, alias0 = (expr.minpoly, expr.root,
+                                           expr.rep, expr.alias)
         else:
             minpoly, root = minimal_polynomial(
                 expr, args.get('gen'), polys=True), expr
@@ -2554,15 +2739,22 @@ def __new__(cls, expr, coeffs=None, alias=None, **args):
                 rep = DMP.from_list(coeffs.to_list(), 0, dom)
                 scoeffs = Tuple(*coeffs.to_list())
 
-            if rep.degree() >= minpoly.degree():
-                rep = rep.rem(minpoly.rep)
-
         else:
             rep = DMP.from_list([1, 0], 0, dom)
             scoeffs = Tuple(1, 0)
 
+        if rep0 is not None:
+            from sympy.polys.densetools import dup_compose
+            c = dup_compose(rep.rep, rep0.rep, dom)
+            rep = DMP.from_list(c, 0, dom)
+            scoeffs = Tuple(*c)
+
+        if rep.degree() >= minpoly.degree():
+            rep = rep.rem(minpoly.rep)
+
         sargs = (root, scoeffs)
 
+        alias = alias or alias0
         if alias is not None:
             from .symbol import Symbol
             if not isinstance(alias, Symbol):
@@ -2576,6 +2768,8 @@ def __new__(cls, expr, coeffs=None, alias=None, **args):
         obj.alias = alias
         obj.minpoly = minpoly
 
+        obj._own_minpoly = None
+
         return obj
 
     def __hash__(self):
@@ -2641,6 +2835,184 @@ def _eval_simplify(self, **kwargs):
                     return AlgebraicNumber(r)
         return self
 
+    @property
+    def is_primitive_element(self):
+        r"""
+        Say whether this algebraic number $\alpha \in \mathbb{Q}(\theta)$ is
+        equal to the primitive element $\theta$ for its field.
+        """
+        c = self.coeffs()
+        # Second case occurs if self.minpoly is linear:
+        return c == [1, 0] or c == [self.root]
+
+    def primitive_element(self):
+        r"""
+        Get the primitive element $\theta$ for the number field
+        $\mathbb{Q}(\theta)$ to which this algebraic number $\alpha$ belongs.
+
+        Returns
+        =======
+
+        AlgebraicNumber
+
+        """
+        if self.is_primitive_element:
+            return self
+        return AlgebraicNumber((self.minpoly, self.root), coeffs=[1, 0])
+
+    def to_primitive_element(self, radicals=True):
+        r"""
+        Convert ``self`` to an :py:class:`~.AlgebraicNumber` instance that is
+        equal to its own primitive element.
+
+        Explanation
+        ===========
+
+        If we represent $\alpha \in \mathbb{Q}(\theta)$, $\alpha \neq \theta$,
+        construct a new :py:class:`~.AlgebraicNumber` that represents
+        $\alpha \in \mathbb{Q}(\alpha)$.
+
+        Examples
+        ========
+
+        >>> from sympy import sqrt, to_number_field
+        >>> from sympy.abc import x
+        >>> a = to_number_field(sqrt(2), sqrt(2) + sqrt(3))
+
+        The :py:class:`~.AlgebraicNumber` ``a`` represents the number
+        $\sqrt{2}$ in the field $\mathbb{Q}(\sqrt{2} + \sqrt{3})$. Rendering
+        ``a`` as a polynomial,
+
+        >>> a.as_poly().as_expr(x)
+        x**3/2 - 9*x/2
+
+        reflects the fact that $\sqrt{2} = \theta^3/2 - 9 \theta/2$, where
+        $\theta = \sqrt{2} + \sqrt{3}$.
+
+        ``a`` is not equal to its own primitive element. Its minpoly
+
+        >>> a.minpoly.as_poly().as_expr(x)
+        x**4 - 10*x**2 + 1
+
+        is that of $\theta$.
+
+        Converting to a primitive element,
+
+        >>> a_prim = a.to_primitive_element()
+        >>> a_prim.minpoly.as_poly().as_expr(x)
+        x**2 - 2
+
+        we obtain an :py:class:`~.AlgebraicNumber` whose ``minpoly`` is that of
+        the number itself.
+
+        Parameters
+        ==========
+
+        radicals : boolean, optional (default=True)
+            If ``True``, then we will try to return an
+            :py:class:`~.AlgebraicNumber` whose ``root`` is an expression
+            in radicals. If that is not possible (or if *radicals* is
+            ``False``), ``root`` will be a :py:class:`~.ComplexRootOf`.
+
+        Returns
+        =======
+
+        AlgebraicNumber
+
+        See Also
+        ========
+
+        is_primitive_element
+
+        """
+        if self.is_primitive_element:
+            return self
+        m = self.minpoly_of_element()
+        r = self.to_root(radicals=radicals)
+        return AlgebraicNumber((m, r))
+
+    def minpoly_of_element(self):
+        r"""
+        Compute the minimal polynomial for this algebraic number.
+
+        Explanation
+        ===========
+
+        Recall that we represent an element $\alpha \in \mathbb{Q}(\theta)$.
+        Our instance attribute ``self.minpoly`` is the minimal polynomial for
+        our primitive element $\theta$. This method computes the minimal
+        polynomial for $\alpha$.
+
+        """
+        if self._own_minpoly is None:
+            if self.is_primitive_element:
+                self._own_minpoly = self.minpoly
+            else:
+                from sympy.polys.numberfields.minpoly import minpoly
+                theta = self.primitive_element()
+                self._own_minpoly = minpoly(self.as_expr(theta), polys=True)
+        return self._own_minpoly
+
+    def to_root(self, radicals=True, minpoly=None):
+        """
+        Convert to an :py:class:`~.Expr` that is not an
+        :py:class:`~.AlgebraicNumber`, specifically, either a
+        :py:class:`~.ComplexRootOf`, or, optionally and where possible, an
+        expression in radicals.
+
+        Parameters
+        ==========
+
+        radicals : boolean, optional (default=True)
+            If ``True``, then we will try to return the root as an expression
+            in radicals. If that is not possible, we will return a
+            :py:class:`~.ComplexRootOf`.
+
+        minpoly : :py:class:`~.Poly`
+            If the minimal polynomial for `self` has been pre-computed, it can
+            be passed in order to save time.
+
+        """
+        if self.is_primitive_element and not isinstance(self.root, AlgebraicNumber):
+            return self.root
+        m = minpoly or self.minpoly_of_element()
+        roots = m.all_roots(radicals=radicals)
+        if len(roots) == 1:
+            return roots[0]
+        root = None
+        if all(hasattr(r, "_get_interval") for r in roots):
+            root = self._to_root_by_intervals(roots)
+        if root is not None:
+            return root
+        return self._to_root_by_distance(roots)
+
+    def _to_root_by_intervals(self, roots):
+        intervals = [r._get_interval() for r in roots]
+        D0 = int(max(i.max_denom for i in intervals))
+        # Make n more than the number of decimal places in D0. This is to
+        # eliminate false positives, i.e. cases where we appear to belong to
+        # an interval but only due to rounding errors.
+        n = math.ceil(D0.bit_length()/3.3) + 2
+        c = self.evalf(n).as_real_imag()
+        for j, i in enumerate(intervals):
+            if c in i:
+                return roots[j]
+        return None
+
+    def _to_root_by_distance(self, roots, max_prec=160):
+        # Compare sympy.polys.numberfields.minpoly._choose_factor()
+        prec1 = 10
+        while prec1 <= max_prec:
+            r0 = self.evalf(prec1)
+            candidates = [(abs(r0 - r.evalf(prec1)), j)
+                          for j, r in enumerate(roots)]
+            can = sorted(candidates)
+            (a, ix), (b, _) = can[:2]
+            if b > a * 10 ** 6:
+                return roots[ix]
+            prec1 *= 2
+        raise NotImplementedError("Could not locate root.")
+
 
 class RationalConstant(Rational):
     """
diff --git a/sympy/polys/numberfields/minpoly.py b/sympy/polys/numberfields/minpoly.py
index 797a2af205..fa0265a0a5 100644
--- a/sympy/polys/numberfields/minpoly.py
+++ b/sympy/polys/numberfields/minpoly.py
@@ -808,10 +808,10 @@ def bottom_up_scan(ex):
                 else:
                     return symbols[expr]
         elif ex.is_AlgebraicNumber:
-            if ex.root not in mapping:
-                return update_mapping(ex.root, ex.minpoly)
+            if ex not in mapping:
+                return update_mapping(ex, ex.minpoly_of_element())
             else:
-                return symbols[ex.root]
+                return symbols[ex]
 
         raise NotAlgebraic("%s doesn't seem to be an algebraic number" % ex)
 
@@ -840,7 +840,7 @@ def simpler_inverse(ex):
     inverted = False
     ex = expand_multinomial(ex)
     if ex.is_AlgebraicNumber:
-        return ex.minpoly.as_expr(x)
+        return ex.minpoly_of_element().as_expr(x)
     elif ex.is_Rational:
         result = ex.q*x - ex.p
     else:
diff --git a/sympy/polys/numberfields/subfield.py b/sympy/polys/numberfields/subfield.py
index cc5e383a3a..dc091c5953 100644
--- a/sympy/polys/numberfields/subfield.py
+++ b/sympy/polys/numberfields/subfield.py
@@ -209,6 +209,9 @@ def field_isomorphism(a, b, *, fast=True):
     if not b.is_AlgebraicNumber:
         b = AlgebraicNumber(b)
 
+    a = a.to_primitive_element()
+    b = b.to_primitive_element()
+
     if a == b:
         return a.coeffs()
 
diff --git a/sympy/polys/rootisolation.py b/sympy/polys/rootisolation.py
index 9bb0647235..cedfbf6dd0 100644
--- a/sympy/polys/rootisolation.py
+++ b/sympy/polys/rootisolation.py
@@ -1732,6 +1732,11 @@ def center(self):
         """Return the center of the real isolating interval. """
         return (self.a + self.b)/2
 
+    @property
+    def max_denom(self):
+        """Return the largest denominator occurring in either endpoint. """
+        return max(self.a.denominator, self.b.denominator)
+
     def as_tuple(self):
         """Return tuple representation of real isolating interval. """
         return (self.a, self.b)
@@ -1739,6 +1744,24 @@ def as_tuple(self):
     def __repr__(self):
         return "(%s, %s)" % (self.a, self.b)
 
+    def __contains__(self, item):
+        """
+        Say whether a complex number belongs to this real interval.
+
+        Parameters
+        ==========
+
+        item : pair (re, im) or number re
+            Either a pair giving the real and imaginary parts of the number,
+            or else a real number.
+
+        """
+        if isinstance(item, tuple):
+            re, im = item
+        else:
+            re, im = item, 0
+        return im == 0 and self.a <= re <= self.b
+
     def is_disjoint(self, other):
         """Return ``True`` if two isolation intervals are disjoint. """
         if isinstance(other, RealInterval):
@@ -1989,6 +2012,12 @@ def center(self):
         """Return the center of the complex isolating interval. """
         return ((self.ax + self.bx)/2, (self.ay + self.by)/2)
 
+    @property
+    def max_denom(self):
+        """Return the largest denominator occurring in either endpoint. """
+        return max(self.ax.denominator, self.bx.denominator,
+                   self.ay.denominator, self.by.denominator)
+
     def as_tuple(self):
         """Return tuple representation of the complex isolating
         interval's SW and NE corners, respectively. """
@@ -2002,6 +2031,25 @@ def conjugate(self):
         return ComplexInterval(self.a, self.b, self.I, self.Q,
             self.F1, self.F2, self.f1, self.f2, self.dom, conj=True)
 
+    def __contains__(self, item):
+        """
+        Say whether a complex number belongs to this complex rectangular
+        region.
+
+        Parameters
+        ==========
+
+        item : pair (re, im) or number re
+            Either a pair giving the real and imaginary parts of the number,
+            or else a real number.
+
+        """
+        if isinstance(item, tuple):
+            re, im = item
+        else:
+            re, im = item, 0
+        return self.ax <= re <= self.bx and self.ay <= im <= self.by
+
     def is_disjoint(self, other):
         """Return ``True`` if two isolation intervals are disjoint. """
         if isinstance(other, RealInterval):
diff --git a/sympy/polys/rootoftools.py b/sympy/polys/rootoftools.py
index e2ae271cfd..e015aec9d1 100644
--- a/sympy/polys/rootoftools.py
+++ b/sympy/polys/rootoftools.py
@@ -167,7 +167,12 @@ class ComplexRootOf(RootOf):
     """Represents an indexed complex root of a polynomial.
 
     Roots of a univariate polynomial separated into disjoint
-    real or complex intervals and indexed in a fixed order.
+    real or complex intervals and indexed in a fixed order:
+
+    * real roots come first and are sorted in increasing order;
+    * complex roots come next and are sorted primarily by increasing
+      real part, secondarily by increasing imaginary part.
+
     Currently only rational coefficients are allowed.
     Can be imported as ``CRootOf``. To avoid confusion, the
     generator must be a Symbol.


Test patch:
diff --git a/sympy/polys/numberfields/tests/test_minpoly.py b/sympy/polys/numberfields/tests/test_minpoly.py
index ee202b9146..39fd55e2ee 100644
--- a/sympy/polys/numberfields/tests/test_minpoly.py
+++ b/sympy/polys/numberfields/tests/test_minpoly.py
@@ -427,6 +427,13 @@ def test_issue_22559():
     assert minimal_polynomial(alpha**3, x) == x**2 - 8
 
 
+def test_issue_22561():
+    a = AlgebraicNumber(sqrt(2) + sqrt(3), [S(1) / 2, 0, S(-9) / 2, 0], gen=x)
+    assert a.as_expr() == sqrt(2)
+    assert minimal_polynomial(a, x) == x**2 - 2
+    assert minimal_polynomial(a**3, x) == x**2 - 8
+
+
 def test_separate_sq_not_impl():
     raises(NotImplementedError, lambda: _separate_sq(x**(S(1)/3) + x))
 
diff --git a/sympy/polys/numberfields/tests/test_numbers.py b/sympy/polys/numberfields/tests/test_numbers.py
index 5292957a2c..f8f350719c 100644
--- a/sympy/polys/numberfields/tests/test_numbers.py
+++ b/sympy/polys/numberfields/tests/test_numbers.py
@@ -1,7 +1,7 @@
 """Tests on algebraic numbers. """
 
 from sympy.core.containers import Tuple
-from sympy.core.numbers import (AlgebraicNumber, Rational)
+from sympy.core.numbers import (AlgebraicNumber, I, Rational)
 from sympy.core.singleton import S
 from sympy.core.symbol import Symbol
 from sympy.functions.elementary.miscellaneous import sqrt
@@ -9,6 +9,7 @@
 from sympy.polys.numberfields.subfield import to_number_field
 from sympy.polys.polyclasses import DMP
 from sympy.polys.domains import QQ
+from sympy.polys.rootoftools import CRootOf
 from sympy.abc import x, y
 
 
@@ -149,6 +150,21 @@ def test_AlgebraicNumber():
     a = AlgebraicNumber(sqrt(2), [1, 2, 3])
     assert a.args == (sqrt(2), Tuple(1, 2, 3))
 
+    a = AlgebraicNumber(sqrt(2), [1, 2], "alpha")
+    b = AlgebraicNumber(a)
+    c = AlgebraicNumber(a, alias="gamma")
+    assert a == b
+    assert c.alias.name == "gamma"
+
+    a = AlgebraicNumber(sqrt(2) + sqrt(3), [S(1)/2, 0, S(-9)/2, 0])
+    b = AlgebraicNumber(a, [1, 0, 0])
+    assert b.root == a.root
+    assert a.to_root() == sqrt(2)
+    assert b.to_root() == 2
+
+    a = AlgebraicNumber(2)
+    assert a.is_primitive_element is True
+
 
 def test_to_algebraic_integer():
     a = AlgebraicNumber(sqrt(3), gen=x).to_algebraic_integer()
@@ -173,3 +189,14 @@ def test_to_algebraic_integer():
     assert a.minpoly == x**2 - 12
     assert a.root == 2*sqrt(3)
     assert a.rep == DMP([QQ(7, 19), QQ(3)], QQ)
+
+
+def test_AlgebraicNumber_to_root():
+    assert AlgebraicNumber(sqrt(2)).to_root() == sqrt(2)
+
+    zeta5_squared = AlgebraicNumber(CRootOf(x**5 - 1, 4), coeffs=[1, 0, 0])
+    assert zeta5_squared.to_root() == CRootOf(x**4 + x**3 + x**2 + x + 1, 1)
+
+    zeta3_squared = AlgebraicNumber(CRootOf(x**3 - 1, 2), coeffs=[1, 0, 0])
+    assert zeta3_squared.to_root() == -S(1)/2 - sqrt(3)*I/2
+    assert zeta3_squared.to_root(radicals=False) == CRootOf(x**2 + x + 1, 0)
diff --git a/sympy/polys/numberfields/tests/test_subfield.py b/sympy/polys/numberfields/tests/test_subfield.py
index 68cf230f51..918e19a694 100644
--- a/sympy/polys/numberfields/tests/test_subfield.py
+++ b/sympy/polys/numberfields/tests/test_subfield.py
@@ -281,3 +281,9 @@ def test_to_number_field():
     assert to_number_field(sqrt(2), AlgebraicNumber(sqrt(2) + sqrt(3))) == a
 
     raises(IsomorphismFailed, lambda: to_number_field(sqrt(2), sqrt(3)))
+
+
+def test_issue_22561():
+    a = to_number_field(sqrt(2), sqrt(2) + sqrt(3))
+    b = to_number_field(sqrt(2), sqrt(2) + sqrt(5))
+    assert field_isomorphism(a, b) == [1, 0]