Repository: sympy/sympy
Pull request: Bugfix in `minimal_polynomial()` and `field_isomorphism()`.
<!-- 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 #22561 

#### Brief description of what is fixed or changed

The `minimal_polynomial()` and `field_isomorphism()` functions were not always doing what they advertised to do.
This was due to a misapprehension of what an instance of the `AlgebraicNumber` class represents.

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 were written as if an `AlgebraicNumber`
instance represented not `alpha` but `theta`. In other words, these functions operated on the primitive
element of the number field, not on the represented element itself.

We repair this as follows:
* Provide new docstrings for the `AlgebraicNumber` class, explaining what it represents and how.
* Provide new methods for the `AlgebraicNumber` class, making it easy to:
  - get the minimal polynomial for alpha
  - test whether alpha == theta
  - convert into an instance were alpha does equal theta
* Update the `minimal_polynomial()` and `field_isomorphism()` functions to operate on alpha, not theta.
* Provide new unit tests that failed before the fix.

#### 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
  * Fixed related issues in `minimal_polynomial()` and `field_isomorphism()`.

<!-- END RELEASE NOTES -->


Problem statement:
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).

Hints from issue discussion:

Commit message: Bugfix in `minimal_polynomial()` and `field_isomorphism()`.

The `minimal_polynomial()` and `field_isomorphism()` functions were not always doing what they advertised to do.
This was due to a misapprehension of what an instance of the `AlgebraicNumber` class represents.

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 were written as if an `AlgebraicNumber`
instance represented not `alpha` but `theta`. In other words, these functions operated on the primitive
element of the number field, not on the represented element itself.

We repair this as follows:
* Provide new docstrings for the `AlgebraicNumber` class, explaining what it represents and how.
* Provide new methods for the `AlgebraicNumber` class, making it easy to:
  - get the minimal polynomial for alpha
  - test whether alpha == theta
  - convert into an instance were alpha does equal theta
* Update the `minimal_polynomial()` and `field_isomorphism()` functions to operate on alpha, not theta.
* Provide new unit tests that failed before the fix.

