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.