14

2

How do I find points on the line segment joining `{-4, 11}`

and `{16, -1}`

whose coordinates are positive integers?

14

2

How do I find points on the line segment joining `{-4, 11}`

and `{16, -1}`

whose coordinates are positive integers?

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

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

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

**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 – 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}}
```

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

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

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