Cython helper methods to compute integral points in polyhedra¶

class sage.geometry.integral_points.InequalityCollection[source]¶

Bases: object

A collection of inequalities.

INPUT:

polyhedron – a polyhedron defining the inequalities
permutation – list; a 0-based permutation of the coordinates Will be used to permute the coordinates of the inequality
box_min, box_max – the (not permuted) minimal and maximal coordinates of the bounding box; used for bounds checking

EXAMPLES:

Sage

sage: from sage.geometry.integral_points import InequalityCollection
sage: P_QQ = Polyhedron(identity_matrix(3).columns() + [(-2, -1,-1)], base_ring=QQ)
sage: ieq = InequalityCollection(P_QQ, [0,1,2], [0]*3,[1]*3); ieq
The collection of inequalities
integer: (3, -2, -2) x + 2 >= 0
integer: (-1, 4, -1) x + 1 >= 0
integer: (-1, -1, 4) x + 1 >= 0
integer: (-1, -1, -1) x + 1 >= 0

sage: P_RR = Polyhedron(identity_matrix(2).columns() + [(-2.7, -1)], base_ring=RDF)
sage: InequalityCollection(P_RR, [0,1], [0]*2, [1]*2)
The collection of inequalities
integer: (-1, -1) x + 1 >= 0
generic: (-1.0, 3.7) x + 1.0 >= 0
generic: (1.0, -1.35) x + 1.35 >= 0

sage: line = Polyhedron(eqns=[(2,3,7)])
sage: InequalityCollection(line, [0,1], [0]*2, [1]*2 )
The collection of inequalities
integer: (3, 7) x + 2 >= 0
integer: (-3, -7) x + -2 >= 0

Python

>>> from sage.all import *
>>> from sage.geometry.integral_points import InequalityCollection
>>> P_QQ = Polyhedron(identity_matrix(Integer(3)).columns() + [(-Integer(2), -Integer(1),-Integer(1))], base_ring=QQ)
>>> ieq = InequalityCollection(P_QQ, [Integer(0),Integer(1),Integer(2)], [Integer(0)]*Integer(3),[Integer(1)]*Integer(3)); ieq
The collection of inequalities
integer: (3, -2, -2) x + 2 >= 0
integer: (-1, 4, -1) x + 1 >= 0
integer: (-1, -1, 4) x + 1 >= 0
integer: (-1, -1, -1) x + 1 >= 0

>>> P_RR = Polyhedron(identity_matrix(Integer(2)).columns() + [(-RealNumber('2.7'), -Integer(1))], base_ring=RDF)
>>> InequalityCollection(P_RR, [Integer(0),Integer(1)], [Integer(0)]*Integer(2), [Integer(1)]*Integer(2))
The collection of inequalities
integer: (-1, -1) x + 1 >= 0
generic: (-1.0, 3.7) x + 1.0 >= 0
generic: (1.0, -1.35) x + 1.35 >= 0

>>> line = Polyhedron(eqns=[(Integer(2),Integer(3),Integer(7))])
>>> InequalityCollection(line, [Integer(0),Integer(1)], [Integer(0)]*Integer(2), [Integer(1)]*Integer(2) )
The collection of inequalities
integer: (3, 7) x + 2 >= 0
integer: (-3, -7) x + -2 >= 0

Sage Live

from sage.geometry.integral_points import InequalityCollection
P_QQ = Polyhedron(identity_matrix(3).columns() + [(-2, -1,-1)], base_ring=QQ)
ieq = InequalityCollection(P_QQ, [0,1,2], [0]*3,[1]*3); ieq
P_RR = Polyhedron(identity_matrix(2).columns() + [(-2.7, -1)], base_ring=RDF)
InequalityCollection(P_RR, [0,1], [0]*2, [1]*2)
line = Polyhedron(eqns=[(2,3,7)])
InequalityCollection(line, [0,1], [0]*2, [1]*2 )

are_satisfied(inner_loop_variable)[source]¶

Return whether all inequalities are satisfied.

You must call prepare_inner_loop() before calling this method.

INPUT:

inner_loop_variable – integer; the 0th coordinate of the lattice point

OUTPUT: boolean; whether the lattice point is in the polyhedron

EXAMPLES:

Sage

sage: from sage.geometry.integral_points import InequalityCollection
sage: line = Polyhedron(eqns=[(2,3,7)])
sage: ieq = InequalityCollection(line, [0,1], [0]*2, [1]*2 )
sage: ieq.prepare_next_to_inner_loop([3,4])
sage: ieq.prepare_inner_loop([3,4])
sage: ieq.are_satisfied(3)
False

Python

>>> from sage.all import *
>>> from sage.geometry.integral_points import InequalityCollection
>>> line = Polyhedron(eqns=[(Integer(2),Integer(3),Integer(7))])
>>> ieq = InequalityCollection(line, [Integer(0),Integer(1)], [Integer(0)]*Integer(2), [Integer(1)]*Integer(2) )
>>> ieq.prepare_next_to_inner_loop([Integer(3),Integer(4)])
>>> ieq.prepare_inner_loop([Integer(3),Integer(4)])
>>> ieq.are_satisfied(Integer(3))
False

Sage Live

from sage.geometry.integral_points import InequalityCollection
line = Polyhedron(eqns=[(2,3,7)])
ieq = InequalityCollection(line, [0,1], [0]*2, [1]*2 )
ieq.prepare_next_to_inner_loop([3,4])
ieq.prepare_inner_loop([3,4])
ieq.are_satisfied(3)

prepare_inner_loop(p)[source]¶

Peel off the inner loop.

In the inner loop of rectangular_box_points(), we have to repeatedly evaluate $A x + b \geq 0$ . To speed up computation, we pre-evaluate

c = A x - A_{0} x_{0} + b = b + \sum_{i = 1} A_{i} x_{i}

and only test $A_{0} x_{0} + c \geq 0$ in the inner loop.

You must call prepare_next_to_inner_loop() before calling this method.

INPUT:

p – the coordinates of the point to loop over. Only the p[1:] entries are used

EXAMPLES:

Sage

sage: from sage.geometry.integral_points import InequalityCollection, print_cache
sage: P = Polyhedron(ieqs=[(2,3,7,11)])
sage: ieq = InequalityCollection(P, [0,1,2], [0]*3,[1]*3); ieq
The collection of inequalities
integer: (3, 7, 11) x + 2 >= 0
sage: ieq.prepare_next_to_inner_loop([2,1,3])
sage: ieq.prepare_inner_loop([2,1,3])
sage: print_cache(ieq)
Cached inner loop: 3 * x_0 + 42 >= 0
Cached next-to-inner loop: 3 * x_0 + 7 * x_1 + 35 >= 0

Python

>>> from sage.all import *
>>> from sage.geometry.integral_points import InequalityCollection, print_cache
>>> P = Polyhedron(ieqs=[(Integer(2),Integer(3),Integer(7),Integer(11))])
>>> ieq = InequalityCollection(P, [Integer(0),Integer(1),Integer(2)], [Integer(0)]*Integer(3),[Integer(1)]*Integer(3)); ieq
The collection of inequalities
integer: (3, 7, 11) x + 2 >= 0
>>> ieq.prepare_next_to_inner_loop([Integer(2),Integer(1),Integer(3)])
>>> ieq.prepare_inner_loop([Integer(2),Integer(1),Integer(3)])
>>> print_cache(ieq)
Cached inner loop: 3 * x_0 + 42 >= 0
Cached next-to-inner loop: 3 * x_0 + 7 * x_1 + 35 >= 0

Sage Live

from sage.geometry.integral_points import InequalityCollection, print_cache
P = Polyhedron(ieqs=[(2,3,7,11)])
ieq = InequalityCollection(P, [0,1,2], [0]*3,[1]*3); ieq
ieq.prepare_next_to_inner_loop([2,1,3])
ieq.prepare_inner_loop([2,1,3])
print_cache(ieq)

prepare_next_to_inner_loop(p)[source]¶

Peel off the next-to-inner loop.

In the next-to-inner loop of rectangular_box_points(), we have to repeatedly evaluate $A x - A_{0} x_{0} + b$ . To speed up computation, we pre-evaluate

c = b + \sum_{i = 2} A_{i} x_{i}

and only compute $A x - A_{0} x_{0} + b = A_{1} x_{1} + c \geq 0$ in the next-to-inner loop.

INPUT:

p – the point coordinates. Only p[2:] coordinates are potentially used by this method

EXAMPLES:

Sage

sage: from sage.geometry.integral_points import InequalityCollection, print_cache
sage: P = Polyhedron(ieqs=[(2,3,7,11)])
sage: ieq = InequalityCollection(P, [0,1,2], [0]*3,[1]*3); ieq
The collection of inequalities
integer: (3, 7, 11) x + 2 >= 0
sage: ieq.prepare_next_to_inner_loop([2,1,3])
sage: ieq.prepare_inner_loop([2,1,3])
sage: print_cache(ieq)
Cached inner loop: 3 * x_0 + 42 >= 0
Cached next-to-inner loop: 3 * x_0 + 7 * x_1 + 35 >= 0

Python

>>> from sage.all import *
>>> from sage.geometry.integral_points import InequalityCollection, print_cache
>>> P = Polyhedron(ieqs=[(Integer(2),Integer(3),Integer(7),Integer(11))])
>>> ieq = InequalityCollection(P, [Integer(0),Integer(1),Integer(2)], [Integer(0)]*Integer(3),[Integer(1)]*Integer(3)); ieq
The collection of inequalities
integer: (3, 7, 11) x + 2 >= 0
>>> ieq.prepare_next_to_inner_loop([Integer(2),Integer(1),Integer(3)])
>>> ieq.prepare_inner_loop([Integer(2),Integer(1),Integer(3)])
>>> print_cache(ieq)
Cached inner loop: 3 * x_0 + 42 >= 0
Cached next-to-inner loop: 3 * x_0 + 7 * x_1 + 35 >= 0

Sage Live

from sage.geometry.integral_points import InequalityCollection, print_cache
P = Polyhedron(ieqs=[(2,3,7,11)])
ieq = InequalityCollection(P, [0,1,2], [0]*3,[1]*3); ieq
ieq.prepare_next_to_inner_loop([2,1,3])
ieq.prepare_inner_loop([2,1,3])
print_cache(ieq)

satisfied_as_equalities(inner_loop_variable)[source]¶

Return the inequalities (by their index) that are satisfied as equalities.

INPUT:

inner_loop_variable – integer; the 0th coordinate of the lattice point

OUTPUT:

A set of integers in ascending order. Each integer is the index of a H-representation object of the polyhedron (either a inequality or an equation).

EXAMPLES:

Sage

sage: from sage.geometry.integral_points import InequalityCollection
sage: quadrant = Polyhedron(rays=[(1,0), (0,1)])
sage: ieqs = InequalityCollection(quadrant, [0,1], [-1]*2, [1]*2 )
sage: ieqs.prepare_next_to_inner_loop([-1,0])
sage: ieqs.prepare_inner_loop([-1,0])
sage: ieqs.satisfied_as_equalities(-1)
frozenset({1})
sage: ieqs.satisfied_as_equalities(0)
frozenset({0, 1})
sage: ieqs.satisfied_as_equalities(1)
frozenset({1})

Python