Patch to review:
diff --git a/doc/src/modules/core.rst b/doc/src/modules/core.rst
index c6fbe3ec8f59..0ade07ed6df3 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 3b0bca77f5e7..accff0a8e23c 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,7 +2541,97 @@ 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)$.
+
+        Parameters
+        ==========
+
+        expr : :py:class:`~.Expr`, or pair (m, r)
+            This defines the primitive element $\theta$ of the number field
+            $k$. If *expr* is an :py:class:`~.AlgebraicNumber`, then our
+            primitive element $\theta$ will be the same as that of *expr*.
+            If it is any other type of :py:class:`~.Expr`, then it itself will
+            be our primitive element. Therefore it must express an algebraic
+            quantity, and we will compute its minimal polynomial.
+            Otherwise *expr* must be an ordered pair
+            $(m, r)$ giving the minimal polynomial $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, and instances of
+            :py:class:`~.ComplexRootOf`.
+
+        coeffs : list, :py:class:`~.ANP`, None, optional (default=None)
+            This defines the algebraic number $\alpha$ as an element of $k$,
+            as a linear combination of falling powers of $\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.
+
+        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 rendered as a polynomial, using the
+            :py:meth:`~.as_poly()` method.
+
+        Examples
+        ========
+
+        >>> from sympy import AlgebraicNumber, sqrt, CRootOf, S
+        >>> from sympy.abc import x, theta
+
+        *expr* an explicit algebraic number, *coeffs* ``None``:
+
+        >>> a0 = AlgebraicNumber(sqrt(2) + sqrt(3))
+        >>> a0.minpoly_of_elt().as_expr(x)
+        x**4 - 10*x**2 + 1
+        >>> a0.n(10)
+        3.146264370
+
+        *expr* an explicit algebraic number, *coeffs* given:
+
+        >>> a1 = AlgebraicNumber(sqrt(2) + sqrt(3), [S(1)/2, 0, S(-9)/2, 0])
+        >>> a1.minpoly_of_elt().as_expr(x)
+        x**2 - 2
+        >>> a1.n(10)
+        1.414213562
+        >>> a1.primitive_elt()
+        sqrt(2) + sqrt(3)
+
+        *expr* an :py:class:`~.AlgebraicNumber` instance, *alias* provided:
+
+        >>> a2 = AlgebraicNumber(a0, [S(1)/2, 0, S(-9)/2, 0], alias=theta)
+        >>> a2.primitive_elt() == a0
+        True
+        >>> a2.as_expr() == a1.as_expr()
+        True
+        >>> a1.as_poly().as_expr()
+        _x**3/2 - 9*_x/2
+        >>> a2.as_poly().as_expr()
+        theta**3/2 - 9*theta/2
+
+        *expr* a pair (poly, explicit root):
+
+        >>> f = x**2 - x - 1
+        >>> a3 = AlgebraicNumber((f, (1 + sqrt(5))/2))
+        >>> a3.primitive_elt().n(10)
+        1.618033989
+
+        *expr* a pair (poly, implicit root):
+
+        >>> a4 = AlgebraicNumber((f, CRootOf(f, -1)))
+        >>> a4.primitive_elt().n(10)
+        1.618033989
+
+        """
         from sympy.polys.polyclasses import ANP, DMP
         from sympy.polys.numberfields import minimal_polynomial
 
@@ -2576,6 +2681,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 +2748,119 @@ def _eval_simplify(self, **kwargs):
                     return AlgebraicNumber(r)
         return self
 
+    @property
+    def is_primitive_elt(self):
+        r"""
+        Say whether this algebraic number $\alpha \in \mathbb{Q}(\theta)$ is
+        equal to the primitive element $\theta$ for its field.
+        """
+        return self.coeffs() == [1, 0]
+
+    def primitive_elt(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_elt:
+            return self
+        return AlgebraicNumber(self, coeffs=[1, 0])
+
+    def to_primitive_elt(self, prec=15):
+        r"""
+        Convert ``self`` to an :py:class:`~.AlgebraicNumber` instance that is
+        equal to its own primitive element.
+
+        Explanation
+        ===========
+
+        Since an :py:class:`~.AlgebraicNumber` stores both the minimal
+        polynomial, and a particular root value, for its primitive element,
+        it is sometimes more convenient to work with instances that equal their
+        own primitive element.
+
+        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_elt()
+        >>> 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
+        ==========
+
+        prec : int, optional (default=15)
+            Decimal places of precision for determining to which root of its
+            minimal polynomial this number is equal.
+
+        Returns
+        =======
+
+        AlgebraicNumber
+
+        See Also
+        ========
+
+        is_primitive_elt
+
+        """
+        if self.is_primitive_elt:
+            return self
+        m = self.minpoly_of_elt()
+        r0 = self.n(prec)
+        closest_root, min_distance = None, oo
+        for r in m.all_roots():
+            d = abs(r0 - r.n(prec))
+            if d < min_distance:
+                min_distance = d
+                closest_root = r
+        return AlgebraicNumber((m, closest_root))
+
+    def minpoly_of_elt(self):
+        r"""
+        Compute the minimal polynomial for the element
+        $\alpha \in \mathbb{Q}(\theta)$ we represent.
+        """
+        if self._own_minpoly is None:
+            if self.is_primitive_elt:
+                self._own_minpoly = self.minpoly
+            else:
+                from sympy.polys.numberfields.minpoly import minpoly
+                theta = self.primitive_elt()
+                self._own_minpoly = minpoly(self.as_expr(theta), polys=True)
+        return self._own_minpoly
+
 
 class RationalConstant(Rational):
     """
diff --git a/sympy/polys/numberfields/minpoly.py b/sympy/polys/numberfields/minpoly.py
index 797a2af20521..607365bc3510 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_elt())
             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_elt().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 cc5e383a3ad1..3705cfae98d3 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_elt()
+    b = b.to_primitive_elt()
+
     if a == b:
         return a.coeffs()
 
diff --git a/sympy/polys/numberfields/tests/test_minpoly.py b/sympy/polys/numberfields/tests/test_minpoly.py
index ee202b914646..39fd55e2ee7d 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_subfield.py b/sympy/polys/numberfields/tests/test_subfield.py
index 68cf230f5157..918e19a694a9 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]


Code review comments:
diff_hunk: @@ -2641,6 +2748,119 @@ def _eval_simplify(self, **kwargs):
                     return AlgebraicNumber(r)
         return self
 
