Repository: sympy/sympy

Issue:
Bug in `_minpoly_groebner()`
## Example

```python
>>> from sympy import AlgebraicNumber, minimal_polynomial, sqrt
>>> alpha = AlgebraicNumber(sqrt(2))
>>> minimal_polynomial(alpha**3, x)
Traceback (most recent call last):
...
sympy.polys.polyerrors.PolynomialError: _a1**(1/3) contains an element of the set of generators.

```


## Background

`minimal_polynomial()` has two methods that it can use internally: `_minpoly_compose()` and `_minpoly_groebner()`.
The former uses resultants, the latter computes a Groebner basis.

Some time ago (commit c4c9d96f0, May 2013) it was decided that only `_minpoly_groebner()` would be used on expressions
involving one or more `AlgebraicNumber`s. (No comment as to why. Maybe it's faster? Either method _could_ be used.)

Aside from a couple of optimizations, what `_minpoly_groebner()` does is:

* Transform the given algebraic expression `e` (whose minimal polynomial is desired) into a multivariate polynomial `f`.

* In order to do this, we must replace any and all root extractions or `AlgebraicNumber`s occurring within `e` with new
  variables `a_i`. For each new variable `a_i`, we record a new equation that defines it. For example, if `a_0**(1/3)` occurs, we will replace it with `a_1`, and record the new defining equation `a_1**3 - a_0` (with RHS implicitly 0).

* Finally we record one more equation, `x - f`.

* Then compute a Groebner basis for the set of all equations, with `x` being the final free symbol.

* Then the last element of the Groebner basis is a polynomial `g` having `e` as a root.

* Factor `g` and choose the right irreducible factor.


## The Bug

The bug is in the scanning process, where we introduce new variables `a_i`, in particular where we handle powers `Pow(b, e)`.

First of all, `e` must be rational, or we give up, which is the correct behavior.

The code handles properly the case where `e` is _not_ an integer. This is the root extraction case, where e.g. we need to
replace `a_0**(1/3)` with `a_1`.

The bug is that the case where `e` _is_ an integer is treated in the same way. In other words, if `a_0**3` occurs -- which should not actually necessitate any new variable at all -- we introduce `a_1` with defining "polynomial" `a_1**(1/3) - a_0`. This of course is not a polynomial, so the attempt to compute a Groebner basis crashes.

Until a very recent commit (88a18df58a of Nov 10) the integer power case had no unit test coverage, so it's not surprising this issue was missed.

Discussion:

Patch:
diff --git a/sympy/polys/numberfields/minpoly.py b/sympy/polys/numberfields/minpoly.py
index 26b031da0d..797a2af205 100644
--- a/sympy/polys/numberfields/minpoly.py
+++ b/sympy/polys/numberfields/minpoly.py
@@ -742,6 +742,33 @@ def update_mapping(ex, exp, base=None):
         return a
 
     def bottom_up_scan(ex):
+        """
+        Transform a given algebraic expression *ex* into a multivariate
+        polynomial, by introducing fresh variables with defining equations.
+
+        Explanation
+        ===========
+
+        The critical elements of the algebraic expression *ex* are root
+        extractions, instances of :py:class:`~.AlgebraicNumber`, and negative
+        powers.
+
+        When we encounter a root extraction or an :py:class:`~.AlgebraicNumber`
+        we replace this expression with a fresh variable ``a_i``, and record
+        the defining polynomial for ``a_i``. For example, if ``a_0**(1/3)``
+        occurs, we will replace it with ``a_1``, and record the new defining
+        polynomial ``a_1**3 - a_0``.
+
+        When we encounter a negative power we transform it into a positive
+        power by algebraically inverting the base. This means computing the
+        minimal polynomial in ``x`` for the base, inverting ``x`` modulo this
+        poly (which generates a new polynomial) and then substituting the
+        original base expression for ``x`` in this last polynomial.
+
+        We return the transformed expression, and we record the defining
+        equations for new symbols using the ``update_mapping()`` function.
+
+        """
         if ex.is_Atom:
             if ex is S.ImaginaryUnit:
                 if ex not in mapping:
@@ -774,7 +801,10 @@ def bottom_up_scan(ex):
                 expr = base**exp
 
                 if expr not in mapping:
-                    return update_mapping(expr, 1/exp, -base)
+                    if exp.is_Integer:
+                        return expr.expand()
+                    else:
+                        return update_mapping(expr, 1 / exp, -base)
                 else:
                     return symbols[expr]
         elif ex.is_AlgebraicNumber:


Test patch:
diff --git a/sympy/polys/numberfields/tests/test_minpoly.py b/sympy/polys/numberfields/tests/test_minpoly.py
index a4fad8b8c8..ee202b9146 100644
--- a/sympy/polys/numberfields/tests/test_minpoly.py
+++ b/sympy/polys/numberfields/tests/test_minpoly.py
@@ -422,6 +422,11 @@ def test_issue_20163():
         I/(x + I)/6 - I/(x - I)/6
 
 
+def test_issue_22559():
+    alpha = AlgebraicNumber(sqrt(2))
+    assert minimal_polynomial(alpha**3, x) == x**2 - 8
+
+
 def test_separate_sq_not_impl():
     raises(NotImplementedError, lambda: _separate_sq(x**(S(1)/3) + x))