>>> from sage.all import *
>>> from sage.geometry.integral_points import InequalityCollection
>>> quadrant = Polyhedron(rays=[(Integer(1),Integer(0)), (Integer(0),Integer(1))])
>>> ieqs = InequalityCollection(quadrant, [Integer(0),Integer(1)], [-Integer(1)]*Integer(2), [Integer(1)]*Integer(2) )
>>> ieqs.prepare_next_to_inner_loop([-Integer(1),Integer(0)])
>>> ieqs.prepare_inner_loop([-Integer(1),Integer(0)])
>>> ieqs.satisfied_as_equalities(-Integer(1))
frozenset({1})
>>> ieqs.satisfied_as_equalities(Integer(0))
frozenset({0, 1})
>>> ieqs.satisfied_as_equalities(Integer(1))
frozenset({1})

Sage Live

from sage.geometry.integral_points import InequalityCollection
quadrant = Polyhedron(rays=[(1,0), (0,1)])
ieqs = InequalityCollection(quadrant, [0,1], [-1]*2, [1]*2 )
ieqs.prepare_next_to_inner_loop([-1,0])
ieqs.prepare_inner_loop([-1,0])
ieqs.satisfied_as_equalities(-1)
ieqs.satisfied_as_equalities(0)
ieqs.satisfied_as_equalities(1)

swap_ineq_to_front(i)[source]¶

Swap the i-th entry of the list to the front of the list of inequalities.

INPUT:

i – integer; the Inequality_int to swap to the beginning of the list of integral inequalities

EXAMPLES:

Sage

sage: from sage.geometry.integral_points import InequalityCollection
sage: P_QQ = Polyhedron(identity_matrix(3).columns() + [(-2, -1,-1)], base_ring=QQ)
sage: iec = InequalityCollection(P_QQ, [0,1,2], [0]*3,[1]*3)
sage: iec
The collection of inequalities
integer: (3, -2, -2) x + 2 >= 0
integer: (-1, 4, -1) x + 1 >= 0
integer: (-1, -1, 4) x + 1 >= 0
integer: (-1, -1, -1) x + 1 >= 0
sage: iec.swap_ineq_to_front(3)
sage: iec
The collection of inequalities
integer: (-1, -1, -1) x + 1 >= 0
integer: (3, -2, -2) x + 2 >= 0
integer: (-1, 4, -1) x + 1 >= 0
integer: (-1, -1, 4) x + 1 >= 0

Python

>>> from sage.all import *
>>> from sage.geometry.integral_points import InequalityCollection
>>> P_QQ = Polyhedron(identity_matrix(Integer(3)).columns() + [(-Integer(2), -Integer(1),-Integer(1))], base_ring=QQ)
>>> iec = InequalityCollection(P_QQ, [Integer(0),Integer(1),Integer(2)], [Integer(0)]*Integer(3),[Integer(1)]*Integer(3))
>>> iec
The collection of inequalities
integer: (3, -2, -2) x + 2 >= 0
integer: (-1, 4, -1) x + 1 >= 0
integer: (-1, -1, 4) x + 1 >= 0
integer: (-1, -1, -1) x + 1 >= 0
>>> iec.swap_ineq_to_front(Integer(3))
>>> iec
The collection of inequalities
integer: (-1, -1, -1) x + 1 >= 0
integer: (3, -2, -2) x + 2 >= 0
integer: (-1, 4, -1) x + 1 >= 0
integer: (-1, -1, 4) x + 1 >= 0

Sage Live

from sage.geometry.integral_points import InequalityCollection
P_QQ = Polyhedron(identity_matrix(3).columns() + [(-2, -1,-1)], base_ring=QQ)
iec = InequalityCollection(P_QQ, [0,1,2], [0]*3,[1]*3)
iec
iec.swap_ineq_to_front(3)
iec

class sage.geometry.integral_points.Inequality_generic[source]¶

Bases: object

An inequality whose coefficients are arbitrary Python/Sage objects

INPUT:

A – list of coefficients
b – element

OUTPUT: inequality $A x + b \geq 0$

EXAMPLES:

Sage

sage: from sage.geometry.integral_points import Inequality_generic
sage: Inequality_generic([2 * pi, sqrt(3), 7/2], -5.5)                          # needs sage.symbolic
generic: (2*pi, sqrt(3), 7/2) x + -5.50000000000000 >= 0

Python

>>> from sage.all import *
>>> from sage.geometry.integral_points import Inequality_generic
>>> Inequality_generic([Integer(2) * pi, sqrt(Integer(3)), Integer(7)/Integer(2)], -RealNumber('5.5'))                          # needs sage.symbolic
generic: (2*pi, sqrt(3), 7/2) x + -5.50000000000000 >= 0

Sage Live

from sage.geometry.integral_points import Inequality_generic
Inequality_generic([2 * pi, sqrt(3), 7/2], -5.5)                          # needs sage.symbolic

class sage.geometry.integral_points.Inequality_int[source]¶

Bases: object

Fast version of inequality in the case that all coefficients fit into machine ints.

INPUT:

A – list of integers
b – integer
max_abs_coordinates – the maximum of the coordinates that one wants to evaluate the coordinates on; used for overflow checking

OUTPUT:

Inequality $A x + b \geq 0$ . A OverflowError is raised if a machine integer is not long enough to hold the results. A ValueError is raised if some of the input is not integral.

EXAMPLES:

Sage

sage: from sage.geometry.integral_points import Inequality_int
sage: Inequality_int([2,3,7], -5, [10]*3)
integer: (2, 3, 7) x + -5 >= 0

sage: Inequality_int([1]*21, -5, [10]*21)
Traceback (most recent call last):
...
OverflowError: Dimension limit exceeded.

sage: Inequality_int([2,3/2,7], -5, [10]*3)
Traceback (most recent call last):
...
ValueError: Not integral.

sage: Inequality_int([2,3,7], -5.2, [10]*3)
Traceback (most recent call last):
...
ValueError: Not integral.

sage: Inequality_int([2,3,7], -5*10^50, [10]*3)  # actual error message can differ between 32 and 64 bit
Traceback (most recent call last):
...
OverflowError: ...

Python

