How to find lattice points on a line segment?

Question

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

Answer 1

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]

Answer 2

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}}

Answer 3

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}]

Answer 4

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 /.
    FindInstance[Thread[sols[k] > 0] ~Join~
                 Thread[Min /@ {{p, q}, {r, s}} <= sols[k] <= Max /@ {{p, q}, {r, s}}],
                 k, Integers, 3]]
   {{11, 2}, {6, 5}, {1, 8}}

Answer 5

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}]

Answer 6

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. ;)