Changeset 7504:5908f80abccc

Timestamp:

12/01/07 21:17:45 (5 years ago)

Author:

Carl Witty <cwitty@…>

Branch:

default

Message:

Add more cases for polynomial .roots() (most importantly, ring=QQbar)

Location:

sage/rings

Files:

• polynomial/complex_roots.py (modified) (6 diffs)
• polynomial/polynomial_element.pyx (modified) (8 diffs)
• qqbar.py (modified) (2 diffs)

sage/rings/polynomial/complex_roots.py

@@ -25 +25 @@
 from copy import copy
 
+from sage.rings.real_mpfi import RealIntervalField
 from sage.rings.complex_field import ComplexField
 from sage.rings.complex_interval_field import ComplexIntervalField
+from sage.rings.qqbar import AA, QQbar
 
 def refine_root(ip, ipd, irt, fld):
@@ -111 +113 @@
 
         nirt = center - val / slope
-        if nirt in irt and nirt.diameter() >= irt.diameter() >> 1:
+        if nirt in irt and nirt.diameter() >= irt.diameter() >> 3:
             # If the new diameter is much less than the original diameter,
             # then we have not yet converged.  (Perhaps we were asked
@@ -122 +124 @@
         else:
             irt = irt.union(nirt)
+            # If we don't find a root after a while, try (approximately)
+            # tripling the size of the region.
+            if i >= 6:
+                rd = irt.real().absolute_diameter()
+                id = irt.imag().absolute_diameter()
+                md = max(rd, id)
+                md_intv = RealIntervalField(rd.prec())(-md, md)
+                md_cintv = ComplexIntervalField(rd.prec())(md_intv, md_intv)
+                irt = irt + md_cintv
 
         if not smashed_real and irt.real().contains_zero():
             irt = irt.parent()(0, irt.imag())
+            smashed_real = True
         if not smashed_imag and irt.imag().contains_zero():
             irt = irt.parent()(irt.real(), 0)
+            smashed_imag = True
 
     return None
@@ -244 +257 @@
     has a repeated root, then the algorithm will loop forever!)
 
+    You can specify retval='interval' (the default) to get roots as
+    complex intervals.  The other options are retval='algebraic' to
+    get elements of QQbar, or retval='algebraic_real' to get only
+    the real roots, and to get them as elements of AA.
+
     EXAMPLES:
         sage: from sage.rings.polynomial.complex_roots import complex_roots
@@ -257 +275 @@
     This polynomial actually has all-real coefficients, and is very, very
     close to (x-1)^5:
-        sage: [RR(QQ(x)) for x in list(p - (x-1)^5)]
+        sage: [RR(QQ(a)) for a in list(p - (x-1)^5)]
         [3.87259191484932e-121, -3.87259191484932e-121]
         sage: rts = complex_roots(p)
         sage: [ComplexIntervalField(10)(rt[0] - 1) for rt in rts]
         [0, [7.8886e-31 .. 7.8887e-31]*I, [7.8886e-31 .. 7.8887e-31], [-7.8887e-31 .. -7.8886e-31]*I, [-7.8887e-31 .. -7.8886e-31]]