>>> from sage.all import *
>>> from sage.geometry.integral_points import Inequality_int
>>> Inequality_int([Integer(2),Integer(3),Integer(7)], -Integer(5), [Integer(10)]*Integer(3))
integer: (2, 3, 7) x + -5 >= 0

>>> Inequality_int([Integer(1)]*Integer(21), -Integer(5), [Integer(10)]*Integer(21))
Traceback (most recent call last):
...
OverflowError: Dimension limit exceeded.

>>> Inequality_int([Integer(2),Integer(3)/Integer(2),Integer(7)], -Integer(5), [Integer(10)]*Integer(3))
Traceback (most recent call last):
...
ValueError: Not integral.

>>> Inequality_int([Integer(2),Integer(3),Integer(7)], -RealNumber('5.2'), [Integer(10)]*Integer(3))
Traceback (most recent call last):
...
ValueError: Not integral.

>>> Inequality_int([Integer(2),Integer(3),Integer(7)], -Integer(5)*Integer(10)**Integer(50), [Integer(10)]*Integer(3))  # actual error message can differ between 32 and 64 bit
Traceback (most recent call last):
...
OverflowError: ...

Sage Live

from sage.geometry.integral_points import Inequality_int
Inequality_int([2,3,7], -5, [10]*3)
Inequality_int([1]*21, -5, [10]*21)
Inequality_int([2,3/2,7], -5, [10]*3)
Inequality_int([2,3,7], -5.2, [10]*3)
Inequality_int([2,3,7], -5*10^50, [10]*3)  # actual error message can differ between 32 and 64 bit

sage.geometry.integral_points.loop_over_parallelotope_points(e, d, VDinv, R, lattice, A=None, b=None)[source]¶

The inner loop of parallelotope_points().

INPUT:

See parallelotope_points() for e, d, VDinv, R, lattice.

A, b – either both None or a vector and number. If present, only the parallelotope points satisfying $A x \leq b$ are returned.

OUTPUT:

The points of the half-open parallelotope as a tuple of lattice points.

EXAMPLES:

Sage

sage: e = [3]
sage: d = prod(e)
sage: VDinv = matrix(ZZ, [[1]])
sage: R = column_matrix(ZZ,[3,3,3])
sage: lattice = ZZ^3
sage: from sage.geometry.integral_points import loop_over_parallelotope_points
sage: loop_over_parallelotope_points(e, d, VDinv, R, lattice)
((0, 0, 0), (1, 1, 1), (2, 2, 2))

sage: A = vector(ZZ, [1,0,0])
sage: b = 1
sage: loop_over_parallelotope_points(e, d, VDinv, R, lattice, A, b)
((0, 0, 0), (1, 1, 1))

Python

>>> from sage.all import *
>>> e = [Integer(3)]
>>> d = prod(e)
>>> VDinv = matrix(ZZ, [[Integer(1)]])
>>> R = column_matrix(ZZ,[Integer(3),Integer(3),Integer(3)])
>>> lattice = ZZ**Integer(3)
>>> from sage.geometry.integral_points import loop_over_parallelotope_points
>>> loop_over_parallelotope_points(e, d, VDinv, R, lattice)
((0, 0, 0), (1, 1, 1), (2, 2, 2))

>>> A = vector(ZZ, [Integer(1),Integer(0),Integer(0)])
>>> b = Integer(1)
>>> loop_over_parallelotope_points(e, d, VDinv, R, lattice, A, b)
((0, 0, 0), (1, 1, 1))

Sage Live

e = [3]
d = prod(e)
VDinv = matrix(ZZ, [[1]])
R = column_matrix(ZZ,[3,3,3])
lattice = ZZ^3
from sage.geometry.integral_points import loop_over_parallelotope_points
loop_over_parallelotope_points(e, d, VDinv, R, lattice)
A = vector(ZZ, [1,0,0])
b = 1
loop_over_parallelotope_points(e, d, VDinv, R, lattice, A, b)

sage.geometry.integral_points.parallelotope_points(spanning_points, lattice)[source]¶

Return integral points in the parallelotope starting at the origin and spanned by the spanning_points.

See semigroup_generators() for a description of the algorithm.

INPUT:

spanning_points – a non-empty list of linearly independent rays ( $Z$ -vectors or toric lattice elements), not necessarily primitive lattice points.

OUTPUT:

The tuple of all lattice points in the half-open parallelotope spanned by the rays $r_{i}$ ,

p a r ({r_{i}}) = \sum_{0 \leq a_{i} < 1} a_{i} r_{i}

By half-open parallelotope, we mean that the points in the facets not meeting the origin are omitted.

EXAMPLES:

Note how the points on the outward-facing factes are omitted:

Sage

sage: from sage.geometry.integral_points import parallelotope_points
sage: rays = list(map(vector, [(2,0), (0,2)]))
sage: parallelotope_points(rays, ZZ^2)
((0, 0), (0, 1), (1, 0), (1, 1))

Python

>>> from sage.all import *
>>> from sage.geometry.integral_points import parallelotope_points
>>> rays = list(map(vector, [(Integer(2),Integer(0)), (Integer(0),Integer(2))]))
>>> parallelotope_points(rays, ZZ**Integer(2))
((0, 0), (0, 1), (1, 0), (1, 1))

Sage Live

from sage.geometry.integral_points import parallelotope_points
rays = list(map(vector, [(2,0), (0,2)]))
parallelotope_points(rays, ZZ^2)

The rays can also be toric lattice points:

Sage

sage: rays = list(map(ToricLattice(2), [(2,0), (0,2)]))
sage: parallelotope_points(rays, ToricLattice(2))
(N(0, 0), N(0, 1), N(1, 0), N(1, 1))

Python

>>> from sage.all import *
>>> rays = list(map(ToricLattice(Integer(2)), [(Integer(2),Integer(0)), (Integer(0),Integer(2))]))
>>> parallelotope_points(rays, ToricLattice(Integer(2)))
(N(0, 0), N(0, 1), N(1, 0), N(1, 1))

