Minkowski sum of regular polygons fails for backend field

[kliem]

Duplicate of #28599.

sage: polytopes.regular_polygon(3) + polytopes.regular_polygon(8)

give assertion error.

Full traceback below.

It appears to be a rather new bug, maybe related to python3. It stills seems to work on 8.9.

Component: geometry

Keywords: polytopes, field, regular polygons

Issue created by migration from https://trac.sagemath.org/ticket/28804

[kliem]

comment:1

sage: polytopes.regular_polygon(3) + polytopes.regular_polygon(8)
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-1-bfd7b97fb44d> in <module>()
----> 1 polytopes.regular_polygon(Integer(3)) + polytopes.regular_polygon(Integer(8))

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/structure/element.pyx in sage.structure.element.Element.__add__ (build/cythonized/sage/structure/element.c:10798)()
   1229         cdef int cl = classify_elements(left, right)
   1230         if HAVE_SAME_PARENT(cl):
-> 1231             return (<Element>left)._add_(right)
   1232         # Left and right are Sage elements => use coercion model
   1233         if BOTH_ARE_ELEMENT(cl):

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/structure/element.pyx in sage.structure.element.Element._add_ (build/cythonized/sage/structure/element.c:11148)()
   1281             raise bin_op_exception('+', self, other)
   1282         else:
-> 1283             return python_op(other)
   1284 
   1285     cdef _add_long(self, long n):

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/structure/element.pyx in sage.structure.element.coerce_binop.new_method (build/cythonized/sage/structure/element.c:26671)()
   4274     def new_method(self, other, *args, **kwargs):
   4275         if have_same_parent(self, other):
-> 4276             return method(self, other, *args, **kwargs)
   4277         else:
   4278             a, b = coercion_model.canonical_coercion(self, other)

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/geometry/polyhedron/base.py in minkowski_sum(self, other)
   3653             new_rays = self.rays() + other.rays()
   3654             new_lines = self.lines() + other.lines()
-> 3655             return self.parent().element_class(self.parent(), [new_vertices, new_rays, new_lines], None)
   3656         else:
   3657             return self.parent().element_class(self.parent(), None, None)

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/geometry/polyhedron/backend_field.py in __init__(self, parent, Vrep, Hrep, Vrep_minimal, Hrep_minimal, **kwds)
    174             self._init_Hrepresentation(*Hrep)
    175         else:
