Trac #27973: Implement wedge over a face of Polyhedron

From ​​https://www.csun.edu/~ctoth/Handbook/chap15.pdf:

The wedge over a facet `F` of a polytope `P` is defined as:

`(P \times \mathbb{R}) \cap \{a^\top x +|x_{d+1}| \leq b\}`

where `F` is a facet defined by `a^\top x leq b`.

It  has  dimension `d+1`, `m+1` facets, and `2n-n_F` vertices, if `F`
has `n_F` vertices. More generally,  the wedge construction can be
performed (defined by the same formula) for a face `F`.

URL: https://trac.sagemath.org/27973
Reported by: jipilab
Ticket author(s): Laith Rastanawi, Jonathan Kliem
Reviewer(s): Jean-Philippe Labbé

src/doc/en/reference/references/index.rst

@@ -2219,6 +2219,10 @@ REFERENCES:
             and Genocchi numbers <https://www.lri.fr/~hivert/PAPER/kshapes.pdf>`_,
             in FPSAC 2011, Reykjav´k, Iceland DMTCS proc. AO, 2011, 493-504.
 
+.. [HoDaCG17] Toth, Csaba D., Joseph O'Rourke, and Jacob E. Goodman.
+              Handbook of Discrete and Computational Geometry (3rd Edition).
+              Chapman and Hall/CRC, 2017.
+
 .. [Hoc] Winfried Hochstaettler, "About the Tic-Tac-Toe Matroid",
          preprint.
 
@@ -2633,7 +2637,7 @@ REFERENCES:
              Compositio Mathematica, **149** (2013), no. 10.
              :arxiv:`1111.3660`.
 
-.. [Kly1990] Klyachko, Aleksandr Anatolevich. 
+.. [Kly1990] Klyachko, Aleksandr Anatolevich.
              Equivariant Bundles on Toral Varieties,
              Math USSR Izv. 35 (1990), 337-375.
              http://iopscience.iop.org/0025-5726/35/2/A04/pdf/0025-5726_35_2_A04.pdf

src/doc/en/thematic_tutorials/geometry/polyhedra_quickref.rst

@@ -114,6 +114,7 @@ List of Polyhedron methods
     :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.truncation` | truncates all vertices simultaneously
     :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.lawrence_extension` | returns the Lawrence extension of self on a given point
     :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.lawrence_polytope` | returns the Lawrence polytope of self
+    :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.wedge` | returns the wedge over a face of self
 
 **Combinatorics**
 

src/sage/geometry/polyhedron/base.py

