Trac #32152: Set up pyramid with both Vrep and Hrep

We set up the pyramid over a polyhedron with the double description.

Before:

{{{
sage: P = polytopes.permutahedron(6)
sage: %time Q = P.pyramid()
CPU times: user 99.4 ms, sys: 0 ns, total: 99.4 ms
Wall time: 98 ms
sage: P = polytopes.hypercube(8)
sage: %time Q = P.pyramid()
CPU times: user 34.2 ms, sys: 24 µs, total: 34.2 ms
Wall time: 33.3 ms
sage: P = polytopes.cross_polytope(8)
sage: %time Q = P.pyramid()
CPU times: user 40.5 ms, sys: 3.96 ms, total: 44.5 ms
Wall time: 43.3 ms
sage: P = polytopes.hypercube(6, backend='field')
sage: %time Q = P.pyramid()
CPU times: user 974 ms, sys: 0 ns, total: 974 ms
Wall time: 973 ms
}}}

After:
{{{
sage: P = polytopes.permutahedron(6)
sage: %time Q = P.pyramid()
CPU times: user 119 ms, sys: 197 µs, total: 120 ms
Wall time: 118 ms
sage: P = polytopes.hypercube(8)
sage: %time Q = P.pyramid()
CPU times: user 30.5 ms, sys: 3.92 ms, total: 34.4 ms
Wall time: 33.4 ms
sage: P = polytopes.cross_polytope(8)
sage: %time Q = P.pyramid()
CPU times: user 40.4 ms, sys: 0 ns, total: 40.4 ms
Wall time: 39.7 ms
sage: P = polytopes.hypercube(6, backend='field')
sage: %time Q = P.pyramid()
CPU times: user 8.04 ms, sys: 0 ns, total: 8.04 ms
Wall time: 7.99 ms
sage: P = polytopes.hypercube(10, backend='field')
sage: %time Q = P.pyramid()
CPU times: user 46.4 ms, sys: 98 µs, total: 46.5 ms
Wall time: 45.6 ms
}}}

URL: https://trac.sagemath.org/32152
Reported by: gh-kliem
Ticket author(s): Jonathan Kliem
Reviewer(s): Travis Scrimshaw

src/sage/geometry/polyhedron/base.py

@@ -7819,13 +7819,25 @@ def pyramid(self):
             'cdd'
         """
         assert self.is_compact(), "Not a polytope."
+        c = self.center()
 
-        new_verts = \
-            [[0] + x for x in self.Vrep_generator()] + \
-            [[1] + list(self.center())]
+        from itertools import chain
+        new_verts = chain(([0] + x for x in self.Vrep_generator()),
+                          [[1] + list(c)])
+        new_ieqs = chain(([i.b()] + [-c*i.A() - i.b()] + list(i.A()) for i in self.inequalities()),
+                         [[0, 1] + [0]*self.ambient_dim()])
+        new_eqns = ([e.b()] + [0] + list(e.A()) for e in self.equations())
 
+        pref_rep = 'Hrep' if self.n_vertices() > self.n_inequalities() else 'Vrep'
         parent = self.parent().base_extend(self.center().parent(), ambient_dim=self.ambient_dim()+1)
-        return parent.element_class(parent, [new_verts, [], []], None)
+
+        if self.n_vertices() == 1:
+            # Fix the polyhedron with one vertex.
+            return parent.element_class(parent, [new_verts, [], []], None)
+
+        return parent.element_class(parent, [new_verts, [], []],
+                                    [new_ieqs, new_eqns],
+                                    Vrep_minimal=True, Hrep_minimal=True, pref_rep=pref_rep)
 
     def _test_pyramid(self, tester=None, **options):
         """
@@ -7875,6 +7887,18 @@ def check_pyramid_certificate(P, cert):
 
             tester.assertTrue(pyr_polar.is_combinatorially_isomorphic(pyr_polar))
 
+            # Basic properties of the pyramid.
+
+            # Check that the prism preserves the backend.
+            tester.assertEqual(pyr.backend(), self.backend())
+
+            tester.assertEqual(1 + self.n_vertices(), pyr.n_vertices())
+            tester.assertEqual(self.n_equations(), pyr.n_equations())
+            tester.assertEqual(1 + self.n_inequalities(), pyr.n_inequalities())
+
+            if self.n_vertices() < 15 and self.n_facets() < 15:
+                pyr._test_basic_properties()
+
     def bipyramid(self):
         """
         Return a polyhedron that is a bipyramid over the original.