Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am trying to generate list of coordinates that forms a subset in the integer lattice $\mathbb{Z}^n$. This set is given by set of inequalities. Right now, I am using the following:

baseLattice = Tuples[Range[0, k], {n}];
coords = DeleteCases[DeleteCases[baseGrid, 
    x_ /; x[[2]] > (1.9*x[[1]] + 3)],
    x_ /; x[[2]] < (1.9*x[[1]] - 4)];

This approach has the obvious disadvantage that for larger n or k the Tuples command generates all points in the first quadrant. I realize could also build my set in some loop, however I feel like there could be a straightforward solution which I have not yet figured out.

For illustration, I am trying to get a set of coordinates similar to the set of red dots in this picture:

Example of set of points given by two lines of same slope

I did not post my original problem in full in the OP. For the slopes of the lines, I am using numerical values of algebraic numbers. I used the following code to compare performance of suggested brute-force solutions:

Solve by Daniel Lichtblau:

f = #^2 - # - 1 &;
q = AlgebraicNumber[Root[f, 2], {0, 1}];
q = N[q];
coords1 = {x, y} /. 
  Solve[{0 <= x <= y - 1, 0 <= y <= 100, q*x - q^2 <= y <= q* x + 2*q}, {x, y}, Integers]

And vectorization by Szabolcz:

lattice = Tuples[Range[0, 100], 2];
{x, y} = Transpose[lattice];
coords2 = 
   Pick[lattice, UnitStep[q*x + 2*q - y] UnitStep[y - (q*x - 2 q^2)], 1];

I run the code with k = 100 and n = 2. For q = N[algebraic number]:
Solve ... 50.11s
Vectorization ... 0.02s
My original way ... 0.17s.

Turns out this is somewhat unfair to solve, as when I take q = algebraic number itself (I guess not surprisingly):
Solve ... 50.28s
Vectorization ... 81.43s
My original way ... 78.91s.

share|improve this question
To generate all coordinates, use Tuples instead of Subsets. This will still be slow though, and also memory intensive. Finally try to vectorize the inequality testing. –  Szabolcs Feb 10 '14 at 22:26
coords = {x, y} /. Solve[{0 <= x <= y - 1, 0 <= y <= 100, 19/10*x - 4 <= y <= 19/10 x + 3}, {x, y}, Integers]; –  Daniel Lichtblau Feb 10 '14 at 22:35
I tried using "the way of solve". It is a neat and I did not think of it, but it is very slow on my machine, therefore it is not useful for me. –  BoZenKhaa Feb 10 '14 at 22:56
What example(s) did you try it on wherein brute force was notably better? –  Daniel Lichtblau Feb 11 '14 at 2:55
(1) Thanks. I now see what are some of the issues. (1) If you do RootReduce or ToRadicals and use that as q then it seems faster. Around 2 seconds for me. The problem with using N[q] is that q is as far from any rational as possible (if you know what that means then fine, else ignore it). –  Daniel Lichtblau Feb 12 '14 at 18:34

2 Answers 2

up vote 1 down vote accepted

Here's a way that might be more computationally tractable:

lin[x_] := (x[[2]] > (1.9*x[[1]] - 4)) && (x[[2]] < (1.9*x[[1]] + 3)); 
line = Table[{i, 1.9 i + 0.5}, {i, -35, 35}];
points = Flatten[DeleteCases[
   Table[If[lin[Round[line[[i]] + {0, j}]], Round[line[[i]] + {0, j}]], 
         {i, 1, Length[line]}, {j, -50, 50}], Null, 2], 1]

This builds a center line (sort of the average of the two lines) and then tests points to the left and right, rejecting them if they are outside the bounds and accepting them if within.

share|improve this answer

Here's a more efficient way. Note that this way is not any better algorithmically and it's still going to be very slow for large problems! It will still brute force check everything. It just makes a better use of Mathematica's functions to try to be faster ...

k = 100;
n = 2;

First, use Tuples to get the complete lattice, not just the upper half.

lattice = Tuples[Range[0, k], {n}];

Then use vectorization for the testing.

Suppose the inequalities are $1.9 x+ 3 \ge y$ and $1.9x-4 \le y$. Then do

{x, y} = Transpose[lattice];
points = Pick[lattice, UnitStep[1.9 x + 3 - y] UnitStep[y - (1.9 x - 4)], 1];

to obtain the points which pass the test.

I was lazy with the inequalities and I used $\ge$ instead of $>$. You'll need to work extra to get this right with UnitStep. Note that UnitStep[0] is 1.


Using the tools I describe here, you can do it as

lattice = Tuples[Range[0, k], {n}];
{x, y} = Transpose[lattice];

<< BoolEval`

ListPlot[{lattice, BoolPick[lattice, 1.9 x + 3 >= y >= 1.9 x - 4]}, PlotStyle -> {Black, Red}]

enter image description here

share|improve this answer
Thank you, even though it is still brute forcing the solution, it is incomparable in terms of speed to the code I was using. I take it there is probably no built-in function for generating a subset of Tuples[] that would cut the number of tuples down? –  BoZenKhaa Feb 10 '14 at 23:05
@BoZenKhaa I don't know of any. –  Szabolcs Feb 11 '14 at 20:00

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.