+    @property
+    def is_primitive_elt(self):
+        r"""
+        Say whether this algebraic number $\alpha \in \mathbb{Q}(\theta)$ is
+        equal to the primitive element $\theta$ for its field.
+        """
+        return self.coeffs() == [1, 0]
+
+    def primitive_elt(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_elt:
+            return self
+        return AlgebraicNumber(self, coeffs=[1, 0])
+
+    def to_primitive_elt(self, prec=15):
+        r"""
+        Convert ``self`` to an :py:class:`~.AlgebraicNumber` instance that is
+        equal to its own primitive element.
+
+        Explanation
+        ===========
+
+        Since an :py:class:`~.AlgebraicNumber` stores both the minimal
+        polynomial, and a particular root value, for its primitive element,
+        it is sometimes more convenient to work with instances that equal their
+        own primitive element.
+
+        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_elt()
+        >>> 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
+        ==========
+
+        prec : int, optional (default=15)
+            Decimal places of precision for determining to which root of its
+            minimal polynomial this number is equal.
+
+        Returns
+        =======
+
+        AlgebraicNumber
+
+        See Also
+        ========
+
+        is_primitive_elt
+
+        """
+        if self.is_primitive_elt:
+            return self
+        m = self.minpoly_of_elt()
+        r0 = self.n(prec)
+        closest_root, min_distance = None, oo
+        for r in m.all_roots():
+            d = abs(r0 - r.n(prec))
+            if d < min_distance:
+                min_distance = d
+                closest_root = r
+        return AlgebraicNumber((m, closest_root))
path: sympy/core/numbers.py
text: @author:
@user1 Do you know if there's an existing function for this problem? I couldn't find one, so wrote an ad hoc solution.

You have a `Poly` `m`, and you have a number (in this case `self`), and you know that the number is equal to one of the roots of `m`. You want to determine which of the roots it is, and then you want to obtain an exact expression for that root (such as a `CRootOf`, or an explicit expression using radicals etc.).

@user1:
I'm not sure there is a function but one could be added to rootisolation.py. The key function is `dup_isolate_all_roots` which gives bounding intervals for each of the roots. Then evalf (i.e. .n) can tell us if our expression is in the bounding interval.

@user1:
This function is the only thing I'm not sure about

@user2:
Maybe this could be based on the same idea that is used in `_choose_factor`:
https://github.com/sympy/sympy/blob/551c16facb2081c7b7a085aeaf20316244a1da01/sympy/polys/numberfields/minpoly.py#L80-L86
All distances are computed and the two least ones are compared. It they have different orders of magnitude (probably even less than `10**6` would suffice), then the minimum is chosen, otherwise the precision in incremented in a loop. There is no need to give a precision parameter as an argument.

@author:
Wait...... First, thanks @user1 and @user2 for your ideas. I think you're on the right track for the problem I asked to solve. However, maybe it's not needed here.

What we're talking about could go in a method called `to_CRootOf()`, which would turn an `AlgebraicNumber` into a `CRootOf` equal to the same complex number.

But I think all we need here is any `Expr` that can be `evalf`-ed to arbitrary precision, and `self` already is that.
So the whole method becomes:

```python
def to_primitive_elt(self):
    if self.is_primitive_elt:
        return self
    m = self.minpoly_of_elt()
    return AlgebraicNumber((m, self))
```

Sorry for the confusion! I didn't see this earlier.

@author:
Cancel that. This turns out to be quite tricky.

## TL;DR

I think we were on the right track in the first place, and I'll push a new commit that does a better job of determining the root.

## Details

It seems that using one `AlgebraicNumber` as the `root` arg to another `AlgebraicNumber` is going to be very error prone and should probably be avoided.

As a demonstration of what can happen, let's start with an `AlgebraicNumber` `a` and its minpoly `m`:

```python
>>> from sympy import sqrt, to_number_field, AlgebraicNumber
>>> a = to_number_field(sqrt(2), sqrt(2) + sqrt(3))
>>> a.as_expr()
sqrt(2)
>>> m = a.minpoly_of_elt()
```

and use these as the `(m, r)` when constructing a new `AlgebraicNumber` `b`:

```python
>>> b = AlgebraicNumber((m, a))
>>> b.as_expr()
sqrt(2)
```

So far, so good.

But now, due to the funny way `AlgebraicNumber.__new__()` treats its first arg, when that arg is an `AlgebraicNumber`,

https://github.com/sympy/sympy/blob/519f1b75eda8830189b1385197f2c4fb49d333aa/sympy/core/numbers.py#L2541-L2542

`b` is not going to be picklable:

```python
>>> import pickle
>>> c = pickle.loads(pickle.dumps(b))
>>> c.as_expr()
sqrt(2) + sqrt(3)
```

Yipes!

But maybe we can repair pickling by changing this line:

https://github.com/sympy/sympy/blob/519f1b75eda8830189b1385197f2c4fb49d333aa/sympy/core/numbers.py#L2564

to

```python
        sargs = (expr, scoeffs)
```

which seems totally reasonable. Let our args be exactly as they were passed to us.

But this turns out to have too many far-flung ramifications.

For just one example, [this unit test](https://github.com/sympy/sympy/blob/519f1b75eda8830189b1385197f2c4fb49d333aa/sympy/polys/numberfields/tests/test_minpoly.py#L416-L422) is broken by it, because this method:

https://github.com/sympy/sympy/blob/519f1b75eda8830189b1385197f2c4fb49d333aa/sympy/polys/domains/gaussiandomains.py#L320-L323

expects the first arg of an `AlgebraicNumber` to always equal its root. In how many other places might similar expectations be embedded? This approach begins to look like much more trouble than it's worth, so I'll revert to the original approach.

@user1:
I have to say that I've never quite understood how `AlgebraicNumber` works. Maybe it is not used consistently and should be fixed.

@user2:
It seems that an instance of `AlgebraicNumber` is practically always the same as the primitive element except in some tests. Maybe `coeffs` could be deprecated (and the tests removed).

@author:
The possibility did cross my mind, but I think there are strong reasons for keeping it.

* A computer algebra system should have a class that represents "an element of a number field".
Such elements are usually thought of as polynomials in the primitive element, and the class should
provide things like printing support for displaying the element in this form. `AlgebraicNumber`
currently provides all this, but it will cease to if we give up `coeffs`.

* The `ANP` class does things the same way -- i.e. it represents the minimal polynomial of the
primitive element, plus the coeffs of the given element -- but this is on the `polys` side of things.
It is not an `Expr`.

* An `AlgebraicNumber` with `coeffs` is needed as the return value of the public `to_number_field()`
function.

* If we give up `coeffs`, then we just have a `minpoly` and a `root`, in which case we're essentially
the same thing as a `CRootOf`. (Which, to be honest, is a fine representation of an algebraic number for
certain applications -- but then it would be redundant.)

I guess those are the main reasons that I have considered.

Also, I'm still hoping to implement more things from Cohen. So, you're probably right that right now the library mostly treats these things as their primitive elements, but I'm hoping to add some cool stuff which will make an "element of a number field" class of this kind more important!

@user1:
I think that `AlgebraicNumber` should just be fixed to work properly.

@user1:
> But now, due to the funny way AlgebraicNumber.__new__() treats its first arg, when that arg is an AlgebraicNumber,

There `expr.coeffs` should be checked as well. That should be combined with the newly provided coeffs using e.g. `dup_compose`

@author:
But isn't it a backwards compatibility issue? Otherwise I'd be all for changing the behavior here. (I think the special case should just be eliminated.)

@user1:
Can you give an example of what would be incompatible if this was just fixed?

@author:
I guess the real question I was raising was where we draw the line between breaking change and bugfix. But it dawns on me now that this entire special case of `AlgebraicNumber.__new__(expr, coeffs, ...)` where `expr` is an `AlgebraicNumber` could be regarded as a part of the same bug this PR sets out to fix; namely, places in the code where the `coeffs` of an `AlgebraicNumber` were forgotten/ignored.

So, I think the composing of polynomials here is clever, it brings the coeffs back into play, and it means that on unpickling we recover a number which is the same as a complex number, all of which is good. But unfortunately the internal structure still may differ. Example:

Adding this in `.__new__()`:

```python
        if expr.is_AlgebraicNumber:
            from sympy.polys.densetools import dup_compose
            c = dup_compose(rep.rep, expr.rep.rep, dom)
            rep = DMP.from_list(c, 0, dom)
```

we now get:

```python
>>> import pickle
>>> a = to_number_field(sqrt(2), sqrt(2) + sqrt(3))
>>> m = a.minpoly_of_elt()
>>> b = AlgebraicNumber((m, a))
>>> c = pickle.loads(pickle.dumps(b))
>>> for e in [b, c]:
>>>     print('--------')
>>>     print('As a complex number:         ', e.as_expr())
>>>     print('Our coeffs:                  ', e.rep)
>>>     print('Minpoly of primitive element:', e.minpoly)
>>>     print('Primitive element:           ', e.root)
>>> print('--------')
>>> print('Equal after pickling:', b == c)
--------
As a complex number:          sqrt(2)
Our coeffs:                   DMP([mpq(1,1), mpq(0,1)], QQ, None)
Minpoly of primitive element: PurePoly(_x**2 - 2, _x, domain='QQ')
Primitive element:            sqrt(2)
--------
As a complex number:          sqrt(2)
Our coeffs:                   DMP([mpq(1,2), mpq(0,1), mpq(-9,2), mpq(0,1)], QQ, None)
Minpoly of primitive element: PurePoly(_x**4 - 10*_x**2 + 1, _x, domain='QQ')
Primitive element:            sqrt(2) + sqrt(3)
--------
Equal after pickling: False
```

I'm thinking the fix could instead be as simple as just deleting the special case,

https://github.com/sympy/sympy/blob/52f606a503cea5e9588de14150ccb9f7f9ed4752/sympy/core/numbers.py#L2541-L2542

altogether. Then we do test as equal after pickling:

```python
--------
As a complex number:          -9*sqrt(2) + sqrt(3)/2 + sqrt(2) + sqrt(3)**3/2
Our coeffs:                   DMP([mpq(1,1), mpq(0,1)], QQ, None)
Minpoly of primitive element: PurePoly(_x**2 - 2, _x, domain='QQ')
Primitive element:            -9*sqrt(2) + sqrt(3)/2 + sqrt(2) + sqrt(3)**3/2
--------
As a complex number:          -9*sqrt(2) + sqrt(3)/2 + sqrt(2) + sqrt(3)**3/2
Our coeffs:                   DMP([mpq(1,1), mpq(0,1)], QQ, None)
Minpoly of primitive element: PurePoly(_x**2 - 2, _x, domain='QQ')
Primitive element:            -9*sqrt(2) + sqrt(3)/2 + sqrt(2) + sqrt(3)**3/2
--------
Equal after pickling: True

```

Note that the printing under `.as_expr()` is also changed in this case. This is because we are now succeeding in restoring `self.root` as an instance of `AlgebraicNumber`, which stays atomic in simplifications.

I think this is a good thing. If the old behavior is desired, then `b` can just be constructed a little differently:

```python
>>> b = AlgebraicNumber((m, a.to_root()))
>>> print(b.as_expr())
sqrt(2)
```

Here I'm using the `to_root()` method which I've now implemented, and have just pushed. In this case `a` is transformed into a `Pow`, which does participate in expression simplifications.

It's complex. I want to continue thinking about it a bit more...

@user1:
Okay, I see the problem now. Your examples use lots of complicated methods but it's actually just this:
```python
In [15]: a = AlgebraicNumber(sqrt(2), [1, 2])

In [16]: a
Out[16]: √2 + 2

In [17]: AlgebraicNumber(a)
Out[17]: √2
```
That's a clear bug to be fixed. Does the fix need to be much more complicated than this:
```diff
diff --git a/sympy/core/numbers.py b/sympy/core/numbers.py
index e87b039..93c07c8 100644
--- a/sympy/core/numbers.py
+++ b/sympy/core/numbers.py
@@ -2536,7 +2536,7 @@ 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, coeffs = expr.minpoly, expr.root, expr.coeffs()
         else:
             minpoly, root = minimal_polynomial(
                 expr, args.get('gen'), polys=True), expr
```
Could also check if coeffs was provided as well as an AlgebraicNumber and if so use `dup_compose`.

@author:
Okay, this has the virtue of sticking closer to existing behavior, so I think it's good.

But while we're making `AlgebraicNumber(a)` behave more like a "copy constructor" when `a` is an `AlgebraicNumber`, we should probably also copy the `self.alias` from `a`, allowing this to be overridden by any `alias` that may have been passed.

Also, we need to modify the case where `expr` is an ordered pair `(m, r)`. If `r` is an `AlgebraicNumber` we have to convert it with `r = r.to_root()`. This (a) repairs pickling, and (b) makes this case consistent with the others.

I've pushed some code that implements this, and adapts the docs.

@author:
Scratch the bit about `r = r.to_root()`. I tried it, but reverted it since it broke some things.

Anyway this was intended to address a separate issue -- pickling of `AlgebraicNumber` and `AlgebraicField` -- and I'd like to handle this in a separate PR.