Sage Live

rays = list(map(ToricLattice(2), [(2,0), (0,2)]))
parallelotope_points(rays, ToricLattice(2))

A non-smooth cone:

Sage

sage: c = Cone([ (1,0), (1,2) ])
sage: parallelotope_points(c.rays(), c.lattice())
(N(0, 0), N(1, 1))

Python

>>> from sage.all import *
>>> c = Cone([ (Integer(1),Integer(0)), (Integer(1),Integer(2)) ])
>>> parallelotope_points(c.rays(), c.lattice())
(N(0, 0), N(1, 1))

Sage Live

c = Cone([ (1,0), (1,2) ])
parallelotope_points(c.rays(), c.lattice())

A ValueError is raised if the spanning_points are not linearly independent:

Sage

sage: rays = list(map(ToricLattice(2), [(1,1)]*2))
sage: parallelotope_points(rays, ToricLattice(2))
Traceback (most recent call last):
...
ValueError: The spanning points are not linearly independent!

Python

>>> from sage.all import *
>>> rays = list(map(ToricLattice(Integer(2)), [(Integer(1),Integer(1))]*Integer(2)))
>>> parallelotope_points(rays, ToricLattice(Integer(2)))
Traceback (most recent call last):
...
ValueError: The spanning points are not linearly independent!

Sage Live

rays = list(map(ToricLattice(2), [(1,1)]*2))
parallelotope_points(rays, ToricLattice(2))

sage.geometry.integral_points.print_cache(inequality_collection)[source]¶

Print the cached values in Inequality_int (for debugging/doctesting only).

EXAMPLES:

Sage

sage: from sage.geometry.integral_points import InequalityCollection, print_cache
sage: P = Polyhedron(ieqs=[(2,3,7)])
sage: ieq = InequalityCollection(P, [0,1], [0]*2,[1]*2); ieq
The collection of inequalities
integer: (3, 7) x + 2 >= 0
sage: ieq.prepare_next_to_inner_loop([3,5])
sage: ieq.prepare_inner_loop([3,5])
sage: print_cache(ieq)
Cached inner loop: 3 * x_0 + 37 >= 0
Cached next-to-inner loop: 3 * x_0 + 7 * x_1 + 2 >= 0

Python

>>> from sage.all import *
>>> from sage.geometry.integral_points import InequalityCollection, print_cache
>>> P = Polyhedron(ieqs=[(Integer(2),Integer(3),Integer(7))])
>>> ieq = InequalityCollection(P, [Integer(0),Integer(1)], [Integer(0)]*Integer(2),[Integer(1)]*Integer(2)); ieq
The collection of inequalities
integer: (3, 7) x + 2 >= 0
>>> ieq.prepare_next_to_inner_loop([Integer(3),Integer(5)])
>>> ieq.prepare_inner_loop([Integer(3),Integer(5)])
>>> print_cache(ieq)
Cached inner loop: 3 * x_0 + 37 >= 0
Cached next-to-inner loop: 3 * x_0 + 7 * x_1 + 2 >= 0

Sage Live

from sage.geometry.integral_points import InequalityCollection, print_cache
P = Polyhedron(ieqs=[(2,3,7)])
ieq = InequalityCollection(P, [0,1], [0]*2,[1]*2); ieq
ieq.prepare_next_to_inner_loop([3,5])
ieq.prepare_inner_loop([3,5])
print_cache(ieq)

sage.geometry.integral_points.ray_matrix_normal_form(R)[source]¶

Compute the Smith normal form of the ray matrix for parallelotope_points().

INPUT:

R – $Z$ -matrix whose columns are the rays spanning the parallelotope

OUTPUT: a tuple containing e, d, and VDinv

EXAMPLES:

Sage

sage: from sage.geometry.integral_points import ray_matrix_normal_form
sage: R = column_matrix(ZZ,[3,3,3])
sage: ray_matrix_normal_form(R)
([3], 3, [1])

Python

>>> from sage.all import *
>>> from sage.geometry.integral_points import ray_matrix_normal_form
>>> R = column_matrix(ZZ,[Integer(3),Integer(3),Integer(3)])
>>> ray_matrix_normal_form(R)
([3], 3, [1])

Sage Live

from sage.geometry.integral_points import ray_matrix_normal_form
R = column_matrix(ZZ,[3,3,3])
ray_matrix_normal_form(R)

sage.geometry.integral_points.rectangular_box_points(box_min, box_max, polyhedron=None, count_only=False, return_saturated=False)[source]¶

Return the integral points in the lattice bounding box that are also contained in the given polyhedron.

INPUT:

box_min – list of integers; the minimal value for each coordinate of the rectangular bounding box
box_max – list of integers; the maximal value for each coordinate of the rectangular bounding box
polyhedron – a Polyhedron_base, a PPL C_Polyhedron, or None (default)
count_only – boolean (default: False); whether to return only the total number of vertices, and not their coordinates. Enabling this option speeds up the enumeration. Cannot be combined with the return_saturated option.
return_saturated – boolean (default: False); whether to also return which inequalities are saturated for each point of the polyhedron. Enabling this slows down the enumeration. Cannot be combined with the count_only option.

OUTPUT:

By default, this function returns a tuple containing the integral points of the rectangular box spanned by box_min and box_max and that lie inside the polyhedron. For sufficiently large bounding boxes, this are all integral points of the polyhedron.

If no polyhedron is specified, all integral points of the rectangular box are returned.

If count_only is specified, only the total number (an integer) of found lattice points is returned.

If return_saturated is enabled, then for each integral point a pair (point, Hrep) is returned where point is the point and Hrep is the set of indices of the H-representation objects that are saturated at the point.

ALGORITHM:

This function implements the naive algorithm towards counting integral points. Given min and max of vertex coordinates, it iterates over all points in the bounding box and checks whether they lie in the polyhedron. The following optimizations are implemented:

Cython: Use machine integers and optimizing C/C++ compiler where possible, arbitrary precision integers where necessary. Bounds checking, no compile time limits.
Unwind inner loop (and next-to-inner loop):

$A x \leq b \Leftrightarrow a_{1} x_{1} \leq b - \sum_{i = 2}^{d} a_{i} x_{i}$

so we only have to evaluate $a_{1} * x_{1}$ in the inner loop.
Coordinates are permuted to make the longest box edge the inner loop. The inner loop is optimized to run very fast, so its best to do as much work as possible there.
Continuously reorder inequalities and test the most restrictive inequalities first.
Use convexity and only find first and last allowed point in the inner loop. The points in-between must be points of the polyhedron, too.

EXAMPLES:

Sage

sage: from sage.geometry.integral_points import rectangular_box_points
sage: rectangular_box_points([0,0,0],[1,2,3])
((0, 0, 0), (0, 0, 1), (0, 0, 2), (0, 0, 3),
 (0, 1, 0), (0, 1, 1), (0, 1, 2), (0, 1, 3),
 (0, 2, 0), (0, 2, 1), (0, 2, 2), (0, 2, 3),
 (1, 0, 0), (1, 0, 1), (1, 0, 2), (1, 0, 3),
 (1, 1, 0), (1, 1, 1), (1, 1, 2), (1, 1, 3),
 (1, 2, 0), (1, 2, 1), (1, 2, 2), (1, 2, 3))

sage: from sage.geometry.integral_points import rectangular_box_points
sage: rectangular_box_points([0,0,0],[1,2,3], count_only=True)
24

sage: cell24 = polytopes.twenty_four_cell()
sage: rectangular_box_points([-1]*4, [1]*4, cell24)
((-1, 0, 0, 0), (0, -1, 0, 0), (0, 0, -1, 0), (0, 0, 0, -1),
 (0, 0, 0, 0),
 (0, 0, 0, 1), (0, 0, 1, 0), (0, 1, 0, 0), (1, 0, 0, 0))
sage: d = 3
sage: dilated_cell24 = d*cell24
sage: len( rectangular_box_points([-d]*4, [d]*4, dilated_cell24) )
305

sage: d = 6
sage: dilated_cell24 = d*cell24
sage: len( rectangular_box_points([-d]*4, [d]*4, dilated_cell24) )
3625

sage: rectangular_box_points([-d]*4, [d]*4, dilated_cell24, count_only=True)
3625

sage: polytope = Polyhedron([(-4,-3,-2,-1),(3,1,1,1),(1,2,1,1),(1,1,3,0),(1,3,2,4)])
sage: pts = rectangular_box_points([-4]*4, [4]*4, polytope); pts
((-4, -3, -2, -1), (-1, 0, 0, 1), (0, 1, 1, 1), (1, 1, 1, 1), (1, 1, 3, 0),
 (1, 2, 1, 1), (1, 2, 2, 2), (1, 3, 2, 4), (2, 1, 1, 1), (3, 1, 1, 1))
sage: all(polytope.contains(p) for p in pts)
True

sage: set(map(tuple,pts)) == \
....: set([(-4,-3,-2,-1),(3,1,1,1),(1,2,1,1),(1,1,3,0),(1,3,2,4),
....:      (0,1,1,1),(1,2,2,2),(-1,0,0,1),(1,1,1,1),(2,1,1,1)])   # computed with PALP
True

Python

>>> from sage.all import *
>>> from sage.geometry.integral_points import rectangular_box_points
>>> rectangular_box_points([Integer(0),Integer(0),Integer(0)],[Integer(1),Integer(2),Integer(3)])
((0, 0, 0), (0, 0, 1), (0, 0, 2), (0, 0, 3),
 (0, 1, 0), (0, 1, 1), (0, 1, 2), (0, 1, 3),
 (0, 2, 0), (0, 2, 1), (0, 2, 2), (0, 2, 3),
 (1, 0, 0), (1, 0, 1), (1, 0, 2), (1, 0, 3),
 (1, 1, 0), (1, 1, 1), (1, 1, 2), (1, 1, 3),
 (1, 2, 0), (1, 2, 1), (1, 2, 2), (1, 2, 3))

>>> from sage.geometry.integral_points import rectangular_box_points
>>> rectangular_box_points([Integer(0),Integer(0),Integer(0)],[Integer(1),Integer(2),Integer(3)], count_only=True)
24

>>> cell24 = polytopes.twenty_four_cell()
>>> rectangular_box_points([-Integer(1)]*Integer(4), [Integer(1)]*Integer(4), cell24)
((-1, 0, 0, 0), (0, -1, 0, 0), (0, 0, -1, 0), (0, 0, 0, -1),
 (0, 0, 0, 0),
 (0, 0, 0, 1), (0, 0, 1, 0), (0, 1, 0, 0), (1, 0, 0, 0))
>>> d = Integer(3)
>>> dilated_cell24 = d*cell24
>>> len( rectangular_box_points([-d]*Integer(4), [d]*Integer(4), dilated_cell24) )
305

>>> d = Integer(6)
>>> dilated_cell24 = d*cell24
>>> len( rectangular_box_points([-d]*Integer(4), [d]*Integer(4), dilated_cell24) )
3625

>>> rectangular_box_points([-d]*Integer(4), [d]*Integer(4), dilated_cell24, count_only=True)
3625

>>> polytope = Polyhedron([(-Integer(4),-Integer(3),-Integer(2),-Integer(1)),(Integer(3),Integer(1),Integer(1),Integer(1)),(Integer(1),Integer(2),Integer(1),Integer(1)),(Integer(1),Integer(1),Integer(3),Integer(0)),(Integer(1),Integer(3),Integer(2),Integer(4))])
>>> pts = rectangular_box_points([-Integer(4)]*Integer(4), [Integer(4)]*Integer(4), polytope); pts
((-4, -3, -2, -1), (-1, 0, 0, 1), (0, 1, 1, 1), (1, 1, 1, 1), (1, 1, 3, 0),
 (1, 2, 1, 1), (1, 2, 2, 2), (1, 3, 2, 4), (2, 1, 1, 1), (3, 1, 1, 1))
