This post is my attempt to answer the question from   here : how to find all integer points (all points with integer coordinates) in the intersection of two cubes. The following procedure  IntegerPoints  solves a more general problem: it finds all the integer points of a bounded polyhedral region of arbitrary dimension, defined by a system of linear inequalities and / or equations.

Required parameters of the procedure: SN is a set or a list of linear inequalities and/or equations with any number of variables, the Var is the list of variables. The procedure returns the set of all integer points, satisfying the conditions  SN .

Code of the procedure:

restart;

IntegerPoints := proc (SN::{list, set}, Var::list)

local SN1, sn, n, Sol, k, i, s, S, R;

uses PolyhedralSets, SolveTools[Inequality];

SN1 := convert(evalf(SN), fraction);

for sn in SN1 do

if type(sn, `<`) then SN1 := subs(sn = (`<=`(op(sn))), SN1)

end if; end do;

if IsBounded(PolyhedralSet(SN1)) = false then error "The region should be bounded" end if;

n := nops(Var);

Sol := LinearMultivariateSystem(SN, Var);

if Sol = {} then return {} else

k := 0;

for s in Sol do if nops(indets(s[1])) = 1 then

S[0] := [[]];

for i to n do

S[i] := [seq(seq([op(j1), op(j2)], j2 = [isolve(eval(s[i], j1))]), j1 = S[i-1])] end do;

k := k+1; R[k] := op(S[n]);

end if; end do;

convert(R, set);

map(t->rhs~(t), %);

end if;

end proc:

Examples of use:

IntegerPoints({x > 0, y > 0, z > 0, 2*x+3*y+z < 12}, [x, y, z]);

  {[1, 1, 1], [1, 1, 2], [1, 1, 3], [1, 1, 4], [1, 1, 5], [1, 1, 6], [1, 2, 1], [1, 2, 2], [1, 2, 3], [2, 1, 1], [2, 1, 2],

                                   [2, 1, 3], [2, 1, 4], [2, 2, 1], [3, 1, 1], [3, 1, 2]}

IntegerPoints({x > 0, y > 0, z > 0, 2*x+3*y+z = 12}, [x, y, z]);

                                    {[1, 1, 7], [1, 2, 4], [1, 3, 1], [2, 1, 5], [2, 2, 2], [3, 1, 3], [4, 1, 1]}

IntegerPoints([x > 0, y > 0, z > 0, 2*x+3*y+z = 12, x+y+z <= 6], [x, y, z]);

                                                           {[1, 3, 1], [2, 2, 2], [4, 1, 1]}

isolve({x > 0, y > 0, z > 0, 2*x+3*y+z < 12});  #  isolve fails with these examples

              Warning, solutions may have been lost

isolve({x > 0, y > 0, z > 0, 2*x+3*y+z = 12});

              Warning, solutions may have been lost

In the following example (with a visualization) we find all integer point in the intersection of a square and a triangle:

S1 := {x > 0, y > 0, x < 13/2, y < 13/2}:

S2 := {y > (1/4)*x+1, y < 2*x, y+x < 12}:

S := IntegerPoints(`union`(S1, S2), [x, y]):

Region := plots[inequal](`union`(S1, S2), x = 0 .. 7, y = 0 .. 7, color = "LightGreen", nolines):

Points := plot([op(S)], style = point, color = red, symbol = solidcircle):

Square := plottools[curve]([[0, 0], [13/2, 0], [13/2, 13/2], [0, 13/2], [0, 0]], color = blue, thickness = 3):

Triangle := plottools[curve]([[4/7, 8/7], [4, 8], [44/5, 16/5], [4/7, 8/7]], color = blue, thickness = 3):

plots[display](Square, Triangle, Points, Region, scaling = constrained);

In the following example (with a visualization) we find all integer point in the intersection of two cubes. The second cube is obtained from the first cube by rotation with orthogonal matrix  A  and by a translation:

A := <1/3, 2/3, 2/3; -2/3, 2/3, -1/3; -2/3, -1/3, 2/3>:

f := unapply(A^(-1).<x+5, y-4, z-7>, x, y, z):

S1 := {x > 0, y > 0, z > 0, x < 6, y < 6, z < 6}:

S2 := eval(S1, {x = f(x, y, z)[1], y = f(x, y, z)[2], z = f(x, y, z)[3]}):

S := IntegerPoints(`union`(S1, S2), [x, y, z]);

Points := plots[pointplot3d](S, color = red, symbol = box):

Cube := plottools[cuboid]([0, 0, 0], [6, 6, 6], color = blue, linestyle = solid):

F := plottools[transform]((x, y, z)->convert(A.<x, y, z>+<-5, 4, 7>, list)):

plots[display](Cube,  F(Cube), Points, scaling = constrained, linestyle = solid, transparency = 0.7, orientation = [25, 75], axes = normal);

In the example below, all the ways to exchange $ 1 coins of 1, 5, 10, 25 and 50 cents, if the number of coins no more than 8, there is no pennies and there is at least one 50-cent coin:

IntegerPoints({x1 = 0, x2 >= 0, x3 >= 0, x4 >= 0, x5 >= 1,  x1+5*x2+10*x3+25*x4+50*x5 = 100, x1+x2+x3+x4+x5 <= 8}, [x1, x2, x3, x4, x5]);

nops(%);

                              {[0, 0, 0, 0, 2], [0, 0, 0, 2, 1], [0, 0, 5, 0, 1], [0, 1, 2, 1, 1], [0, 2, 4, 0, 1],

                                                 [0, 3, 1, 1, 1], [0, 4, 3, 0, 1], [0, 5, 0, 1, 1]}

                                                                                    8

Integer_points.mw

Addition: Below in my comments another procedure  IntegerPoints1  is presented that solves the same problem.