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.

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

share|improve this question
    
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

6 Answers 6

up vote 15 down vote accepted

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

share|improve this answer

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}}
share|improve this answer
    
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
2  
+1 for casting the problem in a way that never occurred to me. –  Mr.Wizard Dec 21 '12 at 8:03
1  
@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}}]

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}]
share|improve this answer
    
+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 :) –  kguler 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. –  kguler 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 /.
    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}}
share|improve this answer
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. ;)

share|improve this answer

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

share|improve this answer

Your Answer

 
discard

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.