>>> all(polytope.contains(p) for p in pts)
True

>>> set(map(tuple,pts)) == set([(-Integer(4),-Integer(3),-Integer(2),-Integer(1)),(Integer(3),Integer(1),Integer(1),Integer(1)),(Integer(1),Integer(2),Integer(1),Integer(1)),(Integer(1),Integer(1),Integer(3),Integer(0)),(Integer(1),Integer(3),Integer(2),Integer(4)),
...      (Integer(0),Integer(1),Integer(1),Integer(1)),(Integer(1),Integer(2),Integer(2),Integer(2)),(-Integer(1),Integer(0),Integer(0),Integer(1)),(Integer(1),Integer(1),Integer(1),Integer(1)),(Integer(2),Integer(1),Integer(1),Integer(1))])   # computed with PALP
True

Sage Live

from sage.geometry.integral_points import rectangular_box_points
rectangular_box_points([0,0,0],[1,2,3])
from sage.geometry.integral_points import rectangular_box_points
rectangular_box_points([0,0,0],[1,2,3], count_only=True)
cell24 = polytopes.twenty_four_cell()
rectangular_box_points([-1]*4, [1]*4, cell24)
d = 3
dilated_cell24 = d*cell24
len( rectangular_box_points([-d]*4, [d]*4, dilated_cell24) )
d = 6
dilated_cell24 = d*cell24
len( rectangular_box_points([-d]*4, [d]*4, dilated_cell24) )
rectangular_box_points([-d]*4, [d]*4, dilated_cell24, count_only=True)
polytope = Polyhedron([(-4,-3,-2,-1),(3,1,1,1),(1,2,1,1),(1,1,3,0),(1,3,2,4)])
pts = rectangular_box_points([-4]*4, [4]*4, polytope); pts
all(polytope.contains(p) for p in pts)
set(map(tuple,pts)) == \
set([(-4,-3,-2,-1),(3,1,1,1),(1,2,1,1),(1,1,3,0),(1,3,2,4),
     (0,1,1,1),(1,2,2,2),(-1,0,0,1),(1,1,1,1),(2,1,1,1)])   # computed with PALP

Long ints and non-integral polyhedra are explicitly allowed:

Sage

sage: polytope = Polyhedron([[1], [10*pi.n()]], base_ring=RDF)                  # needs sage.symbolic
sage: len(rectangular_box_points([-100], [100], polytope))                      # needs sage.symbolic
31

sage: halfplane = Polyhedron(ieqs=[(-1,1,0)])
sage: rectangular_box_points([0,-1+10^50], [0,1+10^50])
((0, 99999999999999999999999999999999999999999999999999),
 (0, 100000000000000000000000000000000000000000000000000),
 (0, 100000000000000000000000000000000000000000000000001))
sage: len( rectangular_box_points([0,-100+10^50], [1,100+10^50], halfplane) )
201

Python

>>> from sage.all import *
>>> polytope = Polyhedron([[Integer(1)], [Integer(10)*pi.n()]], base_ring=RDF)                  # needs sage.symbolic
>>> len(rectangular_box_points([-Integer(100)], [Integer(100)], polytope))                      # needs sage.symbolic
31

>>> halfplane = Polyhedron(ieqs=[(-Integer(1),Integer(1),Integer(0))])
>>> rectangular_box_points([Integer(0),-Integer(1)+Integer(10)**Integer(50)], [Integer(0),Integer(1)+Integer(10)**Integer(50)])
((0, 99999999999999999999999999999999999999999999999999),
 (0, 100000000000000000000000000000000000000000000000000),
 (0, 100000000000000000000000000000000000000000000000001))
>>> len( rectangular_box_points([Integer(0),-Integer(100)+Integer(10)**Integer(50)], [Integer(1),Integer(100)+Integer(10)**Integer(50)], halfplane) )
201

Sage Live

polytope = Polyhedron([[1], [10*pi.n()]], base_ring=RDF)                  # needs sage.symbolic
len(rectangular_box_points([-100], [100], polytope))                      # needs sage.symbolic
halfplane = Polyhedron(ieqs=[(-1,1,0)])
rectangular_box_points([0,-1+10^50], [0,1+10^50])
len( rectangular_box_points([0,-100+10^50], [1,100+10^50], halfplane) )

Using a PPL polyhedron:

Sage

sage: # needs pplpy
sage: from ppl import Variable, Generator_System, C_Polyhedron, point
sage: gs = Generator_System()
sage: x = Variable(0); y = Variable(1); z = Variable(2)
sage: gs.insert(point(0*x + 1*y + 0*z))
sage: gs.insert(point(0*x + 1*y + 3*z))
sage: gs.insert(point(3*x + 1*y + 0*z))
sage: gs.insert(point(3*x + 1*y + 3*z))
sage: poly = C_Polyhedron(gs)
sage: rectangular_box_points([0]*3, [3]*3, poly)
((0, 1, 0), (0, 1, 1), (0, 1, 2), (0, 1, 3), (1, 1, 0), (1, 1, 1), (1, 1, 2), (1, 1, 3),
 (2, 1, 0), (2, 1, 1), (2, 1, 2), (2, 1, 3), (3, 1, 0), (3, 1, 1), (3, 1, 2), (3, 1, 3))

Python

