From 8c52b0ddc72fd5ceb02b7c460683d4bd08714790 Mon Sep 17 00:00:00 2001 From: Matthias Goerner <1239022+unhyperbolic@users.noreply.github.com> Date: Wed, 10 Jul 2024 19:04:53 -0700 Subject: [PATCH 1/3] Implementing tight intervals in ComplexIntervalFieldElement.__invert__ - thus for complex interval division. The implementation is inspired by Rokne-Lancaster 1971. Some history: - The first implementation of returning tight intervals was commit 84ba655 (see ticket #19964). - That implementation had bugs (returning too small intervals) and was reverted to fix ticket #37941. - This change is the second implementation. --- src/sage/rings/complex_interval.pyx | 429 +++++++++++++++++++++++++++- 1 file changed, 415 insertions(+), 14 deletions(-) diff --git a/src/sage/rings/complex_interval.pyx b/src/sage/rings/complex_interval.pyx index 5b3064688b0..575829d256c 100644 --- a/src/sage/rings/complex_interval.pyx +++ b/src/sage/rings/complex_interval.pyx @@ -1166,24 +1166,161 @@ cdef class ComplexIntervalFieldElement(FieldElement): - [RL1971]_ """ - cdef ComplexIntervalFieldElement x + cdef ComplexIntervalFieldElement result x = self._new() - cdef mpfi_t t0, t1 - mpfi_init2(t0, self._prec) - mpfi_init2(t1, self._prec) - - mpfi_sqr(t0, self.__re) - mpfi_sqr(t1, self.__im) - - mpfi_add(t0, t0, t1) # now t0 is the norm - mpfi_div(x.__re, self.__re, t0) # x.__re = self.__re/norm + if mpfi_nan_p(self.__re) or mpfi_nan_p(self.__im): + # Return NaN if any input is NaN + return x - mpfi_neg(t1, self.__im) - mpfi_div(x.__im, t1, t0) # x.__im = -self.__im/norm + if mpfi_has_zero(self.__re) and mpfi_has_zero(self.__im): + # Return NaN when dividing by (a complex interval containing) zero. + return x - mpfi_clear(t0) - mpfi_clear(t1) + # The algorithm roughly follows [RL1971]. + # + # However, we deviate from [RL1971] significantly and the + # documentation here is self-contained and hopefully easier + # to follow than [RL1971]. There is a visualization at: + # https://www.shadertoy.com/view/M3VXWG + # + # Note that [RL1971] has several mistakes. For example, when analyzing + # the second case where y1 < 0 and y2 > 0 (and x1 >= 0 from earlier + # assumptions), they have a subcase "y_1 >= x_1". That inequality can + # never be true given the earlier assumptions and should be + # "-y_1 >= x_1". + + # To maximize symmetry considerations, we actually compute the + # circle inversion (z |-> conjugate(1/z)) and conjugate at the end. + + # Let the input be the complex interval [xmin, xmax] + [ymin, ymax] * I. + # + # We say that the complex interval is in standard form if either: + # (I) It is contained within the first quadrant. Furthermore its + # left bottom corner is on or below the north east diagonal. + # That is 0 <= ymin <= xmin. Or: + # (II) It is contained within the first and fourth quadrant crossing + # the (positive) x-Axis. Furthermore, its midpoint is above the + # x-Axis. + # That is 0 < xmin; ymin < 0 < ymax and |ymin| <= |ymax|. + # + # Since we already guarded against the input containing zero, we always + # can and will bring it into standard form by negating or swapping the + # real and imaginary part. + + cdef mpfr_t xmin, xmax, ymin, ymax + mpfr_init2(xmin, self._prec) + mpfr_init2(xmax, self._prec) + mpfr_init2(ymin, self._prec) + mpfr_init2(ymax, self._prec) + + mpfi_get_left(xmin, self.__re) + mpfi_get_right(xmax, self.__re) + mpfi_get_left(ymin, self.__im) + mpfi_get_right(ymax, self.__im) + + # Record what mirror symmetries we applied to undo them later. + # + # We assume that we flip about the coordinate axes before flipping + # about the north east diagonal. + cdef bint negated_x = False + cdef bint negated_y = False + cdef bint swapped_xy = False + + # Record whether we are in case (II). + cdef bint crosses_x_axis = False + + ######################################################################## + # Now bring the input into standard form. + + if mpfr_sgn(ymax) <= 0: + # Interval below (and possibly touching) x-Axis, flip about it. + _negate_interval(ymin, ymax) + negated_y = True + elif mpfr_sgn(ymin) < 0: + # Interval crosses x-Axis + crosses_x_axis = True + # else: Interval is above x-Axis + + # Negating interval for y and swapping do not commute, so + # order of the above and below if-block is important. + + if mpfr_sgn(xmax) <= 0: + # Interval left of (and possibly touching) y-Axis, flip about it. + _negate_interval(xmin, xmax) + negated_x = True + elif mpfr_sgn(xmin) < 0: + # Interval crosses y-Axis, swap to make it cross x-Axis instead. + mpfr_swap(xmin, ymin) + mpfr_swap(xmax, ymax) + swapped_xy = True + crosses_x_axis = True + # else: Interval is right of y-Axis + + if crosses_x_axis: + # Case (II). Ensure standard condition |ymax| >= |ymin| + if mpfr_cmpabs(ymin, ymax) > 0: + _negate_interval(ymin, ymax) + # Note that we negate after potentially swapping. + # Determine what we should have negated before swapping instead. + if swapped_xy: + # negated_x cannot be True already. + negated_x = True + else: + # negated_y cannot be True already. + negated_y = True + else: + # Case (I). Ensure standard condition |ymin| <= |xmin|. + if mpfr_cmp(xmin, ymin) < 0: + mpfr_swap(xmin, ymin) + mpfr_swap(xmax, ymax) + swapped_xy = True + + ######################################################################## + # Apply circle inversion + + # Result will be [amin, amax] + [bmin, bmax] * I. + cdef mpfr_t amin, amax, bmin, bmax + + mpfr_init2(amin, self._prec) + mpfr_init2(amax, self._prec) + mpfr_init2(bmin, self._prec) + mpfr_init2(bmax, self._prec) + + _circle_invert_standard( + amin, amax, bmin, bmax, + xmin, xmax, ymin, ymax, + crosses_x_axis, self._prec) + + ######################################################################## + # Undo symmetries applied above. + + if swapped_xy: + mpfr_swap(amin, bmin) + mpfr_swap(amax, bmax) + + if negated_x: + _negate_interval(amin, amax) + + # We sneak in the promised conjugation by + # inverting negated_y. + if not negated_y: + _negate_interval(bmin, bmax) + + ######################################################################## + # Write out result and free memory + + mpfi_interv_fr(x.__re, amin, amax) + mpfi_interv_fr(x.__im, bmin, bmax) + + mpfr_clear(xmin) + mpfr_clear(xmax) + mpfr_clear(ymin) + mpfr_clear(ymax) + mpfr_clear(amin) + mpfr_clear(amax) + mpfr_clear(bmin) + mpfr_clear(bmax) return x @@ -1988,6 +2125,270 @@ cdef class ComplexIntervalFieldElement(FieldElement): return ComplexBallField(self.prec())(self).zeta(a).\ _complex_mpfi_(self._parent) +cdef _negate_interval(mpfr_ptr xmin, mpfr_ptr xmax): + """ + Negate interval with endpoints xmin and xmax in place. + """ + mpfr_swap(xmin, xmax) + mpfr_neg(xmin, xmin, MPFR_RNDD) + mpfr_neg(xmax, xmax, MPFR_RNDU) + +cdef _inversion_coordinate( + mpfr_ptr result, mpfr_ptr tmp, + mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_denom, mpfr_rnd_t rnd): + """ + Compute the x-coordinate of the image of the point (x,y) under inversion. + That is compute x / (x^2 + y^2). + + rnd_denom determines how to round when computing the denominator. + rnd determines how to round when doing the division. + + This function expects the callee to allocate and pass an mpfr number tmp + that will be used for intermediate computations. This way, an allocated + mpfr number can be used in multiple calculations. + """ + + mpfr_sqr(tmp, x, rnd_denom) + mpfr_sqr(result, y, rnd_denom) + mpfr_add(tmp, tmp, result, rnd_denom) + mpfr_div(result, x, tmp, rnd) + +cdef _inversion_coordinate_pos_down( + mpfr_ptr result, mpfr_ptr tmp, + mpfr_srcptr x, mpfr_srcptr y): + """ + Computes a lower bound for x / (x^2 + y^2) assuming that x is non-negative. + + Note that the result will not be NaN as long as x or y is positive. + """ + + # Let us check that we do not get NaN. + # We need to check mpfr_div since it could return NaN if we divide by zero + # (it currently returns +/-inf, but the documentation says this might be + # subject to change). + # We round up when computing the denominator and thus should get something + # positive if x or y is positive. Even if, say, x is so small that squaring + # produces something smaller than the smallest positive mpfr floating point + # number, we get something non-zero because we round up. + + _inversion_coordinate( + result, tmp, x, y, + MPFR_RNDU, # denominator rounding + MPFR_RNDD) # division rounding + +cdef _inversion_coordinate_pos_up( + mpfr_ptr result, mpfr_ptr tmp, + mpfr_srcptr x, mpfr_srcptr y): + """ + Computes an upper bound for x / (x^2 + y^2) assuming that x is non-negative. + """ + _inversion_coordinate( + result, tmp, x, y, + MPFR_RNDD, # denominator rounding + MPFR_RNDU) # division rounding + +cdef _inversion_coordinate_neg_down( + mpfr_ptr result, mpfr_ptr tmp, + mpfr_srcptr x, mpfr_srcptr y): + """ + Computes a lower bound for x / (x^2 + y^2) assuming that x is non-positive. + """ + _inversion_coordinate( + result, tmp, x, y, + MPFR_RNDD, # denominator rounding + MPFR_RNDD) # division rounding + +cdef _circle_invert_standard( + mpfr_ptr amin, mpfr_ptr amax, mpfr_ptr bmin, mpfr_ptr bmax, + mpfr_srcptr xmin, mpfr_srcptr xmax, mpfr_srcptr ymin, mpfr_srcptr ymax, + bint crosses_x_axis, mpfr_prec_t prec): + """ + Assumes that the input [xmin, xmax] + [ymin, ymax] * I is in standard form + as described above in ComplexIntervalFieldElement.__invert__ with + crosses_x_axis saying whether case (II) applies. + + Computes a complex interval [amin, amax] + [bmin, bmax] * I containing + the image of the input under circle inversion f(z) = conjugate(1/z). + """ + + # Note that, by the maximum principle, it is sufficient to consider the + # image of the boundary of the input rect which will be formed by four + # arcs or line segments. + + # It is useful to do a case analysis by considering whether the input + # crosses the x-Axis, the north east or south east diagonal, respectively. + # + # Given standard form, the input also has to cross the north east + # diagonal and x-Axis if it corsses the south east diagonal. + # + # Thus, we are left with five cases: + # + # NE diagonal x-Axis SE diagonal. + # 1. + # 2. crosses + # 3. crosses + # 4. crosses crosses + # 5. crosses crosses crosses + + # The reader can go to https://www.shadertoy.com/view/M3VXWG to explore + # the different cases visually. + + # Case 1 is the easiest (and the generic case for small intervals). + # + # Consider the images + # f(xmin + ymin * I), ..., f(xmax + ymax * I) + # of the four corners of the input rect under inversion f. + # Now consider the the axis-parallel rectangle R that these images span. + # In general, the image of the input rect might not be contained in R. + # In case 1, however, (and only in case 1) it is and we furthermore know + # which image is mapped to which edge of R. Thus, we have: + # + # amin = Re(f(xmax + ymax * I)) # Image of right top corner + # amax = Re(f(xmin + ymin * I)) # left bottom corner + # bmin = Im(f(xmax + ymin * I)) # right bottom corner + # = Re(f(ymin + xmax * I)) + # bmax = Im(f(xmin + ymax * I)) # left top corner + # = Re(f(ymax + xmin * I)) + # + # Re(f(...)) can be computed with the correct rounding using one of the + # helper functions such as _inversion_coordinate_pos_down. + + # For the other cases, we might need to consider the images of two corners + # and take the minimum to compute, say amin. + + # Furthermore, we also need to consider the image of t + xmin * I which is + # a circle touching the y-Axis at the origin. Depending on which of the + # diagonals and x-Axis the input rect crosses, we need to expand R to + # include the lowest, highest or rightmost point on this circle for the + # correct result. + # + # For example, we have + # amax = 1 / xmin for case 2 and bmax = 1 / (2 * ymin) for case 3. + # + + cdef bint crosses_NE_diagonal = mpfr_cmp(ymax, xmin) > 0 + cdef bint crosses_both_diagonals = False + # Using that input can only cross south east diagonal if it crosses + # x-Axis and north east diagonal. + if crosses_x_axis and crosses_NE_diagonal: + crosses_both_diagonals = mpfr_cmpabs(ymin, xmin) > 0 + + # Some temporary variables + cdef mpfr_t tmp, min2 + + mpfr_init2(tmp, prec) + if crosses_NE_diagonal: + # Temporary variable needed only in some cases. + mpfr_init2(min2, prec) + + ######################################################################## + # Compute amin + + # Use image of right top corner + _inversion_coordinate_pos_down(amin, tmp, xmax, ymax) + + if crosses_NE_diagonal: + # Also use image of left top corner. + + # This is because in this case, the image of the left top corner + # can be left of the image of the right top corner. + + _inversion_coordinate_pos_down(min2, tmp, xmin, ymax) + + # Note that mpfr_min does not return NaN if one (but not the other) + # input is NaN. This could be a problem. Luckily, + # _inversion_coordinate_pos_down never produces NaN. + mpfr_min(amin, amin, min2, MPFR_RNDD) + + # Note that we do not need to consider the images of the left or right + # bottom corner here. Not even in case 5. + # This is because in standard form, the top edge is further aways from + # the x-Axis than the bottom edge. Thus, the images of the left and + # right corner of the top edge is left of those of the bottom edge, + # respectively. + + # Potential optimization: if bmax >= amax, we only need to use the + # image of the left top corner and skip computing the image of the + # right top corner. + + ######################################################################## + # Compute amax + + if crosses_x_axis: + # Use rightmost point on the circle that is the image of + # image of t + xmin * I + mpfr_ui_div(amax, 1, xmin, MPFR_RNDU) + else: + # Use image of left bottom corner + _inversion_coordinate_pos_up(amax, tmp, xmin, ymin) + + ######################################################################## + # Compute bmax + + # In case 5, bmin can reuse bmax (up to sign), so we compute bmax first. + + if crosses_NE_diagonal: + # Use highest point on the circle that is the image of + # t + xmin * I. That is bmax = 1 / (2 * xmin). + if crosses_x_axis: + # Re-use amax = 1 / xmin. + # We can just copy the mantissa and decrease the exponent by 1. + mpfr_div_2ui(bmax, amax, 1, MPFR_RNDU) + else: + # Compute bmax from scratch. + mpfr_ui_div(bmax, 1, xmin, MPFR_RNDU) + mpfr_div_2ui(bmax, bmax, 1, MPFR_RNDU) + else: + # Use image of of left top corner + _inversion_coordinate_pos_up(bmax, tmp, ymax, xmin) + + ######################################################################## + # Compute bmin + + # bmin is probably the hardest to compute. + + cdef bint right_edge_crosses_NE_diagonal + + if crosses_x_axis: + # ymin and thus bmin will be negative. + if crosses_both_diagonals: + # We are in case 5. + # + # Use lowest point on the circle that is the image of + # t + xmin * I. That is bmin = -1 / (2 * xmin). + # + # We can reuse bmax = 1 / (2 * xmin). + mpfr_neg(bmin, bmax, MPFR_RNDD) + else: + # Use image of left bottom corner. + _inversion_coordinate_neg_down(bmin, tmp, ymin, xmin) + else: + # ymin and thus bmin will be non-negative. + + # Use image of right bottom corner. + _inversion_coordinate_pos_down(bmin, tmp, ymin, xmax) + + if crosses_NE_diagonal: + right_edge_crosses_NE_diagonal = mpfr_cmp(ymax, xmax) > 0 + if right_edge_crosses_NE_diagonal: + # Also use image of right top corner. + # + # This is similar to the computation of amin which also + # considered a second corner when crosses_NE_diagonal. + # + # In particular, the same comment about NaN applies. + # + _inversion_coordinate_pos_down(min2, tmp, ymax, xmax) + + # See comment about NaN above. + mpfr_min(bmin, bmin, min2, MPFR_RNDD) + + ######################################################################## + # Free memory + + mpfr_clear(tmp) + if crosses_NE_diagonal: + mpfr_clear(min2) def make_ComplexIntervalFieldElement0( fld, re, im ): """ From a39dd3ea78f204ce5ee36934a218aa01e7098a63 Mon Sep 17 00:00:00 2001 From: Matthias Goerner <1239022+unhyperbolic@users.noreply.github.com> Date: Sat, 13 Jul 2024 13:38:41 -0700 Subject: [PATCH 2/3] Updating tests in complex_interval.pyx due to reintroducing tight complex interval division. --- src/sage/rings/complex_interval.pyx | 30 ++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/sage/rings/complex_interval.pyx b/src/sage/rings/complex_interval.pyx index 575829d256c..30b27d44cb5 100644 --- a/src/sage/rings/complex_interval.pyx +++ b/src/sage/rings/complex_interval.pyx @@ -799,16 +799,16 @@ cdef class ComplexIntervalFieldElement(FieldElement): sage: b = CIF(-1, (2, 3)) sage: c = a/b sage: c.endpoints() - (0.200000000000000 - 2.00000000000000*I, - 2.30000000000000 - 0.500000000000000*I, - 0.200000000000000 - 0.500000000000000*I, - 2.30000000000000 - 2.00000000000000*I) + (0.500000000000000 - 1.60000000000000*I, + 1.50000000000000 - 0.600000000000000*I, + 0.500000000000000 - 0.600000000000000*I, + 1.50000000000000 - 1.60000000000000*I) sage: c = b/a sage: c.endpoints() - (0.100000000000000 + 0.250000000000000*I, - 1.15000000000000 + 1.00000000000000*I, - 0.100000000000000 + 1.00000000000000*I, - 1.15000000000000 + 0.250000000000000*I) + (0.246153846153846 + 0.317647058823529*I, + 0.841176470588236 + 0.761538461538462*I, + 0.246153846153846 + 0.761538461538462*I, + 0.841176470588236 + 0.317647058823529*I) """ return self * right.__invert__() @@ -1135,10 +1135,10 @@ cdef class ComplexIntervalFieldElement(FieldElement): sage: a = CIF((1, 2), (3, 4)) sage: c = a.__invert__() sage: c.endpoints() - (0.0500000000000000 - 0.400000000000000*I, - 0.200000000000000 - 0.150000000000000*I, - 0.0500000000000000 - 0.150000000000000*I, - 0.200000000000000 - 0.400000000000000*I) + (0.0588235294117647 - 0.300000000000000*I, + 0.153846153846154 - 0.200000000000000*I, + 0.0588235294117647 - 0.200000000000000*I, + 0.153846153846154 - 0.300000000000000*I) TESTS: @@ -1155,7 +1155,7 @@ cdef class ComplexIntervalFieldElement(FieldElement): Test that the bug reported in :issue:`25414` has been fixed:: sage: 1 / CIF(RIF(-1 ,1), 0) - [-infinity .. +infinity] + [0.0000000000000000 .. +infinity]*I + [.. NaN ..] + [.. NaN ..]*I Test that the bug reported in :issue:`37927` is fixed:: @@ -2013,7 +2013,7 @@ cdef class ComplexIntervalFieldElement(FieldElement): sage: CIF(1,1).tan() 0.27175258531952? + 1.08392332733870?*I sage: CIF(2).tan() - -2.18503986326152? + -2.185039863261519? sage: CIF(0,2).tan() 0.964027580075817?*I """ @@ -2102,7 +2102,7 @@ cdef class ComplexIntervalFieldElement(FieldElement): sage: CIF(2).tanh() 0.964027580075817? sage: CIF(0,2).tanh() - -2.18503986326152?*I + -2.185039863261519?*I """ return self.sinh() / self.cosh() From 32568caef0e486eea6e16c50c00263afd0f63265 Mon Sep 17 00:00:00 2001 From: Matthias Goerner <1239022+unhyperbolic@users.noreply.github.com> Date: Sat, 13 Jul 2024 15:11:45 -0700 Subject: [PATCH 3/3] Updating tests due to reintroducing tight complex interval division. This only touches baselines that were touched in commit 8d59b12a063e9 where the first implementation of tight complex interval division was reverted. This change resets most of the baselines to what they were just before that commit 8d59b12a063e9. --- .../arithmetic_dynamics/projective_ds.py | 2 +- src/sage/matrix/matrix0.pyx | 2 +- src/sage/matrix/matrix2.pyx | 46 ++++++++++--------- .../rings/polynomial/binary_form_reduce.py | 4 +- src/sage/rings/qqbar.py | 14 +++--- src/sage/schemes/curves/zariski_vankampen.py | 2 +- .../schemes/projective/projective_morphism.py | 4 +- .../expression_conversion_algebraic.py | 4 +- 8 files changed, 40 insertions(+), 38 deletions(-) diff --git a/src/sage/dynamics/arithmetic_dynamics/projective_ds.py b/src/sage/dynamics/arithmetic_dynamics/projective_ds.py index a04cf57a28d..223e2b9a3e7 100644 --- a/src/sage/dynamics/arithmetic_dynamics/projective_ds.py +++ b/src/sage/dynamics/arithmetic_dynamics/projective_ds.py @@ -6103,7 +6103,7 @@ def reduced_form(self, **kwds): sage: f.reduced_form(prec=30, smallest_coeffs=False) Traceback (most recent call last): ... - ValueError: accuracy of Newton's root not within tolerance(0.00009... > 1e-06), increase precision + ValueError: accuracy of Newton's root not within tolerance(0.00008... > 1e-06), increase precision sage: f.reduced_form(smallest_coeffs=False) ( Dynamical System of Projective Space of dimension 1 over Rational Field diff --git a/src/sage/matrix/matrix0.pyx b/src/sage/matrix/matrix0.pyx index 83edc74c16d..b792627346b 100644 --- a/src/sage/matrix/matrix0.pyx +++ b/src/sage/matrix/matrix0.pyx @@ -1974,7 +1974,7 @@ cdef class Matrix(sage.structure.element.Matrix): sage: K = (A - e).kernel() sage: P = K.basis_matrix() sage: P.str() - '[ 1.000000000000000? + 0.?e-17*I -2.116651487479748? + 0.0255565807096352?*I -0.2585224251020429? + 0.288602340904754?*I -0.4847545623533090? - 1.871890760086142?*I]' + '[ 1.000000000000000? + 0.?e-17*I -2.116651487479748? + 0.0255565807096352?*I -0.2585224251020429? + 0.2886023409047535?*I -0.4847545623533090? - 1.871890760086142?*I]' Use single-row delimiters where appropriate:: diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 0b51d5c1df7..37042bafa43 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -10589,23 +10589,24 @@ cdef class Matrix(Matrix1): ....: [-1, 1, -6, -6, 5]]) sage: Q, R = A.QR() sage: Q - [ -0.4588314677411235? -0.1260506983326509? 0.3812120831224489? -0.394573711338418? -0.687440062597?] - [ -0.4588314677411235? 0.4726901187474409? -0.05198346588033394? 0.717294125164660? -0.220962877263?] - [ 0.2294157338705618? 0.6617661662464172? 0.6619227988762521? -0.180872093737548? 0.1964114464561?] - [ 0.6882472016116853? 0.1890760474989764? -0.2044682991293135? 0.096630296654307? -0.662888631790?] + [ -0.4588314677411235? -0.1260506983326509? 0.3812120831224489? -0.394573711338418? -0.6874400625964?] + [ -0.4588314677411235? 0.4726901187474409? -0.05198346588033394? 0.7172941251646595? -0.2209628772631?] + [ 0.2294157338705618? 0.6617661662464172? 0.6619227988762521? -0.1808720937375480? 0.1964114464561?] + [ 0.6882472016116853? 0.1890760474989764? -0.2044682991293135? 0.0966302966543065? -0.6628886317894?] [ -0.2294157338705618? 0.5357154679137663? -0.609939332995919? -0.536422031427112? 0.0245514308070?] sage: R [ 4.358898943540674? -0.4588314677411235? 13.07669683062202? 6.194224814505168? 2.982404540317303?] [ 0 1.670171752907625? 0.5987408170800917? -1.292019657909672? 6.207996892883057?] - [ 0 0 5.444401659866974? 5.468660610611130? -0.682716185228386?] - [ 0 0 0 1.027626039419836? -3.61930014968662?] - [ 0 0 0 0 0.02455143080702?] + [ 0 0 5.444401659866974? 5.468660610611130? -0.6827161852283857?] + [ 0 0 0 1.027626039419836? -3.619300149686620?] + [ 0 0 0 0 0.024551430807012?] + sage: Q.conjugate_transpose()*Q - [1.000000000000000? 0.?e-18 0.?e-17 0.?e-15 0.?e-12] - [ 0.?e-18 1.000000000000000? 0.?e-16 0.?e-15 0.?e-12] - [ 0.?e-17 0.?e-16 1.000000000000000? 0.?e-15 0.?e-12] - [ 0.?e-15 0.?e-15 0.?e-15 1.000000000000000? 0.?e-12] - [ 0.?e-12 0.?e-12 0.?e-12 0.?e-12 1.000000000000?] + [1.000000000000000? 0.?e-18 0.?e-17 0.?e-16 0.?e-13] + [ 0.?e-18 1.000000000000000? 0.?e-17 0.?e-16 0.?e-13] + [ 0.?e-17 0.?e-17 1.000000000000000? 0.?e-16 0.?e-13] + [ 0.?e-16 0.?e-16 0.?e-16 1.000000000000000? 0.?e-13] + [ 0.?e-13 0.?e-13 0.?e-13 0.?e-13 1.0000000000000?] sage: Q * R == A True @@ -10621,24 +10622,25 @@ cdef class Matrix(Matrix1): sage: Q, R = A.QR() sage: Q [ -0.7302967433402215? 0.2070566455055649? + 0.5383472783144687?*I 0.2463049809998642? - 0.0764456358723292?*I 0.2381617683194332? - 0.1036596032779695?*I] - [ 0.0912870929175277? -0.2070566455055649? - 0.3778783780476559?*I 0.3786559533863032? - 0.1952221495524667?*I 0.701244450214469? - 0.364371165098660?*I] - [ 0.6390096504226938? + 0.0912870929175277?*I 0.1708217325420910? + 0.6677576817554466?*I -0.03411475806452072? + 0.04090198741767143?*I 0.3140171085506763? - 0.0825191718705412?*I] + [ 0.0912870929175277? -0.2070566455055649? - 0.3778783780476559?*I 0.3786559533863033? - 0.1952221495524667?*I 0.701244450214469? - 0.3643711650986595?*I] + [ 0.6390096504226938? + 0.0912870929175277?*I 0.1708217325420910? + 0.6677576817554466?*I -0.03411475806452072? + 0.04090198741767143?*I 0.3140171085506764? - 0.0825191718705412?*I] [ 0.1825741858350554? + 0.0912870929175277?*I -0.03623491296347385? + 0.0724698259269477?*I 0.8632284069415110? + 0.06322839976356195?*I -0.4499694867611521? - 0.0116119181208918?*I] sage: R [ 10.95445115010333? 0.?e-18 - 1.917028951268082?*I 5.385938482134133? - 2.190890230020665?*I -0.2738612787525831? - 2.190890230020665?*I] - [ 0 4.829596256417300? + 0.?e-17*I -0.869637911123373? - 5.864879483945125?*I 0.993871898426712? - 0.3054085521207082?*I] + [ 0 4.829596256417300? + 0.?e-18*I -0.869637911123373? - 5.864879483945125?*I 0.993871898426712? - 0.3054085521207082?*I] [ 0 0 12.00160760935814? + 0.?e-16*I -0.2709533402297273? + 0.4420629644486323?*I] [ 0 0 0 1.942963944258992? + 0.?e-16*I] sage: Q.conjugate_transpose()*Q [1.000000000000000? + 0.?e-19*I 0.?e-18 + 0.?e-17*I 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I] [ 0.?e-18 + 0.?e-17*I 1.000000000000000? + 0.?e-17*I 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I] - [ 0.?e-17 + 0.?e-17*I 0.?e-17 + 0.?e-17*I 1.000000000000000? + 0.?e-16*I 0.?e-16 + 0.?e-16*I] - [ 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I 1.000000000000000? + 0.?e-15*I] + [ 0.?e-17 + 0.?e-17*I 0.?e-17 + 0.?e-17*I 1.000000000000000? + 0.?e-17*I 0.?e-16 + 0.?e-16*I] + [ 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I 1.000000000000000? + 0.?e-16*I] + sage: Q*R - A [ 0.?e-17 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I] - [ 0.?e-18 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-15 + 0.?e-15*I] + [ 0.?e-18 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I] [0.?e-17 + 0.?e-18*I 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I] - [0.?e-18 + 0.?e-18*I 0.?e-18 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-15 + 0.?e-16*I] + [0.?e-18 + 0.?e-18*I 0.?e-18 + 0.?e-18*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I] A rank-deficient rectangular matrix, with both values of the ``full`` keyword. :: @@ -12328,9 +12330,9 @@ cdef class Matrix(Matrix1): sage: # needs sage.combinat sage.libs.pari sage: _, T = A.is_similar(B, transformation=True) sage: T - [ 1.0000000000000? + 0.?e-13*I 0.?e-13 + 0.?e-13*I 0.?e-13 + 0.?e-13*I] - [-0.6666666666667? + 0.?e-13*I 0.16666666666667? + 0.?e-14*I -0.8333333333334? + 0.?e-13*I] - [ 0.6666666666667? + 0.?e-13*I 0.?e-13 + 0.?e-13*I -0.333333333334? + 0.?e-13*I] + [ 1.00000000000000? + 0.?e-14*I 0.?e-14 + 0.?e-14*I 0.?e-14 + 0.?e-14*I] + [-0.66666666666667? + 0.?e-15*I 0.166666666666667? + 0.?e-15*I -0.83333333333334? + 0.?e-14*I] + [ 0.66666666666667? + 0.?e-14*I 0.?e-14 + 0.?e-14*I -0.33333333333333? + 0.?e-14*I] sage: T.change_ring(QQ) [ 1 0 0] [-2/3 1/6 -5/6] diff --git a/src/sage/rings/polynomial/binary_form_reduce.py b/src/sage/rings/polynomial/binary_form_reduce.py index 6afb11155cf..8b6134ec59b 100644 --- a/src/sage/rings/polynomial/binary_form_reduce.py +++ b/src/sage/rings/polynomial/binary_form_reduce.py @@ -79,8 +79,8 @@ def covariant_z0(F, z0_cov=False, prec=53, emb=None, error_limit=0.000001): sage: F = -x^8 + 6*x^7*y - 7*x^6*y^2 - 12*x^5*y^3 + 27*x^4*y^4\ ....: - 4*x^3*y^5 - 19*x^2*y^6 + 10*x*y^7 - 5*y^8 sage: covariant_z0(F, prec=80) - (0.64189877107807122203369 + 1.1852516565091601348355*I, - 3134.5148284344627168275) + (0.64189877107807122203366 + 1.1852516565091601348355*I, + 3134.5148284344627168276) :: diff --git a/src/sage/rings/qqbar.py b/src/sage/rings/qqbar.py index 73e8878150c..0122ab20ab3 100644 --- a/src/sage/rings/qqbar.py +++ b/src/sage/rings/qqbar.py @@ -543,14 +543,14 @@ ....: (83132288614462248899077910195645998*a + 254420831388756945210526696075302411552001)*y^2) sage: lc = lc.change_ring(QQbar) sage: lc.roots(CIF) - [(-1.0005054922387982573627768714?, 2), - (-1.0000000000000000000000000000?, 2), - (-0.9999999996626050731848036720993?, 1), + [(-1.000505492239?, 2), + (-1.000000000000?, 2), + (-0.999999999662605?, 1), (0, 2), - (1.0000000000000000000000000000?, 2), - (1.0005054922387982573627768714?, 2), - (0.999999586737109168744488? + 0.?e-43*I, 1), - (0.999999999662605073184804? + 0.?e-43*I, 1)] + (1.000000000000?, 2), + (1.000505492239?, 2), + (0.999999587? + 0.?e-11*I, 1), + (0.999999999? + 0.?e-11*I, 1)] Check that issue:`37927` is fixed:: diff --git a/src/sage/schemes/curves/zariski_vankampen.py b/src/sage/schemes/curves/zariski_vankampen.py index 7dec97439db..8260b3db92f 100644 --- a/src/sage/schemes/curves/zariski_vankampen.py +++ b/src/sage/schemes/curves/zariski_vankampen.py @@ -608,7 +608,7 @@ def newton(f, x0, i0): sage: n 0.0? sage: n.real().endpoints() - (-0.0460743801652894, 0.0291454081632654) + (-0.0147727272727274, 0.00982142857142862) sage: n.imag().endpoints() (0.000000000000000, -0.000000000000000) """ diff --git a/src/sage/schemes/projective/projective_morphism.py b/src/sage/schemes/projective/projective_morphism.py index 339406cd9d6..2baab7e5cda 100644 --- a/src/sage/schemes/projective/projective_morphism.py +++ b/src/sage/schemes/projective/projective_morphism.py @@ -1764,10 +1764,10 @@ def _number_field_from_algebraics(self): Scheme morphism: From: Projective Space of dimension 1 over Number Field in a with defining polynomial y^4 + 3*y^2 + 1 - with a = 0.?e-151 + 0.618033988749895?*I + with a = 0.?e-113 + 0.618033988749895?*I To: Projective Space of dimension 2 over Number Field in a with defining polynomial y^4 + 3*y^2 + 1 - with a = 0.?e-151 + 0.618033988749895?*I + with a = 0.?e-113 + 0.618033988749895?*I Defn: Defined on coordinates by sending (x : y) to (x^2 + (a^3 + 2*a)*x*y + 3*y^2 : y^2 : (2*a^2 + 3)*x*y) diff --git a/src/sage/symbolic/expression_conversion_algebraic.py b/src/sage/symbolic/expression_conversion_algebraic.py index 11314aaeb70..a8b932c75f2 100644 --- a/src/sage/symbolic/expression_conversion_algebraic.py +++ b/src/sage/symbolic/expression_conversion_algebraic.py @@ -146,9 +146,9 @@ def composition(self, ex, operator): sage: a.composition(complex_root_of(x^5 - 1, 3), complex_root_of) 0.3090169943749474? - 0.9510565162951536?*I sage: a.composition(complex_root_of(x^2 + 1, 0), complex_root_of) - 1.?e-884 - 0.9999999999999999?*I + 1.?e-683 - 0.9999999999999999?*I sage: a.composition(complex_root_of(x^2 + 1, 1), complex_root_of) - 1.?e-884 + 0.9999999999999999?*I + 1.?e-683 + 0.9999999999999999?*I TESTS::