@@ -3060,14 +3060,14 @@ def normal_fan(self, direction='inner'):
         r"""
         Return the normal fan of a compact full-dimensional rational polyhedron.
 
-        This returns the inner normal fan of ``self``. For the outer normal fan, 
+        This returns the inner normal fan of ``self``. For the outer normal fan,
         use ``direction='outer'``.
 
         INPUT:
 
-        - ``direction`` -- either ``'inner'`` (default) or ``'outer'``; if 
-          set to ``'inner'``, use the inner normal vectors to span the cones of 
-          the fan, if set to ``'outer'``, use the outer normal vectors. 
+        - ``direction`` -- either ``'inner'`` (default) or ``'outer'``; if
+          set to ``'inner'``, use the inner normal vectors to span the cones of
+          the fan, if set to ``'outer'``, use the outer normal vectors.
 
         OUTPUT:
 
@@ -4458,6 +4458,139 @@ def stack(self, face, position=None):
         parent = self.parent().base_extend(new_vertex)
         return parent.element_class(parent, [self.vertices() + (new_vertex,), self.rays(), self.lines()], None)
 
+    def wedge(self, face, width=1):
+        r"""
+        Return the wedge over a ``face`` of the polytope ``self``.
+
+        The wedge over a face `F` of a polytope `P` with width `w \not= 0`
+        is defined as:
+
+        .. MATH::
+
+            (P \times \mathbb{R}) \cap \{a^\top x + |w x_{d+1}| \leq b\}
+
+        where `\{x | a^\top x = b\}` is a supporting hyperplane defining `F`.
+
+        INPUT:
+
+        - ``face`` -- a PolyhedronFace of ``self``, the face which we take
+          the wedge over
+        - ``width`` -- a nonzero number (default: ``1``);
+          specifies how wide the wedge will be
+
+        OUTPUT:
+
+        A (bounded) polyhedron
+
+        EXAMPLES::
+
+            sage: P_4 = polytopes.regular_polygon(4)
+            sage: W1 = P_4.wedge(P_4.faces(1)[0]); W1
+            A 3-dimensional polyhedron in AA^3 defined as the convex hull of 6 vertices
+            sage: triangular_prism = polytopes.regular_polygon(3).prism()
+            sage: W1.is_combinatorially_isomorphic(triangular_prism)
+            True
+
+            sage: Q = polytopes.hypersimplex(4,2)
+            sage: W2 = Q.wedge(Q.faces(2)[0]); W2
+            A 4-dimensional polyhedron in QQ^5 defined as the convex hull of 9 vertices
+            sage: W2.vertices()
+            (A vertex at (0, 1, 0, 1, 0),
+             A vertex at (0, 0, 1, 1, 0),
+             A vertex at (1, 0, 0, 1, -1),
+             A vertex at (1, 0, 0, 1, 1),
+             A vertex at (1, 0, 1, 0, 1),
+             A vertex at (1, 1, 0, 0, -1),
+             A vertex at (0, 1, 1, 0, 0),
+             A vertex at (1, 0, 1, 0, -1),
+             A vertex at (1, 1, 0, 0, 1))
+
+            sage: W3 = Q.wedge(Q.faces(1)[0]); W3
+            A 4-dimensional polyhedron in QQ^5 defined as the convex hull of 10 vertices
+            sage: W3.vertices()
+            (A vertex at (0, 1, 0, 1, 0),
+             A vertex at (0, 0, 1, 1, 0),
+             A vertex at (1, 0, 0, 1, -1),
+             A vertex at (1, 0, 0, 1, 1),
+             A vertex at (1, 0, 1, 0, 2),
+             A vertex at (0, 1, 1, 0, 1),
+             A vertex at (1, 0, 1, 0, -2),
+             A vertex at (1, 1, 0, 0, 2),
+             A vertex at (0, 1, 1, 0, -1),
+             A vertex at (1, 1, 0, 0, -2))
+
+            sage: C_3_7 = polytopes.cyclic_polytope(3,7)
+            sage: P_6 = polytopes.regular_polygon(6)
+            sage: W4 = P_6.wedge(P_6.faces(1)[0])
+            sage: W4.is_combinatorially_isomorphic(C_3_7.polar())
+            True
+
+        REFERENCES:
+
+        For more information, see Chapter 15 of [HoDaCG17]_.
+
+        TESTS::
+
+        The backend should be preserved as long as the value of width permits.
+        The base_ring will change to the field of fractions of the current
+        base_ring, unless width forces a different ring.
+
+            sage: P = polytopes.cyclic_polytope(3,7, base_ring=ZZ, backend='field')
+            sage: W1 = P.wedge(P.faces(2)[0]); W1.base_ring(); W1.backend()
+            Rational Field
+            'field'
+            sage: W2 = P.wedge(P.faces(2)[0], width=5/2); W2.base_ring(); W2.backend()
+            Rational Field
+            'field'
+            sage: W2 = P.wedge(P.faces(2)[0], width=4/2); W2.base_ring(); W2.backend()
+            Rational Field
+            'field'
+            sage: W2.vertices()
+            (A vertex at (3, 9, 27, -1/2),
+             A vertex at (4, 16, 64, -2),
+             A vertex at (6, 36, 216, -10),
+             A vertex at (5, 25, 125, -5),
+             A vertex at (2, 4, 8, 0),
+             A vertex at (1, 1, 1, 0),
+             A vertex at (0, 0, 0, 0),
+             A vertex at (3, 9, 27, 1/2),
+             A vertex at (4, 16, 64, 2),
+             A vertex at (6, 36, 216, 10),
+             A vertex at (5, 25, 125, 5))
+            sage: W2 = P.wedge(P.faces(2)[0], width=1.0); W2.base_ring(); W2.backend()
+            Real Double Field
+            'cdd'
+        """
+        width = width*ZZ.one()
+
+        if not self.is_compact():
+            raise ValueError("polyhedron 'self' must be a polytope")
+
+        if width == 0:
+            raise ValueError("the width should be nonzero")
+
+        from sage.geometry.polyhedron.face import PolyhedronFace
+        if not isinstance(face, PolyhedronFace):
+            raise TypeError("{} should be a PolyhedronFace of {}".format(face, self))
+
+        F_Hrep = vector([0]*(self.ambient_dim()+1))
+        for facet in face.ambient_Hrepresentation():
+            if facet.is_inequality():
+                F_Hrep = F_Hrep + facet.vector()
+        F_Hrep = list(F_Hrep)
+
+        # Preserve the backend, if value of ``width`` permits.
+        backend = None
+        from .parent import does_backend_handle_base_ring
+        if does_backend_handle_base_ring(width.base_ring().fraction_field(), self.backend()):
+            backend = self.backend()
+
+        L = Polyhedron(lines=[[1]])
+        Q = self.product(L)
+        ieqs = [F_Hrep + [width], F_Hrep + [-width]]
+        H = Polyhedron(ieqs=ieqs, backend=backend)
+        return Q.intersection(H)
+
     def lawrence_extension(self, v):
         """
         Return the Lawrence extension of ``self`` on the point ``v``.