Trac #22578: Polyhedron.bounding_box: New keyword argument integral_h…

…ull, use it in .integral_points

An amazing insight is that if one wants to know to bounding box of the
integer points, rather than of all points, then one can round down upper
bounds and round up lower bounds.

This can make a huge difference for integer points enumeration.

URL: https://trac.sagemath.org/22578
Reported by: mkoeppe
Ticket author(s): Matthias Koeppe
Reviewer(s): Travis Scrimshaw

src/doc/en/thematic_tutorials/sandpile.rst

@@ -3987,36 +3987,36 @@ EXAMPLES::
 
     sage: s = sandpiles.Complete(4)
     sage: D = SandpileDivisor(s,[4,2,0,0])
-    sage: D.effective_div()
-    [{0: 0, 1: 6, 2: 0, 3: 0},
+    sage: sorted(D.effective_div(), key=str)
+    [{0: 0, 1: 2, 2: 0, 3: 4},
      {0: 0, 1: 2, 2: 4, 3: 0},
-     {0: 0, 1: 2, 2: 0, 3: 4},
+     {0: 0, 1: 6, 2: 0, 3: 0},
      {0: 1, 1: 3, 2: 1, 3: 1},
      {0: 2, 1: 0, 2: 2, 3: 2},
      {0: 4, 1: 2, 2: 0, 3: 0}]
-    sage: D.effective_div(False)
-    [[0, 6, 0, 0],
+    sage: sorted(D.effective_div(False))
+    [[0, 2, 0, 4],
      [0, 2, 4, 0],
-     [0, 2, 0, 4],
+     [0, 6, 0, 0],
      [1, 3, 1, 1],
      [2, 0, 2, 2],
      [4, 2, 0, 0]]
-    sage: D.effective_div(with_firing_vectors=True)
-    [({0: 0, 1: 6, 2: 0, 3: 0}, (0, -2, -1, -1)),
+    sage: sorted(D.effective_div(with_firing_vectors=True), key=str)
+    [({0: 0, 1: 2, 2: 0, 3: 4}, (0, -1, -1, -2)),
      ({0: 0, 1: 2, 2: 4, 3: 0}, (0, -1, -2, -1)),
-     ({0: 0, 1: 2, 2: 0, 3: 4}, (0, -1, -1, -2)),
+     ({0: 0, 1: 6, 2: 0, 3: 0}, (0, -2, -1, -1)),
      ({0: 1, 1: 3, 2: 1, 3: 1}, (0, -1, -1, -1)),
      ({0: 2, 1: 0, 2: 2, 3: 2}, (0, 0, -1, -1)),
      ({0: 4, 1: 2, 2: 0, 3: 0}, (0, 0, 0, 0))]
-    sage: a = _[0]
+    sage: a = _[2]
     sage: a[0].values()
     [0, 6, 0, 0]
     sage: vector(D.values()) - s.laplacian()*a[1]
     (0, 6, 0, 0)
-    sage: D.effective_div(False, True)
-    [([0, 6, 0, 0], (0, -2, -1, -1)),
+    sage: sorted(D.effective_div(False, True))
+    [([0, 2, 0, 4], (0, -1, -1, -2)),
      ([0, 2, 4, 0], (0, -1, -2, -1)),
-     ([0, 2, 0, 4], (0, -1, -1, -2)),
+     ([0, 6, 0, 0], (0, -2, -1, -1)),
      ([1, 3, 1, 1], (0, -1, -1, -1)),
      ([2, 0, 2, 2], (0, 0, -1, -1)),
      ([4, 2, 0, 0], (0, 0, 0, 0))]
@@ -4377,13 +4377,13 @@ EXAMPLES::
 
     sage: s = sandpiles.Complete(4)
     sage: D = SandpileDivisor(s,[4,2,0,0])
-    sage: D.polytope_integer_pts()
-    ((-2, -1, -1),
+    sage: sorted(D.polytope_integer_pts())
+    [(-2, -1, -1),
      (-1, -2, -1),
      (-1, -1, -2),
      (-1, -1, -1),
      (0, -1, -1),
-     (0, 0, 0))
+     (0, 0, 0)]
     sage: D = SandpileDivisor(s,[-1,0,0,0])
     sage: D.polytope_integer_pts()
     ()

src/sage/geometry/polyhedron/backend_normaliz.py

@@ -566,7 +566,7 @@ def integral_points(self, threshold=10000):
                 return ()
         # for small bounding boxes, it is faster to naively iterate over the points of the box
         if threshold > 1:
-            box_min, box_max = self.bounding_box(integral=True)
+            box_min, box_max = self.bounding_box(integral_hull=True)
             box_points = prod(max_coord-min_coord+1 for min_coord, max_coord in zip(box_min, box_max))
             if  box_points<threshold:
                 from sage.geometry.integral_points import rectangular_box_points

src/sage/geometry/polyhedron/base.py

@@ -4651,7 +4651,7 @@ def _integral_points_PALP(self):
         return [p for p in lp.points() if self.contains(p)]
 
     @cached_method
-    def bounding_box(self, integral=False):
+    def bounding_box(self, integral=False, integral_hull=False):
         r"""
         Return the coordinates of a rectangular box containing the non-empty polytope.
 
@@ -4660,6 +4660,10 @@ def bounding_box(self, integral=False):
         - ``integral`` -- Boolean (default: ``False``). Whether to
           only allow integral coordinates in the bounding box.
 
+        - ``integral_hull`` -- Boolean (default: ``False``). If ``True``, return a
+          box containing the integral points of the polytope, or ``None, None`` if it
+          is known that the polytope has no integral points.
+
         OUTPUT:
 
         A pair of tuples ``(box_min, box_max)`` where ``box_min`` are
@@ -4673,6 +4677,10 @@ def bounding_box(self, integral=False):
             ((1/3, 1/3), (2/3, 2/3))
             sage: Polyhedron([ (1/3,2/3), (2/3, 1/3) ]).bounding_box(integral=True)
             ((0, 0), (1, 1))
+            sage: Polyhedron([ (1/3,2/3), (2/3, 1/3) ]).bounding_box(integral_hull=True)
+            (None, None)
+            sage: Polyhedron([ (1/3,2/3), (3/3, 4/3) ]).bounding_box(integral_hull=True)
+            ((1, 1), (1, 1))
             sage: polytopes.buckyball(exact=False).bounding_box()
             ((-0.8090169944, -0.8090169944, -0.8090169944), (0.8090169944, 0.8090169944, 0.8090169944))
         """
@@ -4686,7 +4694,14 @@ def bounding_box(self, integral=False):
             coords = [ v[i] for v in self.vertex_generator() ]
             max_coord = max(coords)
             min_coord = min(coords)
-            if integral:
+            if integral_hull:
+                a = ceil(min_coord)
+                b = floor(max_coord)
+                if a > b:
+                    return None, None
+                box_max.append(b)
+                box_min.append(a)
+            elif integral:
                 box_max.append(ceil(max_coord))
                 box_min.append(floor(min_coord))
             else:
@@ -4816,6 +4831,12 @@ def integral_points(self, threshold=100000):
             sage: len(simplex.integral_points())
             49
 
+        A case where rounding in the right direction goes a long way::
+
+            sage: P = 1/10*polytopes.hypercube(14)
+            sage: P.integral_points()
+            ((0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),)
+
         Finally, the 3-d reflexive polytope number 4078::
 
             sage: v = [(1,0,0), (0,1,0), (0,0,1), (0,0,-1), (0,-2,1),
@@ -4867,7 +4888,9 @@ def integral_points(self, threshold=100000):
                 return ()
 
         # for small bounding boxes, it is faster to naively iterate over the points of the box
-        box_min, box_max = self.bounding_box(integral=True)
+        box_min, box_max = self.bounding_box(integral_hull=True)
+        if box_min is None:
+            return ()
         box_points = prod(max_coord-min_coord+1 for min_coord, max_coord in zip(box_min, box_max))
         if  not self.is_lattice_polytope() or \
                 (self.is_simplex() and box_points < 1000) or \

src/sage/sandpiles/sandpile.py

@@ -280,14 +280,14 @@
     {0: 4, 1: 0, 2: 0, 3: 1}
     sage: D.rank()
     2
-    sage: D.effective_div()
+    sage: sorted(D.effective_div(), key=str)
     [{0: 0, 1: 0, 2: 0, 3: 5},
-     {0: 0, 1: 4, 2: 0, 3: 1},
      {0: 0, 1: 0, 2: 4, 3: 1},
+     {0: 0, 1: 4, 2: 0, 3: 1},
      {0: 1, 1: 1, 2: 1, 3: 2},
      {0: 4, 1: 0, 2: 0, 3: 1}]
-    sage: D.effective_div(False)
-    [[0, 0, 0, 5], [0, 4, 0, 1], [0, 0, 4, 1], [1, 1, 1, 2], [4, 0, 0, 1]]
+    sage: sorted(D.effective_div(False))
+    [[0, 0, 0, 5], [0, 0, 4, 1], [0, 4, 0, 1], [1, 1, 1, 2], [4, 0, 0, 1]]
     sage: D.rank()
     2
     sage: D.rank(True)
@@ -5345,13 +5345,13 @@ def polytope_integer_pts(self):
 
             sage: s = sandpiles.Complete(4)
             sage: D = SandpileDivisor(s,[4,2,0,0])
-            sage: D.polytope_integer_pts()
-            ((-2, -1, -1),
+            sage: sorted(D.polytope_integer_pts())
+            [(-2, -1, -1),
              (-1, -2, -1),
              (-1, -1, -2),
              (-1, -1, -1),
              (0, -1, -1),
-             (0, 0, 0))
+             (0, 0, 0)]
             sage: D = SandpileDivisor(s,[-1,0,0,0])
             sage: D.polytope_integer_pts()
             ()
@@ -5398,36 +5398,36 @@ def effective_div(self, verbose=True, with_firing_vectors=False):
 
             sage: s = sandpiles.Complete(4)
             sage: D = SandpileDivisor(s,[4,2,0,0])
-            sage: D.effective_div()
-            [{0: 0, 1: 6, 2: 0, 3: 0},
+            sage: sorted(D.effective_div(), key=str)
+            [{0: 0, 1: 2, 2: 0, 3: 4},
              {0: 0, 1: 2, 2: 4, 3: 0},
-             {0: 0, 1: 2, 2: 0, 3: 4},
+             {0: 0, 1: 6, 2: 0, 3: 0},
              {0: 1, 1: 3, 2: 1, 3: 1},
              {0: 2, 1: 0, 2: 2, 3: 2},
              {0: 4, 1: 2, 2: 0, 3: 0}]
-            sage: D.effective_div(False)
-            [[0, 6, 0, 0],
+            sage: sorted(D.effective_div(False))
+            [[0, 2, 0, 4],
              [0, 2, 4, 0],
-             [0, 2, 0, 4],
+             [0, 6, 0, 0],
              [1, 3, 1, 1],
              [2, 0, 2, 2],
              [4, 2, 0, 0]]
-            sage: D.effective_div(with_firing_vectors=True)
-            [({0: 0, 1: 6, 2: 0, 3: 0}, (0, -2, -1, -1)),
+            sage: sorted(D.effective_div(with_firing_vectors=True), key=str)
+            [({0: 0, 1: 2, 2: 0, 3: 4}, (0, -1, -1, -2)),
              ({0: 0, 1: 2, 2: 4, 3: 0}, (0, -1, -2, -1)),
-             ({0: 0, 1: 2, 2: 0, 3: 4}, (0, -1, -1, -2)),
+             ({0: 0, 1: 6, 2: 0, 3: 0}, (0, -2, -1, -1)),
              ({0: 1, 1: 3, 2: 1, 3: 1}, (0, -1, -1, -1)),
              ({0: 2, 1: 0, 2: 2, 3: 2}, (0, 0, -1, -1)),
              ({0: 4, 1: 2, 2: 0, 3: 0}, (0, 0, 0, 0))]
-            sage: a = _[0]
+            sage: a = _[2]
             sage: a[0].values()
             [0, 6, 0, 0]
             sage: vector(D.values()) - s.laplacian()*a[1]
             (0, 6, 0, 0)
-            sage: D.effective_div(False, True)
-            [([0, 6, 0, 0], (0, -2, -1, -1)),
+            sage: sorted(D.effective_div(False, True))
+            [([0, 2, 0, 4], (0, -1, -1, -2)),
              ([0, 2, 4, 0], (0, -1, -2, -1)),
-             ([0, 2, 0, 4], (0, -1, -1, -2)),
+             ([0, 6, 0, 0], (0, -2, -1, -1)),
              ([1, 3, 1, 1], (0, -1, -1, -1)),
              ([2, 0, 2, 2], (0, 0, -1, -1)),
              ([4, 2, 0, 0], (0, 0, 0, 0))]