Implement minpoly() for approximate values

src/sage/arith/misc.py

@@ -185,6 +185,35 @@ def algdep(z, degree, known_bits=None, use_bits=None, known_digits=None,
         sage: from gmpy2 import mpz, mpfr
         sage: algdep(mpfr(1.888888888888888), mpz(1))                                   # needs sage.libs.pari
         9*x - 17
+
+    Test with interval arithmetic::
+
+        sage: algdep(RealIntervalField(100)(sqrt(2)), 2)
+        x^2 - 2
+        sage: algdep(RealIntervalField(100)(sqrt(2) + 2^-99), 2)
+        Traceback (most recent call last):
+        ...
+        ValueError: cannot find a polynomial
+        sage: algdep(RealIntervalField(100)((x^5 + x^3 - 2*x^2 + x + 1).roots(x, ring=QQbar)[0][0]), 5)
+        x^5 + x^3 - 2*x^2 + x + 1
+        sage: algdep(ComplexIntervalField(100)(sqrt(-1)), 2)
+        x^2 + 1
+        sage: algdep(RealBallField(100)(sqrt(2)), 2)
+        x^2 - 2
+        sage: algdep(RealBallField(100)(sqrt(2) + 2^-96), 2)
+        Traceback (most recent call last):
+        ...
+        ValueError: cannot find a polynomial
+        sage: algdep(ComplexBallField(100)(sqrt(-1)), 2)
+        x^2 + 1
+
+    The result might be nonsensical if ``degree`` is too high,
+    or if there is no satisfying polynomial::
+
+        sage: algdep(RR(1.01), 8)
+        x - 1
+
+    In the case above, the found polynomial may be ``(x + 1) * (x - 1)^7``.
     """
     if proof and not height_bound:
         raise ValueError("height_bound must be given for proof=True")
@@ -206,7 +235,17 @@ def algdep(z, degree, known_bits=None, use_bits=None, known_digits=None,
             return None
         return z.denominator() * x - z.numerator()
 
-    if isinstance(z.parent(), (RealField, ComplexField)):
+    from sage.rings.complex_arb import ComplexBallField
+    from sage.rings.complex_interval_field import ComplexIntervalField, ComplexIntervalField_class
+    from sage.rings.real_arb import RealBallField
+    from sage.rings.real_mpfi import RealIntervalField, RealIntervalField_class
+
+    is_interval = isinstance(parent(z), (RealBallField, ComplexBallField, RealIntervalField_class, ComplexIntervalField_class))
+    if is_interval:
+        original_z = z
+        z = z.center()
+
+    if isinstance(parent(z), (RealField, ComplexField)):
 
         log2_10 = math.log(10, 2)
 
@@ -220,7 +259,7 @@ def algdep(z, degree, known_bits=None, use_bits=None, known_digits=None,
         if use_bits is not None:
             prec = int(use_bits)
 
-        is_complex = isinstance(z.parent(), ComplexField)
+        is_complex = isinstance(parent(z), ComplexField)
         n = degree + 1
         from sage.matrix.constructor import matrix
         M = matrix(ZZ, n, n + 1 + int(is_complex))
@@ -269,8 +308,28 @@ def norm(v):
         y = pari(z)
         f = y.algdep(degree)
 
+    f = R(f)
+    if is_interval:
+        P = parent(original_z)
+        is_real = False
+        if isinstance(P, RealBallField):
+            P1 = RealBallField(P.precision() + 10)
+            is_real = True
+        elif isinstance(P, ComplexBallField):
+            P1 = ComplexBallField(P.precision() + 10)
+        elif isinstance(P, RealIntervalField_class):
+            P1 = RealIntervalField(P.precision() + 10)
+            is_real = True
+        else:
+            assert isinstance(P, ComplexIntervalField_class)
+            P1 = ComplexIntervalField(P.precision() + 10)
+
+        from sage.rings.qqbar import QQbar, AA
+        if not any(r in P1(original_z) for r in f.roots(AA if is_real else QQbar, multiplicities=False)):
+            raise ValueError("cannot find a polynomial")
+
     # f might be reducible. Find the best fitting irreducible factor
-    factors = [p for p, e in R(f).factor()]
+    factors = [p for p, e in f.factor()]
     return min(factors, key=lambda f: abs(f(z)))
 
 

src/sage/calculus/calculus.py

@@ -418,13 +418,20 @@
     2
 """
 
+import math
 import re
 from sage.arith.misc import algdep
+from sage.rings.complex_arb import ComplexBallField
+from sage.rings.complex_interval_field import ComplexIntervalField_class
+from sage.rings.complex_mpfr import ComplexField_class
 from sage.rings.integer import Integer
 from sage.rings.rational_field import QQ
+from sage.rings.real_arb import RealBallField
+from sage.rings.real_mpfi import RealIntervalField_class
 from sage.rings.real_double import RealDoubleElement
-from sage.rings.real_mpfr import RR, create_RealNumber
+from sage.rings.real_mpfr import RR, create_RealNumber, RealField_class
 from sage.rings.cc import CC
+from sage.structure.element import parent
 
 from sage.misc.latex import latex
 from sage.misc.parser import Parser, LookupNameMaker
@@ -1106,7 +1113,45 @@
        Of course, failure to produce a minimal polynomial does not
        necessarily indicate that this number is transcendental.
     """
+    if isinstance(ex, Expression):
+        try:
+            ex = ex.pyobject()
+        except TypeError:
+            pass
+
     if algorithm is None or algorithm.startswith('numeric'):
+        def estimated_bit_length(f):
+            return sum(x.numer().bit_length() + x.denom().bit_length() + 1 for x in f)
+
+        if parent(ex) is float or isinstance(parent(ex),
+                                             (RealBallField, ComplexBallField, RealIntervalField_class,
+                                              ComplexIntervalField_class, RealField_class, ComplexField_class)):
+            best_f = None
+            cur_degree = degree or 2  # ``algdep(CC(I), 1)`` just return ``x``
+            prec = parent(ex).prec()
+            while True:
+                if best_f is not None and cur_degree > estimated_bit_length(best_f) / 4:
+                    break
+                if cur_degree > prec / 6:
+                    # experimentally algdep(RealField(500)(2^(1/100)), 100) returns the wrong answer
+                    break
+                try:
+                    f = algdep(ex, cur_degree)
+                except ValueError:
+                    pass
+                else:
+                    f = QQ[var](f)
+                    if best_f is None or estimated_bit_length(f) < estimated_bit_length(best_f):
+                        best_f = f
+                    if degree:
+                        # user-specified, do not try any other degree value
+                        break
+                cur_degree = math.ceil(cur_degree*1.2)
+
+            if best_f is not None:
+                return best_f
+            raise ValueError("Could not find minimal polynomial")
+
         bits_list = [bits] if bits else [100,200,500,1000]
         degree_list = [degree] if degree else [2,4,8,12,24]
 

src/sage/rings/complex_arb.pyx

@@ -2037,6 +2037,8 @@ cdef class ComplexBall(RingElement):
                 max(prec(self), re.prec(), im.prec()))
         return field(re, im)
 
+    center = mid
+
     def squash(self):
         """
         Return an exact ball with the same midpoint as this ball.
@@ -5168,5 +5170,31 @@ cdef class ComplexBall(RingElement):
         if _do_sig(prec(self)): sig_off()
         return result
 
+    def minpoly(self, *args, **kwargs):
+        """
+        Finds some polynomial with integer coefficients that has a root inside the interval.
+
+        TESTS::
+
+            sage: CIF(1).minpoly()
+            x - 1
+            sage: CIF(1/2).minpoly()
+            2*x - 1
+            sage: CIF(sqrt(-1)).minpoly()
+            x^2 + 1
+            sage: CIF(sin(1)).minpoly()
+            Traceback (most recent call last):
+            ...
+            ValueError: Could not find minimal polynomial
+            sage: minpoly(CIF((-1)^(1/50)))  # insufficient precision
+            Traceback (most recent call last):
+            ...
+            ValueError: Could not find minimal polynomial
+            sage: minpoly(ComplexIntervalField(400)((-1)^(1/50)))
+            x^40 - x^30 + x^20 - x^10 + 1
+        """
+        from sage.calculus.calculus import minpoly
+        return minpoly(self, *args, **kwargs)
+
 
 CBF = ComplexBallField()

src/sage/rings/complex_interval.pyx

@@ -2124,6 +2124,32 @@ cdef class ComplexIntervalFieldElement(FieldElement):
         return ComplexBallField(self.prec())(self).zeta(a).\
             _complex_mpfi_(self._parent)
 
+    def minpoly(self, *args, **kwargs):
+        """
+        Finds some polynomial with integer coefficients that has a root inside the interval.
+
+        TESTS::
+
+            sage: CIF(1).minpoly()
+            x - 1
+            sage: CIF(1/2).minpoly()
+            2*x - 1
+            sage: CIF(sqrt(-1)).minpoly()
+            x^2 + 1
+            sage: CIF(sin(1)).minpoly()
+            Traceback (most recent call last):
+            ...
+            ValueError: Could not find minimal polynomial
+            sage: minpoly(CIF((-1)^(1/50)))  # insufficient precision
+            Traceback (most recent call last):
+            ...
+            ValueError: Could not find minimal polynomial
+            sage: minpoly(ComplexIntervalField(400)((-1)^(1/50)))
+            x^40 - x^30 + x^20 - x^10 + 1
+        """
+        from sage.calculus.calculus import minpoly
+        return minpoly(self, *args, **kwargs)
+
 cdef _negate_interval(mpfr_ptr xmin, mpfr_ptr xmax):
     """
     Negate interval with endpoints xmin and xmax in place.

src/sage/rings/complex_mpfr.pyx

@@ -3255,7 +3255,7 @@ cdef class ComplexNumber(sage.structure.element.FieldElement):
         Return an irreducible polynomial of degree at most `n` which is
         approximately satisfied by this complex number.
 
-        ALGORITHM: Uses the PARI C-library :pari:`algdep` command.
+        ALGORITHM: See :func:`~sage.arith.misc.algdep`.
 
         INPUT: Type ``algdep?`` at the top level prompt. All additional
         parameters are passed onto the top-level :func:`algdep` command.
@@ -3276,6 +3276,29 @@ cdef class ComplexNumber(sage.structure.element.FieldElement):
     # Alias
     algdep = algebraic_dependency
 
+    def minpoly(self, *args, **kwargs):
+        """
+        Finds some polynomial with integer coefficients that has a root
+        approximately equal to ``self``.
+
+        TESTS::
+
+            sage: CC(1).minpoly()
+            x - 1
+            sage: CC(1/2).minpoly()
+            2*x - 1
+            sage: CC(sqrt(-1)).minpoly()
+            x^2 + 1
+            sage: CC(pi*(-1)^(1/3)).minpoly()  # transcendental, might return arbitrary polynomial or error
+            793891*x^3 + 24615604
+            sage: minpoly(CC((-1)^(1/50)))  # insufficient precision
+            736594*x^2 - 1470281*x + 736594
+            sage: minpoly(ComplexField(400)((-1)^(1/50)))
+            x^40 - x^30 + x^20 - x^10 + 1
+        """
+        from sage.calculus.calculus import minpoly
+        return minpoly(self, *args, **kwargs)
+
 
 def make_ComplexNumber0(fld, mult_order, real, imag):
     """

src/sage/rings/real_arb.pyx

@@ -4046,4 +4046,30 @@ cdef class RealBall(RingElement):
         if _do_sig(prec(self)): sig_off()
         return res
 
+    def minpoly(self, *args, **kwargs):
+        """
+        Finds some polynomial with integer coefficients that has a root inside the interval.
+
+        TESTS::
+
+            sage: RBF(1).minpoly()
+            x - 1
+            sage: RBF(1/2).minpoly()
+            2*x - 1
+            sage: RBF(sqrt(2)).minpoly()
+            x^2 - 2
+            sage: RBF(sin(1)).minpoly()
+            Traceback (most recent call last):
+            ...
+            ValueError: Could not find minimal polynomial
+            sage: minpoly(RBF(2^(1/50)))  # insufficient precision
+            Traceback (most recent call last):
+            ...
+            ValueError: Could not find minimal polynomial
+            sage: minpoly(RealBallField(400)(2^(1/50)))
+            x^50 - 2
+        """
+        from sage.calculus.calculus import minpoly
+        return minpoly(self, *args, **kwargs)
+
 RBF = RealBallField()

src/sage/rings/real_mpfi.pyx

@@ -5185,6 +5185,34 @@ cdef class RealIntervalFieldElement(RingElement):
         return RealBallField(self.precision())(self).zeta(a).\
             _real_mpfi_(self._parent)
 
+    def minpoly(self, *args, **kwargs):
+        """
+        Finds some polynomial with integer coefficients that has a root inside the interval.
+
+        TESTS::
+
+            sage: RIF(1).minpoly()
+            x - 1
+            sage: RIF(1/2).minpoly()
+            2*x - 1
+            sage: RIF(sqrt(2)).minpoly()
+            x^2 - 2
+            sage: RIF(sin(1)).minpoly()
+            Traceback (most recent call last):
+            ...
+            ValueError: Could not find minimal polynomial
+            sage: minpoly(RIF(2^(1/50)))  # insufficient precision
+            Traceback (most recent call last):
+            ...
+            ValueError: Could not find minimal polynomial
+            sage: minpoly(RealIntervalField(200)(2^(1/50)))  # insufficient precision, returns an arbitrary polynomial
+            17*x^27 + 17*x^26 + 6*x^25 + 25*x^24 + 2*x^23 - 28*x^22 - 72*x^21 - 51*x^20 - 12*x^19 - 82*x^18 + 42*x^17 - 35*x^16 + 42*x^15 - 78*x^14 + 34*x^13 + 33*x^12 + 9*x^11 - 7*x^10 + 77*x^9 - 48*x^8 + 102*x^7 + 53*x^6 + 27*x^5 + 16*x^4 - 32*x^3 - 35*x^2 - 31*x + 30
+            sage: minpoly(RealIntervalField(400)(2^(1/50)))
+            x^50 - 2
+        """
+        from sage.calculus.calculus import minpoly
+        return minpoly(self, *args, **kwargs)
+
 
 def _simplest_rational_test_helper(low, high, low_open=False, high_open=False):
     """

src/sage/rings/real_mpfr.pyx

@@ -3959,6 +3959,29 @@ cdef class RealNumber(sage.structure.element.RingElement):
         """
         return gmpy2.GMPy_MPFR_From_mpfr(self.value)
 
+    def minpoly(self, *args, **kwargs):
+        """
+        Finds some polynomial with integer coefficients that has a root
+        approximately equal to ``self``.
+
+        TESTS::
+
+            sage: RR(1).minpoly()
+            x - 1
+            sage: RR(1/2).minpoly()
+            2*x - 1
+            sage: RR(sqrt(2)).minpoly()
+            x^2 - 2
+            sage: RR(pi).minpoly()  # transcendental, might return arbitrary polynomial or error
+            6201*x^2 + 3606*x - 72530
+            sage: minpoly(RR(2^(1/50)))  # insufficient precision
+            x - 1
+            sage: minpoly(RealField(400)(2^(1/50)))
+            x^50 - 2
+        """
+        from sage.calculus.calculus import minpoly
+        return minpoly(self, *args, **kwargs)
+
     ###########################################
     # Comparisons: ==, !=, <, <=, >, >=
     ###########################################