How to find lattice points on a line segment?

14

2

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

Saket Kumar

Posted 2012-12-20T08:38:13.227

Reputation: 143

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 – 2012-12-20T09:32:42.487

Do you want all the lattice points on the line, or just those on the line segment between the two given points? – murray – 2012-12-20T15:51:39.843

Answers

16

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]

Mathematica graphics

Yves Klett

Posted 2012-12-20T08:38:13.227

Reputation: 14 743

27

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

Artes

Posted 2012-12-20T08:38:13.227

Reputation: 51 831

gotta love FrobeniusSolve – cartonn – 2012-12-20T21:02:59.823

@cartonn This cannot be overcome by anything else in such cases. – Artes – 2012-12-20T21:08:14.470

2+1 for casting the problem in a way that never occurred to me. – Mr.Wizard – 2012-12-21T08:03:03.287

1@Mr.Wizard Thanks, last answers are usually underestimated. – Artes – 2012-12-21T12:04:59.220

13

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

enter image description here

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

kglr

Posted 2012-12-20T08:38:13.227

Reputation: 302 076

+1, never used that before... very useful and another proof that we live to learn. – Yves Klett – 2012-12-20T10:53:52.903

@Yves thanks for the vote. Learned about it just yesterday while struggling with Interpolation :) – kglr – 2012-12-20T11:01:31.450

Nice idea. Incidentally, I noticed, in editing the question, that it stipulates the solutions should have positive coordinates. – whuber – 2012-12-20T16:28:26.967

@whuber, good point - funny I missed 50% of the requirements in a single-line question:) I will update with the added constraints. – kglr – 2012-12-20T16:41:56.493

@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 – 2012-12-20T20:46:23.843

@artes, this sounded like a trick question; so it felt wiser to reply "Of course not ...". But, out of curiousity, I tried the inputs you suggested, and got {{1, 10310}, {6, 10307}, {11, 10304}, {16, 10301}} (almost, in no time). So, then ... "why not?" (or else, I am missing something?) – kglr – 2012-12-20T22:30:57.113

@artes, ... then again, I would not suggest any of these against FrobeniusSolve in a speed and/or elegance contest :) (Btw, I need to train my memory to think of FrobeniusSolve first when I hear the keywords "integer, positive, solution.") – kglr – 2012-12-20T22:41:13.347

@kguler It is not surprising that FindInstance found something, but it doesn't says how many solutions there are. If there are none, you can't be sure working with it. Perhaps I'm prejudiced against it, nevertheless I can't find it very helpful or I'haven't seen really nice examples ot its usage. And at last I understand this task to find all points and when there are thousands it is not constructive. This was my only objection. – Artes – 2012-12-20T22:52:38.160

@artes, agreed. – kglr – 2012-12-20T22:55:24.357

@kguler I haven't tested your solution but shouldn't there be y > 0 instead of y <= 0 ? – Artes – 2012-12-20T22:58:06.420

@Artes, yes ... Thank you. – kglr – 2012-12-20T23:01:10.340

5

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

J. M.'s ennui

Posted 2012-12-20T08:38:13.227

Reputation: 115 520

3

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

enter image description here

BoLe

Posted 2012-12-20T08:38:13.227

Reputation: 5 649

2

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

Aky

Posted 2012-12-20T08:38:13.227

Reputation: 2 547