Ticket #11429: trac_11429_cythonize_lattice_points.patch

File trac_11429_cythonize_lattice_points.patch, 52.0 KB (added by vbraun, 11 years ago)

Updated patch

  • module_list.py

    # HG changeset patch
    # User Volker Braun <vbraun@stp.dias.ie>
    # Date 1323988670 0
    # Node ID 34cc23b585e0eb2958a554f928cd32aed04df6be
    # Parent  29bc5712b1befad488c789e14783f4f9d834ac53
    Trac #11429: Count integral points without PALP
    
    This patch implements the basic algorithms in Cython
    and optimizes the heck out of them. Its fast but sometimes
    not easy to read. No functionality visible to the
    end-user is added, the purpose of this patch is purely
    speed.
    
    diff --git a/module_list.py b/module_list.py
    a b  
    297297               sources = ['sage/geometry/toric_lattice_element.pyx'],
    298298               libraries=['gmp']),
    299299
     300     Extension('sage.geometry.integral_points',
     301               sources = ['sage/geometry/integral_points.pyx']),
     302
    300303     Extension('sage.geometry.triangulation.base',
    301304               sources = ['sage/geometry/triangulation/functions.cc',
    302305                          'sage/geometry/triangulation/data.cc',
  • sage/geometry/cone.py

    diff --git a/sage/geometry/cone.py b/sage/geometry/cone.py
    a b  
    189189from sage.libs.ppl import C_Polyhedron, Generator_System, Constraint_System, \
    190190    Linear_Expression, ray as PPL_ray, point as PPL_point, \
    191191    Poly_Con_Relation
    192 from sage.combinat.cartesian_product import CartesianProduct
     192from sage.geometry.integral_points import parallelotope_points
    193193
    194194
    195195def is_Cone(x):
     
    34993499                gens.update(cone.semigroup_generators())
    35003500            return tuple(gens)
    35013501
    3502         if self.is_trivial():
    3503             return tuple()
    3504         gens = list(parallelotope_points(self.rays())) + list(self.rays())
     3502        gens = list(parallelotope_points(self.rays(), N)) + list(self.rays())
    35053503        gens = filter(lambda v:gcd(v)==1, gens)
    35063504        return tuple(gens)
    35073505
     
    36813679        return vector(ZZ, p.get_values(x))
    36823680
    36833681
    3684 
    3685 #################################################################
    3686 def parallelotope_points(spanning_points):
    3687     r"""
    3688     Return integral points in the parallelotope spanned by a the
    3689     rays.
    3690    
    3691     See :meth:`~ConvexRationalPolyhedralCone.semigroup_generators` for a description of the
    3692     algorithm.
    3693 
    3694     INPUT:
    3695 
    3696     - ``spanning_points`` -- a non-empty list of linearly independent
    3697       rays (`\ZZ`-vectors or :class:`toric lattice
    3698       <sage.geometry.toric_lattice.ToricLatticeFactory>` elements),
    3699       not necessarily primitive lattice points.
    3700 
    3701     OUTPUT:
    3702 
    3703     The tuple of all lattice points in the half-open parallelotope
    3704     spanned by the rays `r_i`,
    3705    
    3706     .. math::
    3707 
    3708         \mathop{par}(\{r_i\}) =
    3709         \sum_{0\leq a_i < 1} a_i r_i
    3710    
    3711     By half-open parallelotope, we mean that the
    3712     points in the facets not meeting the origin are omitted.
    3713    
    3714     EXAMPLES:
    3715    
    3716     Note how the points on the outward-facing factes are omitted::
    3717 
    3718         sage: from sage.geometry.cone import parallelotope_points
    3719         sage: rays = map(vector, [(2,0), (0,2)])
    3720         sage: parallelotope_points(rays)
    3721         ((0, 0), (1, 0), (0, 1), (1, 1))
    3722 
    3723     The rays can also be toric lattice points::
    3724 
    3725         sage: rays = map(ToricLattice(2), [(2,0), (0,2)])
    3726         sage: parallelotope_points(rays)
    3727         (N(0, 0), N(1, 0), N(0, 1), N(1, 1))
    3728 
    3729     A non-smooth cone::
    3730 
    3731         sage: c = Cone([ (1,0), (1,2) ])
    3732         sage: parallelotope_points(c.rays())
    3733         (N(0, 0), N(1, 1))
    3734 
    3735     A ``ValueError`` is raised if the ``spanning_points`` are not
    3736     linearly independent::
    3737 
    3738         sage: rays = map(ToricLattice(2), [(1,1)]*2)
    3739         sage: parallelotope_points(rays)
    3740         Traceback (most recent call last):
    3741         ...
    3742         ValueError: The spanning points are not linearly independent!
    3743     """
    3744     N = spanning_points[0].parent()  # this is why there needs to be at least one ray
    3745 
    3746     # compute points in the open parallelotope, see [BrunsKoch]
    3747     R = matrix(spanning_points).transpose()
    3748     D,U,V = R.smith_form()
    3749     e = D.diagonal()          # the elementary divisors
    3750     d = prod(e)               # the determinant
    3751     if d==0:
    3752         raise ValueError('The spanning points are not linearly independent!')
    3753     u = U.inverse().columns() # generators for gp(semigroup)
    3754     gens = []
    3755    
    3756     # "inverse" of the ray matrix as far as possible over ZZ
    3757     # R*Rinv == diagonal_matrix([d]*D.ncols() + [0]*(D.nrows()-D.ncols())) 
    3758     # If R is full rank, this is Rinv = matrix(ZZ, R.inverse() * d)
    3759     Dinv = D.transpose()
    3760     for i in range(0,D.ncols()):
    3761         Dinv[i,i] = d/D[i,i]   
    3762     Rinv = V * Dinv * U
    3763 
    3764     for b in CartesianProduct(*[ range(0,i) for i in e ]):
    3765         # this is our generator modulo the lattice spanned by the rays
    3766         gen_mod_rays = sum( b_i*u_i for b_i, u_i in zip(b,u) )
    3767         q_times_d = Rinv * gen_mod_rays
    3768         q_times_d = vector(ZZ,[ q_i % d  for q_i in q_times_d ])
    3769         gen = N(R*q_times_d / d)
    3770         gen.set_immutable()
    3771         gens.append(gen)
    3772     assert(len(gens) == d)
    3773     return tuple(gens)
    3774 
  • new file sage/geometry/integral_points.pyx

    diff --git a/sage/geometry/integral_points.pyx b/sage/geometry/integral_points.pyx
    new file mode 100644
    - +  
     1# cython: profile=True
     2
     3r"""
     4Cython helper methods to compute integral points in polyhedra.
     5"""
     6
     7#*****************************************************************************
     8#       Copyright (C) 2010 Volker Braun <vbraun.name@gmail.com>
     9#
     10#  Distributed under the terms of the GNU General Public License (GPL)
     11#
     12#                  http://www.gnu.org/licenses/
     13#*****************************************************************************
     14
     15from sage.matrix.constructor import matrix, column_matrix, vector, diagonal_matrix
     16from sage.rings.all import QQ, RR, ZZ, gcd, lcm
     17from sage.combinat.permutation import Permutation
     18from sage.combinat.cartesian_product import CartesianProduct
     19from sage.misc.all import prod
     20import copy
     21
     22##############################################################################
     23# The basic idea to enumerate the lattice points in the parallelotope
     24# is to use the Smith normal form of the ray coordinate matrix to
     25# bring the lattice into a "nice" basis. Here is the straightforward
     26# implementation. Note that you do not need to reduce to the
     27# full-dimensional case, the Smith normal form takes care of that for
     28# you.
     29#
     30## def parallelotope_points(spanning_points, lattice):
     31##     # compute points in the open parallelotope, see [BrunsKoch]
     32##     R = matrix(spanning_points).transpose()
     33##     D,U,V = R.smith_form()
     34##     e = D.diagonal()          # the elementary divisors
     35##     d = prod(e)               # the determinant
     36##     u = U.inverse().columns() # generators for gp(semigroup)
     37##     
     38##     # "inverse" of the ray matrix as far as possible over ZZ
     39##     # R*Rinv == diagonal_matrix([d]*D.ncols() + [0]*(D.nrows()-D.ncols())) 
     40##     # If R is full rank, this is Rinv = matrix(ZZ, R.inverse() * d)
     41##     Dinv = D.transpose()
     42##     for i in range(0,D.ncols()):
     43##         Dinv[i,i] = d/D[i,i]   
     44##     Rinv = V * Dinv * U
     45##
     46##     gens = []
     47##     for b in CartesianProduct(*[ range(0,i) for i in e ]):
     48##         # this is our generator modulo the lattice spanned by the rays
     49##         gen_mod_rays = sum( b_i*u_i for b_i, u_i in zip(b,u) )
     50##         q_times_d = Rinv * gen_mod_rays
     51##         q_times_d = vector(ZZ,[ q_i % d  for q_i in q_times_d ])
     52##         gen = lattice(R*q_times_d / d)
     53##         gen.set_immutable()
     54##         gens.append(gen)
     55##     assert(len(gens) == d)
     56##     return tuple(gens)
     57#
     58# The problem with the naive implementation is that it is slow:
     59#
     60#   1. You can simplify some of the matrix multiplications
     61#
     62#   2. The inner loop keeps creating new matrices and vectors, which
     63#      is slow. It needs to recycle objects. Instead of creating a new
     64#      lattice point again and again, change the entries of an
     65#      existing lattice point and then copy it!
     66
     67
     68def parallelotope_points(spanning_points, lattice):
     69    r"""
     70    Return integral points in the parallelotope starting at the origin
     71    and spanned by the ``spanning_points``.
     72   
     73    See :meth:`~ConvexRationalPolyhedralCone.semigroup_generators` for a description of the
     74    algorithm.
     75
     76    INPUT:
     77
     78    - ``spanning_points`` -- a non-empty list of linearly independent
     79      rays (`\ZZ`-vectors or :class:`toric lattice
     80      <sage.geometry.toric_lattice.ToricLatticeFactory>` elements),
     81      not necessarily primitive lattice points.
     82
     83    OUTPUT:
     84
     85    The tuple of all lattice points in the half-open parallelotope
     86    spanned by the rays `r_i`,
     87   
     88    .. math::
     89
     90        \mathop{par}(\{r_i\}) =
     91        \sum_{0\leq a_i < 1} a_i r_i
     92   
     93    By half-open parallelotope, we mean that the
     94    points in the facets not meeting the origin are omitted.
     95   
     96    EXAMPLES:
     97   
     98    Note how the points on the outward-facing factes are omitted::
     99
     100        sage: from sage.geometry.integral_points import parallelotope_points
     101        sage: rays = map(vector, [(2,0), (0,2)])
     102        sage: parallelotope_points(rays, ZZ^2)
     103        ((0, 0), (1, 0), (0, 1), (1, 1))
     104
     105    The rays can also be toric lattice points::
     106
     107        sage: rays = map(ToricLattice(2), [(2,0), (0,2)])
     108        sage: parallelotope_points(rays, ToricLattice(2))
     109        (N(0, 0), N(1, 0), N(0, 1), N(1, 1))
     110
     111    A non-smooth cone::
     112
     113        sage: c = Cone([ (1,0), (1,2) ])
     114        sage: parallelotope_points(c.rays(), c.lattice())
     115        (N(0, 0), N(1, 1))
     116
     117    A ``ValueError`` is raised if the ``spanning_points`` are not
     118    linearly independent::
     119
     120        sage: rays = map(ToricLattice(2), [(1,1)]*2)
     121        sage: parallelotope_points(rays, ToricLattice(2))
     122        Traceback (most recent call last):
     123        ...
     124        ValueError: The spanning points are not linearly independent!
     125
     126    TESTS::
     127
     128        sage: rays = map(vector,[(-3, -2, -3, -2), (-2, -1, -8, 5), (1, 9, -7, -4), (-3, -1, -2, 2)])
     129        sage: len(parallelotope_points(rays, ZZ^4))
     130        967
     131    """
     132    R = matrix(spanning_points).transpose()
     133    e, d, VDinv = ray_matrix_normal_form(R)
     134    points = loop_over_parallelotope_points(e, d, VDinv, R, lattice)
     135    for p in points:
     136        p.set_immutable()
     137    return points
     138
     139
     140def ray_matrix_normal_form(R):
     141    r"""
     142    Compute the Smith normal form of the ray matrix for
     143    :func:`parallelotope_points`.
     144
     145    INPUT:
     146
     147    - ``R`` -- `\ZZ`-matrix whose columns are the rays spanning the
     148      parallelotope.
     149
     150    OUTPUT:
     151
     152    A tuple containing ``e``, ``d``, and ``VDinv``.
     153
     154    EXAMPLES::
     155   
     156        sage: from sage.geometry.integral_points import ray_matrix_normal_form
     157        sage: R = column_matrix(ZZ,[3,3,3])
     158        sage: ray_matrix_normal_form(R)
     159        ([3], 3, [1])
     160    """
     161    D,U,V = R.smith_form()
     162    e = D.diagonal()            # the elementary divisors
     163    d = prod(e)                 # the determinant
     164    if d==0:
     165        raise ValueError('The spanning points are not linearly independent!')
     166    Dinv = diagonal_matrix(ZZ,[ d/D[i,i] for i in range(0,D.ncols()) ])
     167    VDinv = V * Dinv
     168    return (e,d,VDinv)
     169
     170
     171
     172# The optimized version avoids constructing new matrices, vectors, and lattice points
     173cpdef loop_over_parallelotope_points(e, d, VDinv, R, lattice, A=None, b=None):
     174    r"""
     175    The inner loop of :func:`parallelotope_points`.
     176   
     177    INPUT:
     178
     179    See :meth:`parallelotope_points` for ``e``, ``d``, ``VDinv``, ``R``, ``lattice``.
     180
     181    - ``A``, ``b``: Either both ``None`` or a vector and number. If
     182      present, only the parallelotope points satisfying `A x \leq b`
     183      are returned.
     184   
     185    OUTPUT:
     186
     187    The points of the half-open parallelotope as a tuple of lattice
     188    points.
     189
     190    EXAMPLES:
     191
     192        sage: e = [3]
     193        sage: d = prod(e)
     194        sage: VDinv = matrix(ZZ, [[1]])
     195        sage: R = column_matrix(ZZ,[3,3,3])
     196        sage: lattice = ZZ^3
     197        sage: from sage.geometry.integral_points import loop_over_parallelotope_points
     198        sage: loop_over_parallelotope_points(e, d, VDinv, R, lattice)
     199        ((0, 0, 0), (1, 1, 1), (2, 2, 2))
     200
     201        sage: A = vector(ZZ, [1,0,0])
     202        sage: b = 1
     203        sage: loop_over_parallelotope_points(e, d, VDinv, R, lattice, A, b)
     204        ((0, 0, 0), (1, 1, 1))
     205    """
     206    cdef int i, j
     207    cdef int dim = VDinv.nrows()
     208    cdef int ambient_dim = R.nrows()
     209    gens = []
     210    s = ZZ.zero() # summation variable
     211    gen = lattice(0)
     212    q_times_d = vector(ZZ, dim)
     213    for base in CartesianProduct(*[ range(0,i) for i in e ]):
     214        for i in range(0, dim):
     215            s = ZZ.zero()
     216            for j in range(0, dim):
     217                s += VDinv[i,j] * base[j]
     218            q_times_d[i] = s % d
     219        for i in range(0, ambient_dim):
     220            s = ZZ.zero()
     221            for j in range(0, dim):
     222                s += R[i,j] * q_times_d[j]
     223            gen[i] = s / d
     224        if A is not None:
     225            s = ZZ.zero()
     226            for i in range(0, ambient_dim):
     227                s += A[i] * gen[i]
     228            if s>b:
     229                continue
     230        gens.append(copy.copy(gen))
     231    if A is None:
     232        assert(len(gens) == d)
     233    return tuple(gens)
     234
     235
     236
     237##############################################################################
     238def simplex_points(vertices):
     239    r"""
     240    Return the integral points in a lattice simplex.
     241   
     242    INPUT:
     243
     244    - ``vertices`` -- an iterable of integer coordinate vectors. The
     245      indices of vertices that span the simplex under
     246      consideration.
     247
     248    OUTPUT:
     249
     250    A tuple containing the integral point coordinates as `\ZZ`-vectors.
     251
     252    EXAMPLES::
     253
     254        sage: from sage.geometry.integral_points import simplex_points
     255        sage: simplex_points([(1,2,3), (2,3,7), (-2,-3,-11)])
     256        ((-2, -3, -11), (0, 0, -2), (1, 2, 3), (2, 3, 7))
     257
     258    The simplex need not be full-dimensional::
     259
     260        sage: simplex = Polyhedron([(1,2,3,5), (2,3,7,5), (-2,-3,-11,5)])
     261        sage: simplex_points(simplex.Vrepresentation())
     262        ((-2, -3, -11, 5), (0, 0, -2, 5), (1, 2, 3, 5), (2, 3, 7, 5))
     263
     264        sage: simplex_points([(2,3,7)])
     265        ((2, 3, 7),)
     266
     267    TESTS::
     268   
     269        sage: v = [(1,0,7,-1), (-2,-2,4,-3), (-1,-1,-1,4), (2,9,0,-5), (-2,-1,5,1)]
     270        sage: simplex = Polyhedron(v); simplex
     271        A 4-dimensional polyhedron in QQ^4 defined as the convex hull of 5 vertices.
     272        sage: pts = simplex_points(simplex.Vrepresentation())
     273        sage: len(pts)
     274        49
     275        sage: for p in pts: p.set_immutable()
     276        sage: len(set(pts))
     277        49
     278
     279        sage: all(simplex.contains(p) for p in pts)
     280        True
     281       
     282        sage: v = [(4,-1,-1,-1), (-1,4,-1,-1), (-1,-1,4,-1), (-1,-1,-1,4), (-1,-1,-1,-1)]
     283        sage: P4mirror = Polyhedron(v); P4mirror
     284        A 4-dimensional polyhedron in QQ^4 defined as the convex hull of 5 vertices.
     285        sage: len( simplex_points(P4mirror.Vrepresentation()) )
     286        126
     287    """
     288    rays = [ vector(ZZ, v) for v in vertices ]
     289    if len(rays)==0:
     290        return tuple()
     291    origin = rays.pop()
     292    origin.set_immutable()
     293    if len(rays)==0:
     294        return tuple([origin])
     295    translate_points(rays, origin)
     296
     297    # Find equation Ax<=b that cuts out simplex from parallelotope
     298    Rt = matrix(rays)
     299    R = Rt.transpose()
     300    if R.is_square():
     301        b = abs(R.det())
     302        A = R.solve_left(vector([b]*len(rays)))
     303    else:
     304        RtR = Rt * R
     305        b = abs(RtR.det())
     306        A = RtR.solve_left(vector([b]*len(rays))) * Rt
     307   
     308    # e, d, VDinv = ray_matrix_normal_form(R)
     309    #    print origin
     310    #    print rays
     311    #    print parallelotope_points(rays, origin.parent())
     312    #    print 'A = ', A
     313    #    print 'b = ', b
     314
     315    e, d, VDinv = ray_matrix_normal_form(R)
     316    lattice = origin.parent()
     317    points = loop_over_parallelotope_points(e, d, VDinv, R, lattice, A, b) + tuple(rays)
     318    translate_points(points, -origin)
     319    return points
     320
     321
     322cdef translate_points(v_list, delta):
     323    r"""
     324    Add ``delta`` to each vector in ``v_list``.
     325    """
     326    cdef int dim = delta.degree()
     327    for v in v_list:
     328        for i in range(0,dim):
     329            v[i] -= delta[i]
     330
     331
     332
     333##############################################################################
     334# For points with "small" coordinates (that is, fitting into a small
     335# rectangular bounding box) it is faster to naively enumerate the
     336# points. This saves the overhead of triangulating the polytope etc.
     337
     338def rectangular_box_points(box_min, box_max, polyhedron=None):
     339    r"""
     340    Return the integral points in the lattice bounding box that are
     341    also contained in the given polyhedron.
     342
     343    INPUT:
     344
     345    - ``box_min`` -- A list of integers. The minimal value for each
     346      coordinate of the rectangular bounding box.
     347
     348    - ``box_max`` -- A list of integers. The maximal value for each
     349      coordinate of the rectangular bounding box.
     350
     351    - ``polyhedron`` -- A :class:`~sage.geometry.polyhedra.Polyhedron`
     352      or ``None`` (default).
     353
     354    OUTPUT:
     355
     356    A tuple containing the integral points of the rectangular box
     357    spanned by `box_min` and `box_max` and that lie inside the
     358    ``polyhedron``. For sufficiently large bounding boxes, this are
     359    all integral points of the polyhedron.
     360
     361    If no polyhedron is specified, all integral points of the
     362    rectangular box are returned.
     363
     364    ALGORITHM:
     365
     366    This function implements the naive algorithm towards counting
     367    integral points. Given min and max of vertex coordinates, it
     368    iterates over all points in the bounding box and checks whether
     369    they lie in the polyhedron. The following optimizations are
     370    implemented:
     371
     372      * Cython: Use machine integers and optimizing C/C++ compiler
     373        where possible, arbitrary precision integers where necessary.
     374        Bounds checking, no compile time limits.
     375       
     376      * Unwind inner loop (and next-to-inner loop):
     377     
     378        .. math::
     379
     380            Ax\leq b
     381            \quad \Leftrightarrow \quad
     382            a_1 x_1 ~\leq~ b - \sum_{i=2}^d a_i x_i
     383
     384        so we only have to evaluate `a_1 * x_1` in the inner loop.
     385
     386      * Coordinates are permuted to make the longest box edge the
     387        inner loop. The inner loop is optimized to run very fast, so
     388        its best to do as much work as possible there.
     389 
     390      * Continuously reorder inequalities and test the most
     391        restrictive inequalities first.
     392
     393      * Use convexity and only find first and last allowed point in
     394        the inner loop. The points in-between must be points of the
     395        polyhedron, too.
     396
     397    EXAMPLES::
     398
     399        sage: from sage.geometry.integral_points import rectangular_box_points
     400        sage: rectangular_box_points([0,0,0],[1,2,3])
     401        ((0, 0, 0), (0, 0, 1), (0, 0, 2), (0, 0, 3),
     402         (0, 1, 0), (0, 1, 1), (0, 1, 2), (0, 1, 3),
     403         (0, 2, 0), (0, 2, 1), (0, 2, 2), (0, 2, 3),
     404         (1, 0, 0), (1, 0, 1), (1, 0, 2), (1, 0, 3),
     405         (1, 1, 0), (1, 1, 1), (1, 1, 2), (1, 1, 3),
     406         (1, 2, 0), (1, 2, 1), (1, 2, 2), (1, 2, 3))
     407
     408        sage: cell24 = polytopes.twenty_four_cell()
     409        sage: rectangular_box_points([-1]*4, [1]*4, cell24)
     410        ((-1, 0, 0, 0), (0, -1, 0, 0), (0, 0, -1, 0), (0, 0, 0, -1),
     411         (0, 0, 0, 0),
     412         (0, 0, 0, 1), (0, 0, 1, 0), (0, 1, 0, 0), (1, 0, 0, 0))
     413        sage: d = 3
     414        sage: dilated_cell24 = d*cell24
     415        sage: len( rectangular_box_points([-d]*4, [d]*4, dilated_cell24) )
     416        305
     417
     418        sage: d = 6
     419        sage: dilated_cell24 = d*cell24
     420        sage: len( rectangular_box_points([-d]*4, [d]*4, dilated_cell24) )
     421        3625
     422
     423        sage: polytope = Polyhedron([(-4,-3,-2,-1),(3,1,1,1),(1,2,1,1),(1,1,3,0),(1,3,2,4)])
     424        sage: pts = rectangular_box_points([-4]*4, [4]*4, polytope); pts
     425        ((-4, -3, -2, -1), (-1, 0, 0, 1), (0, 1, 1, 1), (1, 1, 1, 1), (1, 1, 3, 0),
     426         (1, 2, 1, 1), (1, 2, 2, 2), (1, 3, 2, 4), (2, 1, 1, 1), (3, 1, 1, 1))
     427        sage: all(polytope.contains(p) for p in pts)
     428        True
     429       
     430        sage: set(map(tuple,pts)) == \
     431        ...   set([(-4,-3,-2,-1),(3,1,1,1),(1,2,1,1),(1,1,3,0),(1,3,2,4),
     432        ...        (0,1,1,1),(1,2,2,2),(-1,0,0,1),(1,1,1,1),(2,1,1,1)])   # computed with PALP
     433        True
     434
     435    Long ints and non-integral polyhedra are explictly allowed::
     436   
     437        sage: polytope = Polyhedron([[1], [10*pi.n()]], field=RDF)
     438        sage: len( rectangular_box_points([-100], [100], polytope) )
     439        31
     440
     441        sage: halfplane = Polyhedron(ieqs=[(-1,1,0)])
     442        sage: rectangular_box_points([0,-1+10^50], [0,1+10^50])
     443        ((0, 99999999999999999999999999999999999999999999999999),
     444         (0, 100000000000000000000000000000000000000000000000000),
     445         (0, 100000000000000000000000000000000000000000000000001))
     446        sage: len( rectangular_box_points([0,-100+10^50], [1,100+10^50], halfplane) )
     447        201
     448    """
     449    assert len(box_min)==len(box_max)
     450    cdef int d = len(box_min)
     451    diameter = sorted([ (box_max[i]-box_min[i], i) for i in range(0,d) ], reverse=True)
     452    diameter_value = [ x[0] for x in diameter ]
     453    diameter_index = [ x[1] for x in diameter ]
     454
     455    sort_perm = Permutation([ i+1 for i in diameter_index])
     456    orig_perm = sort_perm.inverse()
     457    orig_perm_list = [ i-1 for i in orig_perm ]
     458    box_min = sort_perm.action(box_min)
     459    box_max = sort_perm.action(box_max)
     460    inequalities = InequalityCollection(polyhedron, sort_perm, box_min, box_max)
     461
     462    v = vector(ZZ, d)
     463    points = []
     464    for p in loop_over_rectangular_box_points(box_min, box_max, inequalities, d):
     465        #  v = vector(ZZ, orig_perm.action(p))   # too slow
     466        for i in range(0,d):
     467            v[i] = p[orig_perm_list[i]]
     468        points.append(copy.copy(v))
     469    return tuple(points)
     470
     471
     472cdef loop_over_rectangular_box_points(box_min, box_max, inequalities, int d):
     473    """
     474    The inner loop of :func:`rectangular_box_points`.
     475   
     476    INPUT:
     477   
     478    - ``box_min``, ``box_max`` -- the bounding box.
     479
     480    - ``inequalities`` -- a :class:`InequalityCollection` containing
     481      the inequalities defining the polyhedron.
     482
     483    - ``d`` -- the ambient space dimension.
     484
     485    OUTPUT:
     486
     487    The integral points in the bounding box satisfying all
     488    inequalities.
     489    """
     490    cdef int inc
     491    points = []
     492    p = copy.copy(box_min)
     493    inequalities.prepare_next_to_inner_loop(p)
     494    while True:
     495        inequalities.prepare_inner_loop(p)
     496        i_min = box_min[0]
     497        i_max = box_max[0]
     498        # Find the lower bound for the allowed region
     499        while i_min <= i_max:
     500            if inequalities.are_satisfied(i_min):
     501                break
     502            i_min += 1
     503        # Find the upper bound for the allowed region
     504        while i_min <= i_max:
     505            if inequalities.are_satisfied(i_max):
     506                break
     507            i_max -= 1
     508        # The points i_min .. i_max are contained in the polyhedron
     509        i = i_min
     510        while i <= i_max:
     511            p[0] = i
     512            points.append(tuple(p))
     513            i += 1
     514        # finally increment the other entries in p to move on to next inner loop
     515        inc = 1
     516        if d==1: return points
     517        while True:
     518            if p[inc]==box_max[inc]:
     519                p[inc] = box_min[inc]
     520                inc += 1
     521                if inc==d:
     522                    return points
     523            else:
     524                p[inc] += 1
     525                break
     526        if inc>1:
     527            inequalities.prepare_next_to_inner_loop(p)
     528
     529
     530
     531cdef class Inequality_generic:
     532    """
     533    An inequality whose coefficients are arbitrary Python/Sage objects
     534
     535    INPUT:
     536
     537    - ``A`` -- list of integers.
     538
     539    - ``b`` -- integer
     540
     541    OUTPUT:
     542
     543    Inequality `A x + b \geq 0`.
     544
     545    EXAMPLES::
     546
     547        sage: from sage.geometry.integral_points import Inequality_generic
     548        sage: Inequality_generic([2*pi,sqrt(3),7/2], -5.5)
     549        generic: (2*pi, sqrt(3), 7/2) x + -5.50000000000000 >= 0
     550    """
     551
     552    cdef object A
     553    cdef object b
     554    cdef object coeff
     555    cdef object cache
     556
     557    def __cinit__(self, A, b):
     558        """
     559        The Cython constructor
     560
     561        INPUT:
     562
     563        See :class:`Inequality_generic`.
     564
     565        EXAMPLES::
     566
     567            sage: from sage.geometry.integral_points import Inequality_generic
     568            sage: Inequality_generic([2*pi,sqrt(3),7/2], -5.5)
     569            generic: (2*pi, sqrt(3), 7/2) x + -5.50000000000000 >= 0
     570        """
     571        self.A = A
     572        self.b = b
     573        self.coeff = 0
     574        self.cache = 0
     575
     576    def __repr__(self):
     577        """
     578        Return a string representation.
     579
     580        OUTPUT:
     581
     582        String.
     583
     584        EXAMPLES::
     585
     586            sage: from sage.geometry.integral_points import Inequality_generic
     587            sage: Inequality_generic([2,3,7], -5).__repr__()
     588            'generic: (2, 3, 7) x + -5 >= 0'
     589        """
     590        s = 'generic: ('
     591        s += ', '.join([str(self.A[i]) for i in range(0,len(self.A))])
     592        s += ') x + ' + str(self.b) + ' >= 0'
     593        return s
     594
     595    cdef prepare_next_to_inner_loop(self, p):
     596        """
     597        In :class:`Inequality_int` this method is used to peel of the
     598        next-to-inner loop.
     599       
     600        See :meth:`InequalityCollection.prepare_inner_loop` for more details.
     601        """
     602        pass
     603
     604    cdef prepare_inner_loop(self, p):
     605        """
     606        Peel off the inner loop.
     607       
     608        See :meth:`InequalityCollection.prepare_inner_loop` for more details.
     609        """
     610        cdef int j
     611        self.coeff = self.A[0]
     612        self.cache = self.b
     613        for j in range(1,len(self.A)):
     614            self.cache += self.A[j] * p[j]
     615
     616    cdef bint is_not_satisfied(self, inner_loop_variable):
     617        r"""
     618        Test the inequality, using the cached value from :meth:`prepare_inner_loop`
     619        """
     620        return inner_loop_variable*self.coeff + self.cache < 0       
     621
     622
     623
     624DEF INEQ_INT_MAX_DIM = 20
     625
     626cdef class Inequality_int:
     627    """
     628    Fast version of inequality in the case that all coefficient fit
     629    into machine ints.
     630
     631    INPUT:
     632
     633    - ``A`` -- list of integers.
     634
     635    - ``b`` -- integer
     636
     637    - ``max_abs_coordinates`` -- the maximum of the coordinates that
     638      one wants to evalate the coordinates on. Used for overflow
     639      checking.
     640
     641    OUTPUT:
     642
     643    Inequality `A x + b \geq 0`. A ``OverflowError`` is raised if a
     644    machine integer is not long enough to hold the results. A
     645    ``ValueError`` is raised if some of the input is not integral.
     646
     647    EXAMPLES::
     648
     649        sage: from sage.geometry.integral_points import Inequality_int
     650        sage: Inequality_int([2,3,7], -5, [10]*3)
     651        integer: (2, 3, 7) x + -5 >= 0
     652
     653        sage: Inequality_int([1]*21, -5, [10]*21)
     654        Traceback (most recent call last):
     655        ...
     656        OverflowError: Dimension limit exceeded.
     657
     658        sage: Inequality_int([2,3/2,7], -5, [10]*3)
     659        Traceback (most recent call last):
     660        ...
     661        ValueError: Not integral.
     662
     663        sage: Inequality_int([2,3,7], -5.2, [10]*3)
     664        Traceback (most recent call last):
     665        ...
     666        ValueError: Not integral.
     667
     668        sage: Inequality_int([2,3,7], -5*10^50, [10]*3)  # actual error message can differ between 32 and 64 bit
     669        Traceback (most recent call last):
     670        ...
     671        OverflowError: ...
     672    """
     673    cdef int A[INEQ_INT_MAX_DIM]
     674    cdef int b
     675    cdef int dim
     676    # the innermost coefficient
     677    cdef int coeff
     678    cdef int cache
     679    # the next-to-innermost coefficient
     680    cdef int coeff_next
     681    cdef int cache_next
     682
     683    cdef int _to_int(self, x) except? -999:
     684        if not x in ZZ: raise ValueError('Not integral.')
     685        return int(x)  # raises OverflowError in Cython if necessary
     686
     687    def __cinit__(self, A, b, max_abs_coordinates):
     688        """
     689        The Cython constructor
     690
     691        See :class:`Inequality_int` for input.
     692
     693        EXAMPLES::
     694
     695            sage: from sage.geometry.integral_points import Inequality_int
     696            sage: Inequality_int([2,3,7], -5, [10]*3)
     697            integer: (2, 3, 7) x + -5 >= 0
     698        """
     699        cdef int i
     700        self.dim = self._to_int(len(A))
     701        if self.dim<1 or self.dim>INEQ_INT_MAX_DIM:
     702            raise OverflowError('Dimension limit exceeded.')
     703        for i in range(0,self.dim):
     704            self.A[i] = self._to_int(A[i])
     705        self.b = self._to_int(b)
     706        self.coeff = self.A[0]
     707        if self.dim>0:
     708            self.coeff_next = self.A[1]
     709        # finally, make sure that there cannot be any overflow during the enumeration
     710        self._to_int(sum( [ZZ(b)]+[ZZ(A[i])*ZZ(max_abs_coordinates[i]) for i in range(0,self.dim)] ))
     711
     712    def __repr__(self):
     713        """
     714        Return a string representation.
     715
     716        OUTPUT:
     717
     718        String.
     719
     720        EXAMPLES::
     721
     722            sage: from sage.geometry.integral_points import Inequality_int
     723            sage: Inequality_int([2,3,7], -5, [10]*3).__repr__()
     724            'integer: (2, 3, 7) x + -5 >= 0'
     725        """
     726        s = 'integer: ('
     727        s += ', '.join([str(self.A[i]) for i in range(0,self.dim)])
     728        s += ') x + ' + str(self.b) + ' >= 0'
     729        return s
     730
     731    cdef prepare_next_to_inner_loop(Inequality_int self, p):
     732        """
     733        Peel off the next-to-inner loop.
     734       
     735        See :meth:`InequalityCollection.prepare_inner_loop` for more details.
     736        """
     737        cdef int j
     738        self.cache_next = self.b
     739        for j in range(2,self.dim):
     740            self.cache_next += self.A[j] * p[j]
     741
     742    cdef prepare_inner_loop(Inequality_int self, p):
     743        """
     744        Peel off the inner loop.
     745       
     746        See :meth:`InequalityCollection.prepare_inner_loop` for more details.
     747        """
     748        cdef int j
     749        if self.dim>1:
     750            self.cache = self.cache_next + self.coeff_next*p[1]
     751        else:
     752            self.cache = self.cache_next
     753
     754    cdef bint is_not_satisfied(Inequality_int self, int inner_loop_variable):
     755        return inner_loop_variable*self.coeff + self.cache < 0       
     756
     757
     758
     759cdef class InequalityCollection:
     760    """
     761    A collection of inequalities.
     762
     763    INPUT:
     764
     765    - ``polyhedron`` -- a polyhedron defining the inequalities.
     766
     767    - ``permutation`` -- a permution of the coordinates. Will be used
     768      to permute the coordinates of the inequality.
     769
     770    - ``box_min``, ``box_max`` -- the (not permuted) minimal and maximal
     771      coordinates of the bounding box. Used for bounds checking.
     772
     773    EXAMPLES::
     774
     775        sage: from sage.geometry.integral_points import InequalityCollection
     776        sage: P_QQ = Polyhedron(identity_matrix(3).columns() + [(-2, -1,-1)], field=QQ)
     777        sage: ieq = InequalityCollection(P_QQ, Permutation([1,2,3]), [0]*3,[1]*3); ieq
     778        The collection of inequalities
     779        integer: (-1, -1, -1) x + 1 >= 0
     780        integer: (-1, -1, 4) x + 1 >= 0
     781        integer: (-1, 4, -1) x + 1 >= 0
     782        integer: (3, -2, -2) x + 2 >= 0
     783
     784        sage: P_RR = Polyhedron(identity_matrix(2).columns() + [(-2.7, -1)], field=RDF)
     785        sage: InequalityCollection(P_RR, Permutation([1,2]), [0]*2, [1]*2)
     786        The collection of inequalities
     787        integer: (-1, -1) x + 1 >= 0
     788        generic: (-1.0, 3.7) x + 1.0 >= 0
     789        generic: (1.0, -1.35) x + 1.35 >= 0
     790
     791        sage: line = Polyhedron(eqns=[(2,3,7)])
     792        sage: InequalityCollection(line, Permutation([1,2]), [0]*2, [1]*2 )
     793        The collection of inequalities
     794        integer: (3, 7) x + 2 >= 0
     795        integer: (-3, -7) x + -2 >= 0
     796
     797    TESTS::
     798
     799        sage: TestSuite(ieq).run(skip='_test_pickling')
     800    """
     801    cdef object ineqs_int
     802    cdef object ineqs_generic
     803
     804    def __repr__(self):
     805        r"""
     806        Return a string representation.
     807       
     808        OUTPUT:
     809
     810        String.
     811
     812        EXAMPLES::
     813
     814            sage: from sage.geometry.integral_points import InequalityCollection
     815            sage: line = Polyhedron(eqns=[(2,3,7)])
     816            sage: InequalityCollection(line, Permutation([1,2]), [0]*2, [1]*2 ).__repr__()
     817            'The collection of inequalities\ninteger: (3, 7) x + 2 >= 0\ninteger: (-3, -7) x + -2 >= 0'
     818        """
     819        s = 'The collection of inequalities\n'
     820        for ineq in self.ineqs_int:
     821            s += str(<Inequality_int>ineq) + '\n'
     822        for ineq in self.ineqs_generic:
     823            s += str(<Inequality_generic>ineq) + '\n'
     824        return s.strip()
     825
     826    def _make_A_b(self, Hrep_obj, permutation):
     827        r"""
     828        Return the coefficients and constant of the H-representation
     829        object.
     830
     831        INPUT:
     832
     833        - ``Hrep_obj`` -- a H-representation object of the polyhedron.
     834
     835        - ``permutation`` -- the permutation of the coordinates to
     836          apply to ``A``.
     837
     838        OUTPUT:
     839
     840        A pair ``(A,b)``.
     841
     842        EXAXMPLES::
     843       
     844            sage: from sage.geometry.integral_points import InequalityCollection
     845            sage: line = Polyhedron(eqns=[(2,3,7)])
     846            sage: ieq = InequalityCollection(line, Permutation([1,2]), [0]*2, [1]*2 )
     847            sage: ieq._make_A_b(line.Hrepresentation(0), Permutation([1,2]))
     848            ([3, 7], 2)
     849            sage: ieq._make_A_b(line.Hrepresentation(0), Permutation([2,1]))
     850            ([7, 3], 2)
     851        """
     852        v = list(Hrep_obj)
     853        A = permutation.action(v[1:])
     854        b = v[0]
     855        try:
     856            x = lcm([a.denominator() for a in A] + [b.denominator()])
     857            A = [ a*x for a in A ]
     858            b = b*x
     859        except AttributeError:
     860            pass
     861        return (A,b)
     862       
     863    def __cinit__(self, polyhedron, permutation, box_min, box_max):
     864        """
     865        The Cython constructor
     866       
     867        See the class documentation for the desrciption of the arguments.
     868
     869        EXAMPLES::
     870
     871            sage: from sage.geometry.integral_points import InequalityCollection
     872            sage: line = Polyhedron(eqns=[(2,3,7)])
     873            sage: InequalityCollection(line, Permutation([1,2]), [0]*2, [1]*2 )
     874            The collection of inequalities
     875            integer: (3, 7) x + 2 >= 0
     876            integer: (-3, -7) x + -2 >= 0
     877        """
     878        max_abs_coordinates = [ max(abs(c_min), abs(c_max))
     879                                for c_min, c_max in zip(box_min, box_max) ]
     880        max_abs_coordinates = permutation.action(max_abs_coordinates)
     881        self.ineqs_int = []
     882        self.ineqs_generic = []
     883        if not polyhedron:
     884            return
     885        for Hrep_obj in polyhedron.inequality_generator():
     886            A, b = self._make_A_b(Hrep_obj, permutation)
     887            try:
     888                H = Inequality_int(A, b, max_abs_coordinates)
     889                self.ineqs_int.append(H)
     890            except (OverflowError, ValueError):
     891                H = Inequality_generic(A, b)
     892                self.ineqs_generic.append(H)
     893        for Hrep_obj in polyhedron.equation_generator():
     894            A, b = self._make_A_b(Hrep_obj, permutation)
     895            # add inequality
     896            try:
     897                H = Inequality_int(A, b, max_abs_coordinates)
     898                self.ineqs_int.append(H)
     899            except (OverflowError, ValueError):
     900                H = Inequality_generic(A, b)
     901                self.ineqs_generic.append(H)
     902            # add sign-reversed inequality
     903            A = [ -a for a in A ]
     904            b = -b
     905            try:
     906                H = Inequality_int(A, b, max_abs_coordinates)
     907                self.ineqs_int.append(H)
     908            except (OverflowError, ValueError):
     909                H = Inequality_generic(A, b)
     910                self.ineqs_generic.append(H)
     911       
     912    def prepare_next_to_inner_loop(self, p):
     913        r"""
     914        Peel off the next-to-inner loop.
     915       
     916        In the next-to-inner loop of :func:`rectangular_box_points`,
     917        we have to repeatedly evaluate `A x-A_0 x_0+b`. To speed up
     918        computation, we pre-evaluate
     919       
     920        .. math::
     921
     922            c = b + sum_{i=2} A_i x_i
     923
     924        and only compute `A x-A_0 x_0+b = A_1 x_1 +c \geq 0` in the
     925        next-to-inner loop.
     926
     927        INPUT:
     928
     929        - ``p`` -- the point coordinates. Only ``p[2:]`` coordinates
     930          are potentially used by this method.
     931
     932        EXAMPLES::
     933
     934             sage: from sage.geometry.integral_points import InequalityCollection, print_cache
     935             sage: P = Polyhedron(ieqs=[(2,3,7,11)])
     936             sage: ieq = InequalityCollection(P, Permutation([1,2,3]), [0]*3,[1]*3); ieq
     937             The collection of inequalities
     938             integer: (3, 7, 11) x + 2 >= 0
     939             sage: ieq.prepare_next_to_inner_loop([2,1,3])
     940             sage: ieq.prepare_inner_loop([2,1,3])
     941             sage: print_cache(ieq)
     942             Cached inner loop: 3 * x_0 + 42 >= 0
     943             Cached next-to-inner loop: 3 * x_0 + 7 * x_1 + 35 >= 0
     944        """
     945        for ineq in self.ineqs_int:
     946            (<Inequality_int>ineq).prepare_next_to_inner_loop(p)
     947        for ineq in self.ineqs_generic:
     948            (<Inequality_generic>ineq).prepare_next_to_inner_loop(p)
     949
     950    def prepare_inner_loop(self, p):
     951        r"""
     952        Peel off the inner loop.
     953       
     954        In the inner loop of :func:`rectangular_box_points`, we have
     955        to repeatedly evaluate `A x+b\geq 0`. To speed up computation, we pre-evaluate
     956       
     957        .. math::
     958
     959            c = A x - A_0 x_0 +b = b + sum_{i=1} A_i x_i
     960
     961        and only test `A_0 x_0 +c \geq 0` in the inner loop.
     962
     963        You must call :meth:`prepare_next_to_inner_loop` before
     964        calling this method.
     965
     966        INPUT:
     967
     968        - ``p`` -- the coordinates of the point to loop over. Only the
     969          ``p[1:]`` entries are used.
     970
     971        EXAMPLES::
     972       
     973             sage: from sage.geometry.integral_points import InequalityCollection, print_cache
     974             sage: P = Polyhedron(ieqs=[(2,3,7,11)])
     975             sage: ieq = InequalityCollection(P, Permutation([1,2,3]), [0]*3,[1]*3); ieq
     976             The collection of inequalities
     977             integer: (3, 7, 11) x + 2 >= 0
     978             sage: ieq.prepare_next_to_inner_loop([2,1,3])
     979             sage: ieq.prepare_inner_loop([2,1,3])
     980             sage: print_cache(ieq)
     981             Cached inner loop: 3 * x_0 + 42 >= 0
     982             Cached next-to-inner loop: 3 * x_0 + 7 * x_1 + 35 >= 0
     983        """
     984        for ineq in self.ineqs_int:
     985            (<Inequality_int>ineq).prepare_inner_loop(p)
     986        for ineq in self.ineqs_generic:
     987            (<Inequality_generic>ineq).prepare_inner_loop(p)
     988
     989    def swap_ineq_to_front(self, i):
     990        r"""
     991        Swap the ``i``-th entry of the list to the front of the list of inequalities.
     992
     993        INPUT:
     994
     995        - ``i`` -- Integer. The :class:`Inequality_int` to swap to the
     996          beginnig of the list of integral inequalities.
     997
     998        EXAMPLES::
     999
     1000            sage: from sage.geometry.integral_points import InequalityCollection
     1001            sage: P_QQ = Polyhedron(identity_matrix(3).columns() + [(-2, -1,-1)], field=QQ)
     1002            sage: iec = InequalityCollection(P_QQ, Permutation([1,2,3]), [0]*3,[1]*3)
     1003            sage: iec
     1004            The collection of inequalities
     1005            integer: (-1, -1, -1) x + 1 >= 0
     1006            integer: (-1, -1, 4) x + 1 >= 0
     1007            integer: (-1, 4, -1) x + 1 >= 0
     1008            integer: (3, -2, -2) x + 2 >= 0
     1009            sage: iec.swap_ineq_to_front(3)
     1010            sage: iec
     1011            The collection of inequalities
     1012            integer: (3, -2, -2) x + 2 >= 0
     1013            integer: (-1, -1, -1) x + 1 >= 0
     1014            integer: (-1, -1, 4) x + 1 >= 0
     1015            integer: (-1, 4, -1) x + 1 >= 0
     1016        """
     1017        i_th_entry = self.ineqs_int.pop(i)
     1018        self.ineqs_int.insert(0, i_th_entry)
     1019
     1020    def are_satisfied(self, inner_loop_variable):
     1021        r"""
     1022        Return whether all inequalities are satisfied.
     1023       
     1024        You must call :meth:`prepare_inner_loop` before calling this
     1025        method.
     1026
     1027        INPUT:
     1028
     1029        - ``inner_loop_variable`` -- Integer. the 0-th coordinate of
     1030          the lattice point.
     1031
     1032        OUTPUT:
     1033
     1034        Boolean. Whether the lattice point is in the polyhedron.
     1035
     1036        EXAMPLES::
     1037
     1038            sage: from sage.geometry.integral_points import InequalityCollection
     1039            sage: line = Polyhedron(eqns=[(2,3,7)])
     1040            sage: ieq = InequalityCollection(line, Permutation([1,2]), [0]*2, [1]*2 )
     1041            sage: ieq.prepare_next_to_inner_loop([3,4])
     1042            sage: ieq.prepare_inner_loop([3,4])
     1043            sage: ieq.are_satisfied(3)
     1044            False
     1045        """
     1046        cdef int i
     1047        for i in range(0,len(self.ineqs_int)):
     1048            ineq = self.ineqs_int[i]
     1049            if (<Inequality_int>ineq).is_not_satisfied(inner_loop_variable):
     1050                if i>0:
     1051                    self.swap_ineq_to_front(i)
     1052                return False
     1053        for i in range(0,len(self.ineqs_generic)):
     1054            ineq = self.ineqs_generic[i]
     1055            if (<Inequality_generic>ineq).is_not_satisfied(inner_loop_variable):
     1056                return False
     1057        return True
     1058   
     1059
     1060
     1061
     1062cpdef print_cache(InequalityCollection inequality_collection):
     1063    r"""
     1064    Print the cached values in :class:`Inequality_int` (for
     1065    debugging/doctesting only).
     1066
     1067    EXAMPLES::
     1068
     1069        sage: from sage.geometry.integral_points import InequalityCollection, print_cache
     1070        sage: P = Polyhedron(ieqs=[(2,3,7)])
     1071        sage: ieq = InequalityCollection(P, Permutation([1,2]), [0]*2,[1]*2); ieq
     1072        The collection of inequalities
     1073        integer: (3, 7) x + 2 >= 0
     1074        sage: ieq.prepare_next_to_inner_loop([3,5])
     1075        sage: ieq.prepare_inner_loop([3,5])
     1076        sage: print_cache(ieq)
     1077        Cached inner loop: 3 * x_0 + 37 >= 0
     1078        Cached next-to-inner loop: 3 * x_0 + 7 * x_1 + 2 >= 0
     1079    """
     1080    cdef Inequality_int ieq = <Inequality_int>(inequality_collection.ineqs_int[0])
     1081    print 'Cached inner loop: ' + \
     1082        str(ieq.coeff) + ' * x_0 + ' + str(ieq.cache) + ' >= 0'
     1083    print 'Cached next-to-inner loop: ' + \
     1084        str(ieq.coeff) + ' * x_0 + ' + \
     1085        str(ieq.coeff_next) + ' * x_1 + ' + str(ieq.cache_next) + ' >= 0'
     1086
     1087
  • sage/geometry/polyhedra.py

    diff --git a/sage/geometry/polyhedra.py b/sage/geometry/polyhedra.py
    a b  
    118118from sage.categories.objects import Objects
    119119
    120120from subprocess import Popen, PIPE
    121 from sage.misc.all import tmp_filename
     121from sage.misc.all import tmp_filename, cached_method, prod
    122122from sage.misc.functional import norm
    123123from sage.misc.package import is_package_installed
    124124
     
    12061206        """
    12071207        Returns a string representation of the vertex.
    12081208
     1209        OUTPUT:
     1210
     1211        String.
     1212
    12091213        TESTS::
    12101214
    12111215            sage: p = Polyhedron(ieqs = [[0,0,1],[0,1,0],[1,-1,0]])
     
    46464650        """
    46474651        if not self.is_compact():
    46484652            raise ValueError, 'Can only enumerate points in a compact polyhedron.'
    4649 
    46504653        lp = self.lattice_polytope(True)
    4651        
    46524654        # remove cached values to get accurate timings
    46534655        try:
    46544656            del lp._points
    46554657            del lp._npoints
    46564658        except AttributeError:
    46574659            pass
    4658 
    46594660        if self.is_lattice_polytope():
    46604661            return lp.points().columns()
    4661 
    46624662        points = filter(lambda p: self.contains(p),
    46634663                        lp.points().columns())
    46644664        return points
    46654665
    4666     def _integral_points_lattice_simplex(self, vertices=None):
     4666    @cached_method
     4667    def bounding_box(self, integral=False):
    46674668        r"""
    4668         Return the integral points in a lattice simplex.
    4669        
     4669        Return the coordinates of a rectangular box containing the non-empty polytope.
     4670
    46704671        INPUT:
    46714672
    4672         - ``vertices`` -- ``None`` (default) or a list of
    4673           integers. The indices of vertices that span the simplex
    4674           under consideration. By default, it is assumed that the
    4675           polytope itself is a simplex and all vertices are used.
    4676 
    4677         OUTPUT:
    4678 
    4679         A tuple containing the integral point coordinates as `\ZZ`-vectors.
    4680 
    4681         EXAMPLES::
     4673        - ``integral`` -- Boolean (default: ``False``). Whether to
     4674          only allow integral coordinates in the bounding box.
     4675       
     4676        OUTPUT:
     4677
     4678        A pair of tuples ``(box_min, box_max)`` where ``box_min`` are
     4679        the coordinates of a point bounding the coordinates of the
     4680        polytope from below and ``box_max`` bounds the coordinates
     4681        from above.
     4682
     4683        EXAMPLES::
     4684
     4685            sage: Polyhedron([ (1/3,2/3), (2/3, 1/3) ]).bounding_box()
     4686            ((1/3, 1/3), (2/3, 2/3))
     4687            sage: Polyhedron([ (1/3,2/3), (2/3, 1/3) ]).bounding_box(integral=True)
     4688            ((0, 0), (1, 1))
     4689            sage: polytopes.buckyball().bounding_box()       
     4690            ((-1059/1309, -1059/1309, -1059/1309), (1059/1309, 1059/1309, 1059/1309))
     4691        """
     4692        box_min = []
     4693        box_max = []
     4694        if self.n_vertices==0:
     4695            raise ValueError('Empty polytope is not allowed')
     4696        for i in range(0,self.ambient_dim()):
     4697            coords = [ v[i] for v in self.Vrep_generator() ]
     4698            max_coord = max(coords)
     4699            min_coord = min(coords)
     4700            if integral:
     4701                box_max.append(ceil(max_coord))
     4702                box_min.append(floor(min_coord))
     4703            else:
     4704                box_max.append(max_coord)
     4705                box_min.append(min_coord)
     4706        return (tuple(box_min), tuple(box_max))
     4707       
     4708    def integral_points(self, threshold=100000):
     4709        r"""
     4710        Return the integral points in the polyhedron.
     4711
     4712        Uses either the naive algorithm (iterate over a rectangular
     4713        bounding box) or triangulation + Smith form.
     4714
     4715        INPUT:
     4716
     4717        - ``threshold`` -- integer (default: 100000). Use the naive
     4718          algorith as long as the bounding box is smaller than this.
     4719
     4720        OUTPUT:
     4721       
     4722        The list of integral points in the polyhedron. If the
     4723        polyhedron is not compact, a ``ValueError`` is raised.
     4724
     4725        EXAMPLES::
     4726       
     4727            sage: Polyhedron(vertices=[(-1,-1),(1,0),(1,1),(0,1)]).integral_points()
     4728            ((-1, -1), (0, 0), (0, 1), (1, 0), (1, 1))
    46824729
    46834730            sage: simplex = Polyhedron([(1,2,3), (2,3,7), (-2,-3,-11)])
    4684             sage: simplex._integral_points_lattice_simplex()
     4731            sage: simplex.integral_points()
    46854732            ((-2, -3, -11), (0, 0, -2), (1, 2, 3), (2, 3, 7))
    46864733
    4687         The simplex need not be full-dimensional::
     4734        The polyhedron need not be full-dimensional::
    46884735
    46894736            sage: simplex = Polyhedron([(1,2,3,5), (2,3,7,5), (-2,-3,-11,5)])
    4690             sage: simplex._integral_points_lattice_simplex()
     4737            sage: simplex.integral_points()
    46914738            ((-2, -3, -11, 5), (0, 0, -2, 5), (1, 2, 3, 5), (2, 3, 7, 5))
    46924739
    46934740            sage: point = Polyhedron([(2,3,7)])
    4694             sage: point._integral_points_lattice_simplex()
     4741            sage: point.integral_points()
    46954742            ((2, 3, 7),)
    46964743
    4697         TESTS::
    4698        
     4744            sage: empty = Polyhedron()
     4745            sage: empty.integral_points()
     4746            ()
     4747
     4748        Here is a simplex where the naive algorithm of running over
     4749        all points in a rectangular bounding box no longer works fast
     4750        enough::
     4751
    46994752            sage: v = [(1,0,7,-1), (-2,-2,4,-3), (-1,-1,-1,4), (2,9,0,-5), (-2,-1,5,1)]
    47004753            sage: simplex = Polyhedron(v); simplex
    47014754            A 4-dimensional polyhedron in QQ^4 defined as the convex hull of 5 vertices.
    4702             sage: len(simplex._integral_points_lattice_simplex())
     4755            sage: len(simplex.integral_points())
    47034756            49
    47044757
    4705             sage: v = [(4,-1,-1,-1), (-1,4,-1,-1), (-1,-1,4,-1), (-1,-1,-1,4), (-1,-1,-1,-1)]
    4706             sage: P4mirror = Polyhedron(v)
    4707             sage: len( P4mirror._integral_points_lattice_simplex() )
    4708             126
    4709         """
    4710         if not vertices:
    4711             assert self.is_simplex() and self.is_lattice_polytope()
    4712             vertices = [ vector(ZZ, v) for v in self.Vrepresentation() ]
    4713         else:
    4714             vertices = [ vector(ZZ, self.Vrepresentation(i)) for i in vertices ]
    4715            
    4716         origin = vertices.pop()
    4717         origin.set_immutable()
    4718         if len(vertices)==0:
    4719             return tuple([origin])
    4720         rays = [ v-origin for v in vertices ]
    4721 
    4722         # Find equation Ax<=b that cuts out simplex from parallelotope
    4723         R = matrix(rays)
    4724         if R.is_square():
    4725             b = 1
    4726             A = R.solve_right(vector([1]*len(rays)))
    4727         else:
    4728             b = 1
    4729             RRt = R * R.transpose()
    4730             A = RRt.solve_right(vector([1]*len(rays))) * R
    4731        
    4732         from sage.geometry.cone import parallelotope_points
    4733         gens = list(parallelotope_points(rays)) + list(rays)
    4734         gens = [ origin+x for x in
    4735                  list(parallelotope_points(rays)) + list(rays)
    4736                  if A*x<=b ]
    4737         for x in gens:
    4738             x.set_immutable()
    4739         return tuple(gens)
    4740 
    4741     def integral_points(self):
    4742         r"""
    4743         Return the integral points in the polyhedron.
    4744 
    4745         OUTPUT:
    4746        
    4747         The list of integral points in the polyhedron. If the
    4748         polyhedron is not compact, a ``ValueError`` is raised.
    4749 
    4750         EXAMPLES::
    4751        
    4752             sage: Polyhedron(vertices=[(-1,-1),(1,0),(1,1),(0,1)]).integral_points()
    4753             ((0, 1), (1, 0), (0, 0), (1, 1), (-1, -1))
    4754 
    4755         Finally, 3-d reflexive polytope number 4078::
     4758        Finally, the 3-d reflexive polytope number 4078::
    47564759
    47574760            sage: v = [(1,0,0), (0,1,0), (0,0,1), (0,0,-1), (0,-2,1),
    47584761            ...        (-1,2,-1), (-1,2,-2), (-1,1,-2), (-1,-1,2), (-1,-3,2)]
    4759             sage: pts1 = set(Polyhedron(v).integral_points())    # Sage's own code
     4762            sage: P = Polyhedron(v)
     4763            sage: pts1 = P.integral_points()                     # Sage's own code
     4764            sage: all(P.contains(p) for p in pts1)
     4765            True
    47604766            sage: pts2 = LatticePolytope(v).points().columns()   # PALP
    4761             sage: for p in pts2:
    4762             ...       p.set_immutable()
    4763             sage: pts2 = set(pts2)
    4764             sage: pts1 == pts2
     4767            sage: for p in pts1: p.set_immutable()
     4768            sage: for p in pts2: p.set_immutable()
     4769            sage: set(pts1) == set(pts2)
    47654770            True
    47664771
    4767         TODO: profile and make faster::
    4768 
    47694772            sage: timeit('Polyhedron(v).integral_points()')   # random output
    47704773            sage: timeit('LatticePolytope(v).points()')       # random output
    47714774        """
    47724775        if not self.is_compact():
    47734776            raise ValueError('Can only enumerate points in a compact polyhedron.')
    4774         if not self.is_lattice_polytope():
    4775             raise NotImplementedError('Only implemented for polyhedra with integral vertices')
    4776        
     4777        if self.n_vertices() == 0:
     4778            return tuple()
     4779
     4780        # for small bounding boxes, it is faster to naively iterate over the points of the box
     4781        box_min, box_max = self.bounding_box(integral=True)
     4782        box_points = prod(max_coord-min_coord+1 for min_coord, max_coord in zip(box_min, box_max))
     4783        if  not self.is_lattice_polytope() or \
     4784                (self.is_simplex() and box_points<1000) or \
     4785                box_points<threshold:
     4786            from sage.geometry.integral_points import rectangular_box_points
     4787            return rectangular_box_points(box_min, box_max, self)
     4788
     4789        # for more complicate polytopes, triangulate & use smith normal form
     4790        from sage.geometry.integral_points import simplex_points
    47774791        if self.is_simplex():
    4778             points = self._integral_points_lattice_simplex()
    4779         else:
    4780             triangulation = self.triangulate()
    4781             points = set()
    4782             for simplex in triangulation:
    4783                 points.update(self._integral_points_lattice_simplex(simplex))
    4784             points = tuple(points)
    4785                
     4792            return simplex_points(self.Vrepresentation())
     4793        triangulation = self.triangulate()
     4794        points = set()
     4795        for simplex in triangulation:
     4796            triang_vertices = [ self.Vrepresentation(i) for i in simplex ]
     4797            new_points = simplex_points(triang_vertices)
     4798            for p in new_points:
     4799                p.set_immutable()
     4800            points.update(new_points)
    47864801        # assert all(self.contains(p) for p in points)   # slow
    4787         return points
     4802        return tuple(points)
    47884803   
    47894804    def combinatorial_automorphism_group(self):
    47904805        """