>>> from sage.all import *
>>> # needs pplpy
>>> from ppl import Variable, Generator_System, C_Polyhedron, point
>>> gs = Generator_System()
>>> x = Variable(Integer(0)); y = Variable(Integer(1)); z = Variable(Integer(2))
>>> gs.insert(point(Integer(0)*x + Integer(1)*y + Integer(0)*z))
>>> gs.insert(point(Integer(0)*x + Integer(1)*y + Integer(3)*z))
>>> gs.insert(point(Integer(3)*x + Integer(1)*y + Integer(0)*z))
>>> gs.insert(point(Integer(3)*x + Integer(1)*y + Integer(3)*z))
>>> poly = C_Polyhedron(gs)
>>> rectangular_box_points([Integer(0)]*Integer(3), [Integer(3)]*Integer(3), poly)
((0, 1, 0), (0, 1, 1), (0, 1, 2), (0, 1, 3), (1, 1, 0), (1, 1, 1), (1, 1, 2), (1, 1, 3),
 (2, 1, 0), (2, 1, 1), (2, 1, 2), (2, 1, 3), (3, 1, 0), (3, 1, 1), (3, 1, 2), (3, 1, 3))

Sage Live

# needs pplpy
from ppl import Variable, Generator_System, C_Polyhedron, point
gs = Generator_System()
x = Variable(0); y = Variable(1); z = Variable(2)
gs.insert(point(0*x + 1*y + 0*z))
gs.insert(point(0*x + 1*y + 3*z))
gs.insert(point(3*x + 1*y + 0*z))
gs.insert(point(3*x + 1*y + 3*z))
poly = C_Polyhedron(gs)
rectangular_box_points([0]*3, [3]*3, poly)

Optionally, return the information about the saturated inequalities as well:

Sage

sage: cube = polytopes.cube()
sage: cube.Hrepresentation(0)
An inequality (-1, 0, 0) x + 1 >= 0
sage: cube.Hrepresentation(1)
An inequality (0, -1, 0) x + 1 >= 0
sage: cube.Hrepresentation(2)
An inequality (0, 0, -1) x + 1 >= 0
sage: rectangular_box_points([0]*3, [1]*3, cube, return_saturated=True)
(((0, 0, 0), frozenset()),
 ((0, 0, 1), frozenset({2})),
 ((0, 1, 0), frozenset({1})),
 ((0, 1, 1), frozenset({1, 2})),
 ((1, 0, 0), frozenset({0})),
 ((1, 0, 1), frozenset({0, 2})),
 ((1, 1, 0), frozenset({0, 1})),
 ((1, 1, 1), frozenset({0, 1, 2})))

Python

>>> from sage.all import *
>>> cube = polytopes.cube()
>>> cube.Hrepresentation(Integer(0))
An inequality (-1, 0, 0) x + 1 >= 0
>>> cube.Hrepresentation(Integer(1))
An inequality (0, -1, 0) x + 1 >= 0
>>> cube.Hrepresentation(Integer(2))
An inequality (0, 0, -1) x + 1 >= 0
>>> rectangular_box_points([Integer(0)]*Integer(3), [Integer(1)]*Integer(3), cube, return_saturated=True)
(((0, 0, 0), frozenset()),
 ((0, 0, 1), frozenset({2})),
 ((0, 1, 0), frozenset({1})),
 ((0, 1, 1), frozenset({1, 2})),
 ((1, 0, 0), frozenset({0})),
 ((1, 0, 1), frozenset({0, 2})),
 ((1, 1, 0), frozenset({0, 1})),
 ((1, 1, 1), frozenset({0, 1, 2})))

Sage Live

cube = polytopes.cube()
cube.Hrepresentation(0)
cube.Hrepresentation(1)
cube.Hrepresentation(2)
rectangular_box_points([0]*3, [1]*3, cube, return_saturated=True)

sage.geometry.integral_points.simplex_points(vertices)[source]¶

Return the integral points in a lattice simplex.

INPUT:

vertices – an iterable of integer coordinate vectors. The indices of vertices that span the simplex under consideration.

OUTPUT:

A tuple containing the integral point coordinates as $Z$ -vectors.

EXAMPLES:

Sage

sage: from sage.geometry.integral_points import simplex_points
sage: simplex_points([(1,2,3), (2,3,7), (-2,-3,-11)])
((-2, -3, -11), (0, 0, -2), (1, 2, 3), (2, 3, 7))

Python

>>> from sage.all import *
>>> from sage.geometry.integral_points import simplex_points
>>> simplex_points([(Integer(1),Integer(2),Integer(3)), (Integer(2),Integer(3),Integer(7)), (-Integer(2),-Integer(3),-Integer(11))])
((-2, -3, -11), (0, 0, -2), (1, 2, 3), (2, 3, 7))

Sage Live

from sage.geometry.integral_points import simplex_points
simplex_points([(1,2,3), (2,3,7), (-2,-3,-11)])

The simplex need not be full-dimensional:

Sage

sage: simplex = Polyhedron([(1,2,3,5), (2,3,7,5), (-2,-3,-11,5)])
sage: simplex_points(simplex.Vrepresentation())
((2, 3, 7, 5), (0, 0, -2, 5), (-2, -3, -11, 5), (1, 2, 3, 5))

sage: simplex_points([(2,3,7)])
((2, 3, 7),)

Python

>>> from sage.all import *
>>> simplex = Polyhedron([(Integer(1),Integer(2),Integer(3),Integer(5)), (Integer(2),Integer(3),Integer(7),Integer(5)), (-Integer(2),-Integer(3),-Integer(11),Integer(5))])
>>> simplex_points(simplex.Vrepresentation())
((2, 3, 7, 5), (0, 0, -2, 5), (-2, -3, -11, 5), (1, 2, 3, 5))

>>> simplex_points([(Integer(2),Integer(3),Integer(7))])
((2, 3, 7),)

Sage Live

simplex = Polyhedron([(1,2,3,5), (2,3,7,5), (-2,-3,-11,5)])
simplex_points(simplex.Vrepresentation())
simplex_points([(2,3,7)])