#28464 : fix the construction of a simplex of vertices in Polyhedron_base.is_inscribed

Author: Thierry Monteil

src/sage/geometry/polyhedron/base.py

@@ -2553,6 +2553,55 @@ def is_inscribed(self, certificate=False):
             Traceback (most recent call last):
             ...
             NotImplementedError: this function is not implemented for unbounded polyhedron
+
+        TESTS:
+
+        We check that :trac:`28464` is fixed::
+
+            sage: P = Polyhedron(vertices=[(-130658298093891402635075/416049251842505144482473,
+            ....: 177469511761879509172000/1248147755527515433447419,
+            ....: 485550543257132133136169/2496295511055030866894838,
+            ....: 2010744967797898733758669/2496295511055030866894838),
+            ....: (-146945725603929909850/706333405676769433081,
+            ....: -84939725782618445000/706333405676769433081,
+            ....: 560600045283000988081/1412666811353538866162,
+            ....: 969778382942371268081/1412666811353538866162),
+            ....: (-46275018824497300/140422338198040641,
+            ....: -5747688262110000/46807446066013547, 1939357556329/7033601552658,
+            ....: 1939357556329/7033601552658), (-17300/59929, -10000/59929, 39929/119858,
+            ....: 39929/119858), (-4700/32209, -10000/32209, 12209/64418, 12209/64418),
+            ....: (QQ(0), QQ(0), QQ(0), QQ(1)), (QQ(0), QQ(0), 1/2, 1/2), (300/10027,
+            ....: -10000/30081, 10081/60162, 10081/60162), (112393975400/1900567733649,
+            ....: 117311600000/633522577883, 43678681/95197362, 43678681/95197362),
+            ....: (6109749955400/133380598418321, 37106807920000/133380598418321,
+            ....: 2677964249/6680888498, 2677964249/6680888498),
+            ....: (29197890764005600/402876806828660641,
+            ....: -2150510776960000/402876806828660641,
+            ....: 398575785274740641/805753613657321282,
+            ....: 398575785274740641/805753613657321282),
+            ....: (5576946899441759759983005325/110078073300232813237456943251,
+            ....: -29071211718677797926570478000/110078073300232813237456943251,
+            ....: 59439312069347378584317232001/220156146600465626474913886502,
+            ....: 181346577228466312205473034501/220156146600465626474913886502),
+            ....: (150040732779124914266530235300/6774574358246204311268446913881,
+            ....: -2813827375989039189507000218000/6774574358246204311268446913881,
+            ....: 1260217414021285074925933133881/13549148716492408622536893827762,
+            ....: 3232518047094242684574253773881/13549148716492408622536893827762),
+            ....: (3816349407976279597850158016285000/88842127448735433741180809504357161,
+            ....: 27965821247423216557301387453968000/88842127448735433741180809504357161,
+            ....: 68546256000224819256028677086357161/177684254897470867482361619008714322,
+            ....: 86062257922545755787315412690197161/177684254897470867482361619008714322)])
+            sage: P.is_inscribed()
+            True
+
+            sage: P = Polyhedron(vertices=[[0, -1, 0, 0],
+            ....:                          [0, 0, -1, 0],
+            ....:                          [0, 0, 0, -1],
+            ....:                          [0, 0, +1, 0],
+            ....:                          [0, 0, 0, +1],
+            ....:                          [+1, 0, 0, 0]])
+            sage: P.is_inscribed()
+            True
         """
 
         if not self.is_compact():
@@ -2563,13 +2612,20 @@ def is_inscribed(self, certificate=False):
 
         dimension = self.dimension()
         vertices = self.vertices()
-        vertex = vertices[0]
-        vertex_neighbors = vertex.neighbors()
 
-        # The following simplex is full-dimensional because `self` is assumed
-        # to be: every vertex has at least `dimension` neighbors and they form
-        # a full simplex with `vertex`.
-        simplex_vertices = [vertex] + [next(vertex_neighbors) for i in range(dimension)]
+        # We greedily construct a full-dimensional simplex made of vertices
+        # around some vertex of self.
+        vertex = vertices[0]
+        v0 = vertex.vector()
+        M = matrix(self.base_ring(), dimension, 0)
+        simplex_vertices = [vertex]
+        for v in vertex.neighbors():
+            MA = M.augment(v.vector()-v0)
+            if MA.rank() > M.rank():
+                simplex_vertices.append(v)
+                M = MA
+            if M.rank() == dimension:
+                break
 
         raw_data = []
         for vertex in simplex_vertices: