# How to find lattice points on a line segment?

How do I find points on the line segment joining {-4, 11} and {16, -1} whose coordinates are positive integers?

-
Thanks for the accept - you could also wait a few days to encourage different answers and then choose the best one. You can also change the accept at any time (not that I mind if you keep it where it is). – Yves Klett Dec 20 '12 at 9:32
Do you want all the lattice points on the line, or just those on the line segment between the two given points? – murray Dec 20 '12 at 15:51

Is this what you are searching for?

a = {-4, 11};
b = {16, -1};

dy = (b[[2]] - a[[2]])/(b[[1]] - a[[1]]);

offset = u /. Solve[a[[2]] == dy*a[[1]] + u, u][[1]];

coords = {x,
y} /. {Reduce[y == dy*x + offset && x > 0 && y > 0, {x, y},
Integers] // ToRules}


(* {{1, 8}, {6, 5}, {11, 2}} *)

Graphics[{PointSize[Large], Point[{a, b}], Red, Point[coords],
Line[{a, b}]}, Axes -> True, GridLines -> {Range[16], Range[16]},
ImageSize -> 640]


-

There are many ways to proceed, the best one uses FrobeniusSolve :

I

Since we know, that

a x + b == y /. Solve[{-4 a + b == 11, 16 a + b == -1}, {a, b}] // Simplify

{3 x + 5 y == 43}


we find

FrobeniusSolve[ {3, 5}, 43]

{{1, 8}, {6, 5}, {11, 2}}


a bit more straightforward way :

II

{x, y} /. Solve[ (a x + b == y /. Solve[ {-4 a + b == 11, 16 a + b == -1}, {a, b}])
~  Join  ~  {x > 0, y > 0}, {x, y}, Integers]

{{1, 8}, {6, 5}, {11, 2}}

-
gotta love FrobeniusSolve – cartonn Dec 20 '12 at 21:02
@cartonn This cannot be overcome by anything else in such cases. – Artes Dec 20 '12 at 21:08
+1 for casting the problem in a way that never occurred to me. – Mr.Wizard Dec 21 '12 at 8:03
@Mr.Wizard Thanks, last answers are usually underestimated. – Artes Dec 21 '12 at 12:04

You can also use InterpolatingPolynomial with Solve, Reduce or Eliminate:

  a = {-4, 11}; b = {16, -1};
coords =  Solve[y == InterpolatingPolynomial[{a, b}, x] && 0 <= x <= 16&&0<=y,
{x, y}, Integers][[All, All, 2]];
(* or *)
coords={ToRules[Reduce[ y == InterpolatingPolynomial[{a, b}, x] &&
0 <= x <= 16&&0<=y,  {x, y}, Integers]]}[[All, All, 2]];
(* or *)
coords = FindInstance[y == InterpolatingPolynomial[{a, b}, x] && 0 <= x <= 16&&0<=y,
{x, y},  Integers, 5][[All, All, 2]]


All three give

  { {1, 8}, {6, 5}, {11, 2}}


To show in a plot:

  Plot[InterpolatingPolynomial[{a, b}, x], {x, -5, 17},
Mesh -> {First /@ coords}, MeshStyle -> PointSize[Large],
PlotRange -> {{-5, 20}, {-2, 15}}]


Update: You can also use the plain old Interpolation in all of the above. For example,

  FindInstance[ y == Quiet@Interpolation[{a, b}][x] && 0 <= x <= 16 && 0 <= y,
{x, y}, Integers, 5][[All, All, 2]]
(* {{1, 8}, {6, 5}, {11, 2}} *)


Update 2: Getting into Cases, Select, Pick ... territory:

  Cases[{#, Interpolation[{a, b}, InterpolationOrder -> 1][#]} & /@
Range[0, 16], {_Integer, _Integer?Positive}]


or

 Cases[{#, InterpolatingPolynomial[{a, b}, #]} & /@
Range[0, 16], {_Integer, _Integer?Positive}]

-
+1, never used that before... very useful and another proof that we live to learn. – Yves Klett Dec 20 '12 at 10:53
@Yves thanks for the vote. Learned about it just yesterday while struggling with Interpolation :) – kglr Dec 20 '12 at 11:01
Nice idea. Incidentally, I noticed, in editing the question, that it stipulates the solutions should have positive coordinates. – whuber Dec 20 '12 at 16:28
@whuber, good point - funny I missed 50% of the requirements in a single-line question:) I will update with the added constraints. – kglr Dec 20 '12 at 16:41
@kguler FindInstance is not a way to go. Could you really solve this problem having initially e.g. these points {{-4, 10313}, {16, 10301}} ? – Artes Dec 20 '12 at 20:46

Artes's solution is the best, I think. If you just want to treat this as an ordinary Diophantine problem, you can do that with Solve[] (making this approach more or less equivalent to Yves's):

{p, q} = {-4, 11};
{r, s} = {16, -1};
{x, y} /. Solve[{(q - s) x - (p - r) y == -Det[{{p, q}, {r, s}}],
x > 0, y > 0, Min[p, r] < x < Max[p, r], Min[q, s] < y < Max[q, s]},
{x, y}, Integers]
{{1, 8}, {6, 5}, {11, 2}}


One could also choose to use Bézout's identity to solve this problem (see for instance this excellent math.SE post by Arturo Magidin).

Luckily, ExtendedGCD[] is a built-in function for performing the extended Euclidean algorithm, so let's use that:

{g, v} = ExtendedGCD[q - s, p - r]
{4, {2, 1}}


We check something first:

w = -Det[{{p, q}, {r, s}}];
Divisible[w, g]
True


So a particular solution is then given by

f = w v/g
{86, -43}


We can derive a parametrized set of solutions like so:

sols[k_] = Simplify[f + (k - Max[Quotient[w v, {p - r, q - s}]]) {p - r, q - s}/g]
{11 - 5 k, 2 + 3 k}


As it turns out, sols[0] gives one of the needed solutions, and stepping forward (i.e. sols[1] and sols[2]) gives the others. If you're lazy, however, then you can use FindInstance[]:

Map[sols, k /.
Thread[Min /@ {{p, q}, {r, s}} <= sols[k] <= Max /@ {{p, q}, {r, s}}],
k, Integers, 3]]
{{11, 2}, {6, 5}, {1, 8}}

-
Quiet[Cases[Outer[List, Range[-4, 11], Range[16, -1, -1]],
{x_, y_} /; (y - 11)/(x + 4) == (y + 1)/(x - 16), {2}]]


This solution shows how to transform linear complexity to quadratic, and provides some relief of the comic variety. ;)

-

Suppose we know the equation of line through the two points, one can generate all points on the line with integer x and of them keep those with integer y. Without invoking solving function.

With[{x1 = -4, y1 = 11, x2 = 16, y2 = -1},
Table[{x, (y2 - y1)/(x2 - x1) (x - x1) + y1},
{x, x1, x2}]] // Cases[#, {_, _Integer}] &

(* {{-4, 11}, {1, 8}, {6, 5}, {11, 2}, {16, -1}} *)


Put this in Manipulate and you have an interactive canvas, showing the points on a segment between the end points which can me moved around the lattice.

Manipulate[
DynamicModule[{x1, y1, x2, y2, pts},
{x1, y1} = Round@p1;
{x2, y2} = Round@p2;
pts = Table[{x, (y2 - y1)/(x2 - x1) (x - x1) + y1},
{x, x1, x2}] // Cases[#, {_, _Integer}] &;
Graphics[{
Gray, Line[{{x1, y1}, {x2, y2}}],
Black, PointSize[.01], Point[pts]},
GridLines -> {Range[-10, 20], Range[-5, 15]},
GridLinesStyle -> LightGray,
AspectRatio -> Automatic,
Frame -> True,
ImageSize -> 500,
PlotRange -> {{-10, 20}, {-5, 15}}]],
{{p1, {-4, 11}}, Locator},
{{p2, {16, -1}}, Locator}]


-