+
+    We can get roots either as intervals, or as elements of QQbar or AA.
+        sage: p = (x^2 + x - 1)
+        sage: p = p * p(x*im)
+        sage: p
+        (-1)*x^4 + (im - 1)*x^3 + im*x^2 + (-im - 1)*x + 1
+        sage: rts = complex_roots(p); type(rts[0][0]), rts
+        (<type 'sage.rings.complex_interval.ComplexIntervalFieldElement'>, [([-0.000000000000000085405681994008930 .. 0.00000000000000018852810564551411] + [1.6180339887498944 .. 1.6180339887498952]*I, 1), ([-1.6180339887498952 .. -1.6180339887498944] + [-0.00000000000000018852810564551411 .. 0.000000000000000085405681994008930]*I, 1), ([0.61803398874989468 .. 0.61803398874989502] + [-0.00000000000000011497311295586537 .. 0.000000000000000052767663576832497]*I, 1), ([-0.000000000000000054022849329007172 .. 0.00000000000000011371792720369062] - [0.61803398874989468 .. 0.61803398874989502]*I, 1)])
+        sage: rts = complex_roots(p, retval='algebraic'); type(rts[0][0]), rts
+        (<class 'sage.rings.qqbar.AlgebraicNumber'>, [([-3.7249532943328317e-20 .. 3.7095114777437122e-20] + [1.6180339887498946 .. 1.6180339887498950]*I, 1), ([-1.6180339887498950 .. -1.6180339887498946] + [-3.7095114777437122e-20 .. 3.7249532943328317e-20]*I, 1), ([0.61803398874989479 .. 0.61803398874989491] + [-1.8145161579843782e-20 .. 3.1750950176743468e-20]*I, 1), ([-2.5399597786669076e-20 .. 2.0313092838036301e-20] - [0.61803398874989479 .. 0.61803398874989491]*I, 1)])
+        sage: rts = complex_roots(p, retval='algebraic_real'); type(rts[0][0]), rts
+        (<class 'sage.rings.qqbar.AlgebraicReal'>, [([-1.6180339887498950 .. -1.6180339887498946], 1), ([0.61803398874989479 .. 0.61803398874989491], 1)])
     """
 
@@ -284 +314 @@
                 ok = False
                 break
+            if retval != 'interval':
+                factor = QQbar.common_polynomial(factor)
             for irt in irts:
-                all_rts.append((irt, exp))
-
-        if ok and intervals_disjoint([rt for (rt, mult) in all_rts]):
-            return all_rts
+                all_rts.append((irt, factor, exp))
+
+        if ok and intervals_disjoint([rt for (rt, fac, mult) in all_rts]):
+            if retval == 'interval':
+                return [(rt, mult) for (rt, fac, mult) in all_rts]
+            elif retval == 'algebraic':
+                return [(QQbar.polynomial_root(fac, rt), mult) for (rt, fac, mult) in all_rts]
+            elif retval == 'algebraic_real':
+                rts = []
+                for (rt, fac, mult) in all_rts:
+                    qqbar_rt = QQbar.polynomial_root(fac, rt)
+                    if qqbar_rt.imag().is_zero():
+                        rts.append((AA(qqbar_rt), mult))
+                return rts
 
         prec = prec * 2

sage/rings/polynomial/polynomial_element.pyx

@@ -74 +74 @@
 
 cdef object is_AlgebraicRealField
+cdef object is_AlgebraicField
+cdef object is_AlgebraicField_common
 cdef object NumberField_quadratic
 cdef object is_ComplexIntervalField
@@ -80 +82 @@
     # A hack to avoid circular imports.
     global is_AlgebraicRealField
+    global is_AlgebraicField
+    global is_AlgebraicField_common
     global NumberField_quadratic
     global is_ComplexIntervalField
@@ -88 +92 @@
     import sage.rings.qqbar
     is_AlgebraicRealField = sage.rings.qqbar.is_AlgebraicRealField
+    is_AlgebraicField = sage.rings.qqbar.is_AlgebraicField
+    is_AlgebraicField_common = sage.rings.qqbar.is_AlgebraicField_common
     import sage.rings.number_field.number_field
     NumberField_quadratic = sage.rings.number_field.number_field.NumberField_quadratic
@@ -2507 +2513 @@
             sage: p.roots(ring=ComplexIntervalField(200))
             [([1.1673039782614186842560458998548421807205603715254890391400816 .. 1.1673039782614186842560458998548421807205603715254890391400829], 1), ([0.18123244446987538390180023778112063996871646618462304743773153 .. 0.18123244446987538390180023778112063996871646618462304743773341] + [1.0839541013177106684303444929807665742736402431551156543011306 .. 1.0839541013177106684303444929807665742736402431551156543011344]*I, 1), ([0.18123244446987538390180023778112063996871646618462304743773153 .. 0.18123244446987538390180023778112063996871646618462304743773341] - [1.0839541013177106684303444929807665742736402431551156543011306 .. 1.0839541013177106684303444929807665742736402431551156543011344]*I, 1), ([-0.76488443360058472602982318770854173032899665194736756700777454 .. -0.76488443360058472602982318770854173032899665194736756700777...] + [0.35247154603172624931794709140258105439420648082424733283769... .. 0.35247154603172624931794709140258105439420648082424733283769...]*I, 1), ([-0.76488443360058472602982318770854173032899665194736756700777454 .. -0.76488443360058472602982318770854173032899665194736756700777204] - [0.35247154603172624931794709140258105439420648082424733283769122 .. 0.35247154603172624931794709140258105439420648082424733283769341]*I, 1)]
+            sage: rts = p.roots(ring=QQbar); rts
+            [([1.1673039782614185 .. 1.1673039782614188], 1), ([0.18123244446987538 .. 0.18123244446987541] + [1.0839541013177105 .. 1.0839541013177108]*I, 1), ([0.18123244446987538 .. 0.18123244446987541] - [1.0839541013177105 .. 1.0839541013177108]*I, 1), ([-0.76488443360058478 .. -0.76488443360058466] + [0.35247154603172620 .. 0.35247154603172626]*I, 1), ([-0.76488443360058478 .. -0.76488443360058466] - [0.35247154603172620 .. 0.35247154603172626]*I, 1)]
+            sage: p.roots(ring=AA)
+            [([1.1673039782614185 .. 1.1673039782614188], 1)]
+            sage: p = (x - rts[1][0])^2 * (x^2 + x + 1)
+            sage: p.roots(ring=QQbar)
+            [([-0.50000000000000012 .. -0.49999999999999994] + [0.86602540378443859 .. 0.86602540378443871]*I, 1), ([-0.50000000000000000 .. -0.50000000000000000] - [0.86602540378443859 .. 0.86602540378443871]*I, 1), ([0.18123244446987538 .. 0.18123244446987541] + [1.0839541013177105 .. 1.0839541013177108]*I, 2)]
+            sage: p.roots(ring=CIF)
+            [([-0.50000000000000012 .. -0.49999999999999994] + [0.86602540378443859 .. 0.86602540378443871]*I, 1), ([-0.50000000000000000 .. -0.50000000000000000] - [0.86602540378443859 .. 0.86602540378443871]*I, 1), ([0.18123244446987538 .. 0.18123244446987541] + [1.0839541013177105 .. 1.0839541013177108]*I, 2)]
 
         Note that coefficients in a number field with defining polynomial
@@ -2551 +2566 @@
 
         Note that we can find the roots of a polynomial with
-        algebraic real coefficients:
+        algebraic coefficients:
         
             sage: rt2 = sqrt(AA(2))
@@ -2591 +2606 @@
         this method gives.)
 
-        If L is CIF, and K is ZZ, QQ, AA, or the Gaussian rationals,
-        then the root isolation algorithm
+        If L is QQbar or CIF, and K is ZZ, QQ, AA, QQbar, or the
+        Gaussian rationals, then the root isolation algorithm
         sage.rings.polynomial.complex_roots.complex_roots() is used.
-        (You can call complex_roots() directly to get more control than
-        this method gives.)
+        (You can call complex_roots() directly to get more control
+        than this method gives.)
+
+        If L is AA and K is QQbar or the Gaussian rationals, then
+        complex_roots() is used (as above) to find roots in QQbar,
+        then these roots are filtered to select only the real roots.
 
         If L is floating-point and K is not, then we attempt to change
@@ -2714 +2733 @@
                 return rts
 
-        if L != K or is_AlgebraicRealField(L):
+        if L != K or is_AlgebraicField_common(L):
             # So far, the only "special" implementations are for real
             # and complex root isolation.
@@ -2748 +2767 @@
 
             if (is_IntegerRing(K) or is_RationalField(K)
-                or is_AlgebraicRealField(K) or input_gaussian) and \
-                (is_ComplexIntervalField(L)):
+                or is_AlgebraicField_common(K) or input_gaussian) and \
+                (is_ComplexIntervalField(L) or is_AlgebraicField_common(L)):
 
                 from sage.rings.polynomial.complex_roots import complex_roots
 
-                rts = complex_roots(self, min_prec=L.prec())
+                if is_ComplexIntervalField(L):
+                    rts = complex_roots(self, min_prec=L.prec())
+                elif is_AlgebraicField(L):
+                    rts = complex_roots(self, retval='algebraic')
+                else:
+                    rts = complex_roots(self, retval='algebraic_real')
 
                 if multiplicities:

sage/rings/qqbar.py

@@ -645 +645 @@
 
 QQbar = AlgebraicField()
+
+def is_AlgebraicField_common(F):
+    return isinstance(F, AlgebraicField_common)
 
 def prec_seq():
@@ -1628 +1631 @@
                  isinstance(x, NumberFieldElement_quadratic) and \
                  list(x.parent().polynomial()) == [1, 0, 1]:
-            self._descr = ANExtensionElement(QQbar_I_generator, QQbar_I_nf(x))
+            self._descr = ANExtensionElement(QQbar_I_generator, QQbar_I_nf(x.list()))
         else:
             raise TypeError, "Illegal initializer for algebraic number"