--> 176             super(Polyhedron_field, self).__init__(parent, Vrep, Hrep, **kwds)
    177 
    178     def _init_from_Vrepresentation(self, vertices, rays, lines,

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/geometry/polyhedron/base.py in __init__(self, parent, Vrep, Hrep, **kwds)
    119         if Vrep is not None:
    120             vertices, rays, lines = Vrep
--> 121             self._init_from_Vrepresentation(vertices, rays, lines, **kwds)
    122         elif Hrep is not None:
    123             ieqs, eqns = Hrep

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/geometry/polyhedron/backend_field.py in _init_from_Vrepresentation(self, vertices, rays, lines, minimize, verbose)
    207         H = Vrep2Hrep(self.base_ring(), self.ambient_dim(), vertices, rays, lines)
    208         V = Hrep2Vrep(self.base_ring(), self.ambient_dim(),
--> 209                       H.inequalities, H.equations)
    210         self._init_Vrepresentation_backend(V)
    211         self._init_Hrepresentation_backend(H)

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/geometry/polyhedron/double_description_inhomogeneous.py in __init__(self, base_ring, dim, inequalities, equations)
    204         equations = [list(x) for x in equations]
    205         A = self._init_Vrep(inequalities, equations)
--> 206         DD = Algorithm(A).run()
    207         self._extract_Vrep(DD)
    208         if VERIFY_RESULT:

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/geometry/polyhedron/double_description.py in run(self)
    762         DD, remaining = self.initial_pair()
    763         for a in remaining:
--> 764             DD.add_inequality(a)
    765             if VERIFY_RESULT:
    766                 DD.verify()

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/geometry/polyhedron/double_description.py in add_inequality(self, a)
    713         R_new = []
    714         for rp, rn in itertools.product(R_pos, R_neg):
--> 715             if not self.are_adjacent(rp, rn):
    716                 continue
    717             r = a.inner_product(rp) * rn - a.inner_product(rn) * rp

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/geometry/polyhedron/double_description.py in are_adjacent(self, r1, r2)
    445             False
    446         """
--> 447         Z = self.zero_set(r1).intersection(self.zero_set(r2))
    448         if not Z:
    449             return self.problem.dim() == 2

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/geometry/polyhedron/double_description.py in zero_set(self, ray)
    374         n, t = self.zero_set_cache[ray]
    375         if n != len(self.A):
--> 376             t.update(self.A[i] for i in range(n,len(self.A)) if self.A[i].inner_product(ray) == self.zero)
    377             self.zero_set_cache[ray] = (len(self.A), t)
    378         return t

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/geometry/polyhedron/double_description.py in <genexpr>(.0)
    374         n, t = self.zero_set_cache[ray]
    375         if n != len(self.A):
--> 376             t.update(self.A[i] for i in range(n,len(self.A)) if self.A[i].inner_product(ray) == self.zero)
    377             self.zero_set_cache[ray] = (len(self.A), t)
    378         return t

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/structure/element.pyx in sage.structure.element.Element.__richcmp__ (build/cythonized/sage/structure/element.c:9939)()
   1087             # an instance of Element. The explicit casts below make
   1088             # Cython generate optimized code for this call.
-> 1089             return (<Element>self)._richcmp_(other, op)
   1090         else:
   1091             return coercion_model.richcmp(self, other, op)

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/structure/element.pyx in sage.structure.element.Element._richcmp_ (build/cythonized/sage/structure/element.c:10046)()
   1091             return coercion_model.richcmp(self, other, op)
   1092 
-> 1093     cpdef _richcmp_(left, right, int op):
   1094         r"""
   1095         Default implementation of rich comparisons for elements with

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/rings/qqbar.py in _richcmp_(self, other, op)
   5156                 return bool(other) == (op == op_NE)
   5157             elif type(od) is ANRational and not od._value:
-> 5158                 return bool(self) == (op == op_NE)
   5159             elif (type(sd) is ANExtensionElement and
   5160                   type(od) is ANExtensionElement and

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/rings/qqbar.py in __bool__(self)
   3844                 return True
   3845 
-> 3846             c = cmp_elements_with_same_minpoly(left, right, left.minpoly())
   3847             if c is not None:
   3848                 if c == 0:

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/rings/qqbar.py in cmp_elements_with_same_minpoly(a, b, p)
   2665     else:
   2666         ring = QQbar
-> 2667     roots = p.roots(ring, False)
   2668 
   2669     real = ar.union(br)

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/rings/polynomial/polynomial_element.pyx in sage.rings.polynomial.polynomial_element.Polynomial.roots (build/cythonized/sage/rings/polynomial/polynomial_element.c:61446)()
   7759 
   7760                 if is_AlgebraicRealField(L):
-> 7761                     rts = real_roots(self, retval='algebraic_real')
   7762                 else:
   7763                     diam = ~(ZZ(1) << L.prec())

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.real_roots (build/cythonized/sage/rings/polynomial/real_roots.c:44224)()
   4052 
   4053             oc = ocean(ctx, bernstein_polynomial_factory_ratlist(b), linear_map(left, right))
-> 4054             oc.find_roots()
   4055             oceans.append((oc, factor, exp))
   4056 

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.ocean.find_roots (build/cythonized/sage/rings/polynomial/real_roots.c:33086)()
   3068         """
   3069         while not self.all_done():
-> 3070             self.refine_all()
   3071             self.increase_precision()
   3072 

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.ocean.refine_all (build/cythonized/sage/rings/polynomial/real_roots.c:33381)()
   3113         while isle is not self.endpoint:
   3114             if not isle.done(self.ctx):
-> 3115                 isle.refine(self.ctx)
   3116             isle = isle.rgap.risland
   3117 

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.island.refine (build/cythonized/sage/rings/polynomial/real_roots.c:35696)()
   3366         self.shrink_bp(ctx)
   3367         try:
-> 3368             self.refine_recurse(ctx, self.bp, self.ancestors, [], True)
   3369         except PrecisionError:
   3370             pass

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.island.refine_recurse (build/cythonized/sage/rings/polynomial/real_roots.c:36635)()
   3470                     rv = bp.try_rand_split(ctx, [])
   3471                 if rv is None:
-> 3472                     (ancestors, bp) = self.more_bits(ctx, ancestors, bp, rightmost)
   3473                     if rv is None:
   3474                         rv = bp.try_split(ctx, [])

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.island.more_bits (build/cythonized/sage/rings/polynomial/real_roots.c:39191)()
   3616                     assert(rel_bounds[1] == 1)
   3617 
-> 3618                 ancestor_val = split_for_targets(ctx, ancestor_val, [(self.lgap, maybe_rgap, target_lsb_h)])[0]
   3619 #                 if rel_lbounds[1] > 0:
   3620 #                     left_split = -exact_rational(simple_wordsize_float(-rel_lbounds[1], -rel_lbounds[0]))

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.split_for_targets (build/cythonized/sage/rings/polynomial/real_roots.c:31921)()
   2909             max_lsb = max([t[2] for t in tl1])
   2910             p1 = p1.down_degree_iter(ctx, max_lsb)
-> 2911         r1 = split_for_targets(ctx, p1, tl1)
   2912     else:
   2913         r1 = []

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.split_for_targets (build/cythonized/sage/rings/polynomial/real_roots.c:31126)()
   2879     split = wordsize_rational(split_targets[best_index][0], split_targets[best_index][1], ctx.wordsize)
   2880 
-> 2881     (p1_, p2_, ok) = bp.de_casteljau(ctx, split, msign=split_targets[best_index][2])
   2882     assert(ok)
   2883 

/srv/public/kliem/sage/local/lib/python3.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.interval_bernstein_polynomial_integer.de_casteljau (build/cythonized/sage/rings/polynomial/real_roots.c:9108)()
    727             msign = sign
    728         elif sign != 0:
--> 729             assert(msign == sign)
    730 
    731         cdef Rational absolute_mid = self.lower + mid * (self.upper - self.lower)

AssertionError: 

[kliem]

comment:3

Sorry. I was even CC of that…