How do I find points on the line segment joining {-4, 11}
and {16, -1}
whose coordinates are positive integers?
-
$\begingroup$ 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). $\endgroup$ – Yves Klett Dec 20 '12 at 9:32
-
$\begingroup$ Do you want all the lattice points on the line, or just those on the line segment between the two given points? $\endgroup$ – 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}}
-
-
$\begingroup$ @cartonn This cannot be overcome by anything else in such cases. $\endgroup$ – Artes Dec 20 '12 at 21:08
-
2$\begingroup$ +1 for casting the problem in a way that never occurred to me. $\endgroup$ – Mr.Wizard Dec 21 '12 at 8:03
-
1$\begingroup$ @Mr.Wizard Thanks, last answers are usually underestimated. $\endgroup$ – 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}]
-
$\begingroup$ +1, never used that before... very useful and another proof that we live to learn. $\endgroup$ – Yves Klett Dec 20 '12 at 10:53
-
$\begingroup$ @Yves thanks for the vote. Learned about it just yesterday while struggling with
Interpolation
:) $\endgroup$ – kglr Dec 20 '12 at 11:01 -
$\begingroup$ Nice idea. Incidentally, I noticed, in editing the question, that it stipulates the solutions should have positive coordinates. $\endgroup$ – whuber Dec 20 '12 at 16:28
-
$\begingroup$ @whuber, good point - funny I missed 50% of the requirements in a single-line question:) I will update with the added constraints. $\endgroup$ – kglr Dec 20 '12 at 16:41
-
$\begingroup$ @kguler
FindInstance
is not a way to go. Could you really solve this problem having initially e.g. these points{{-4, 10313}, {16, 10301}}
? $\endgroup$ – 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}}
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}]
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. ;)