How to find a closest integer point to the intersection of two lines?

Question

Here's a question that originates from StackOverflow.

Given are two lines on a plane, specified by equations ($a x + b y = c$) with integer coefficients. The lines aren't parallel and they don't necessarily pass through any integer points. Given also is an integer point $(x,y)$ (x and y are both in $\mathbb{Z}$) that lies on neither line.

The problem is to find the integer point $(x',y')$ closest to the intersection of the lines that lies in the same quarter as $(x,y)$ (as separated by the lines). It may be $(x,y)$ itself, but we're interested in the nontrivial case when an even closer point exists.

Other than being a way to specify the quarter of interest, the (x,y) point is seemingly needed to ensure the problem is NP or easier. Without this point, it seems, the answer might not be polynomial wrt the line parameters $a_1$, $a_2$, $b_1$ etc.

The problem is claimed to have a polynomial solution, but I doubt that it really exists. Is there a solution or a proof of its NP-completeness?

Answer 1

If the point $(x,y)$ is in an obtuse quadrant between the lines the problem is easily solved by enumerating the lattice points in a closed sphere of radius $\sqrt{2}$ about the intersection point so I will only consider the case that $(x,y)$ lies in a acute quadrant.

The problem can be converted to one that is very similar to inhomogeneous Diophantine approximation. I know of at least one algorithm that will definitely find the correct answer in $O\left(\tfrac{1}{\theta}\right)$ operations where $2\theta$ is the acute angle between the lines. The algorithm is a minor modification of that given on page 19 of Vaughan Clarkson's thesis. I am very confident that Cassels' algorithm (or a minor modification of it) will solve the problem in $O\left(1 + \log\tfrac{1}{\theta}\right)$ operations (page 34 of 1), but I am not quite sure how to show it, or what modification needs to be made. Before I describe this, you need to know a little bit about inhomogeneous Diophantine approximation.

Let $\alpha$ and $\beta$ be real numbers. Define the function

$F(p,q) = |\alpha p - q - \beta|$.

The problem of inhomogeneous Diophantine approximation involves minimising $F(p,q)$ over integers $p$ and $q$ where $q$ is positive. The pair $(p, q)$ is called a best approximation if $F(p,q) < F(p',q')$ for all $q' < q$. The best approximations describe the minima found as $q$ and the magnitude of $p$ increase. There are various algorithms than can enumerate all of the best approximations. Two examples are the 'naive algorithm' (page 19 of 1) and Cassel's algorithm (page 34 of 1). The OP's problem is not exactly the same as this, but it is so similar that the algorithms (at least the naive algorithm) carry over.

Let our two lines be $\ell_1$ and $\ell_2$ and let $2\theta < \pi/2$ be the angle between them. It will be easier to describe the approach if we set the intersection of these lines to be at the origin and we look for the nearest point in the translated lattice $\mathbb{Z}^2 + t$ where $t$ is the appropriate translation. Define $\ell$ to be the unique line that passes through the origin and bisects the acute angle between $\ell_1$ and $\ell_2$. The angle between $\ell$ and $\ell_1$ and $\ell$ and $\ell_2$ is $\theta$. The problem can now be stated as:

Find the point $x \in \mathbb{Z}^2 + > t$ that is nearest to the origin such that the angle between $x$ and $\ell$ is less than or equal to $\theta$

Let $A(x)$ denote the angle between $\ell$ and $x$. Our motivation is now very similar to that in Diophantine approximation. That is, find all of the best approximations for $A(x)$, our problem is solved by the best approximation that first yields $A(x) \leq \theta$. It so happens that $A(x)$ is a very similar function to $F(p,q)$. To give this some context I will say that $F(p,q)$ is, in a sense, computing an inner product between two vectors, whereas $A(x)$ is computing an angle. In this context it is not surprising that the algorithms for Diophantine approximation can be used.

I will only consider the 'naive algorithm' and I'll just give some geometric insight as to its functionality, this should be enough to convince most people. Working through this rigorously is really beyond a typed answer on MO, but all the required machinery is in 1. The 'naive algorithm' enumerates every point in $\mathbb{Z}^2 + t$ that is a nearest lattice point to any point in the line $\ell$. In other words it consecutively locates (starting from the origin) every lattice point in $\mathbb{Z}^2 + t$ whose Voronoi cell (in this case squares) intersect $\ell$. A picture might be useful

It is not difficult to devise an algorithm which does this, just start at the origin and check where $\ell$ next crosses a boundary of a Voronoi cell. It is also easy to see that the points it locates are a super set of the best approximations for $A(x)$. The first point that the algorithm finds such that $A(x) < \theta$ is the solution to the problem (the blue circle).

This algorithm is called 'naive' because it checks a lot of lattice points that are not best approximations. Cassels' algorithm improves this substantially for the function $F(p,q)$. It's likely that a similar improvement is possible for $A(x)$ and someone might wish to work it out.

The OP (particularly on stack overflow, but also here) seems to have thrown quite a number of red herrings into the problem statement (rather annoying). For example, knowledge of the point (x,y) does nothing other than tell you which quadrant you are looking in. The statement about it converting the problem to NP-complete rather than NP-hard doesn't make any sense. Also, the fact that the lines pass through integer points appears to be irrelevant.

Answer 2

I don't believe that the problem is NP-complete, because you're working in a fixed dimension. Since most of us believe that the hardest case is when the angle between the two lines is very small, then you can take the line L half way in between the two lines, and orthogonally project onto it for the objective function. Then we have an integer programming problem in 2 dimensions -- the two inequalities specify the proper side of the two lines. Hendrik Lenstra in "Integer Programming with a fixed number of variables" in Mathematics of Operations Research, showed that when the dimension is fixed there is a polynomial time algorithm for IP (using a variant of the L^3 lattice basis reduction). There's also the paper http://www.math.uni-klu.ac.at/or/doctoralschool/deloera.pdf "Integer Polynomial Optimization in Fixed Dimension" which mentions that a convex polynomial objective function also has a polynomial time algorithm in fixed dimension, so that should do it for this problem.

[Added Comments] Using the two papers that I mentioned in my comments http://mpi-inf.mpg.de/~soeren/pubs/2ip.ps http://homepages.cwi.nl/~aardal/journal_rev.ps one can proceed as follows: We're going to have upper and lower supporting lines orthogonal to the midpoint line L. These will always have the property that we know that there is at least one lattice point in the quadrilateral bounded by the the original lines and the supporting line (though at the beginning, it has degenerated into a triangle). At the beginning the lower supporting line is the one passing through the intersection of the two original lines. The upper supporting line can be calculated by noticing that any disk of radius > sqrt(2)/2 must contain a lattice point, so you can find the smallest distance along the midline where you can place the center of such a disk -- it will be where the line from center orthogonal to each of the original lines has distance sqrt(2)/2. By simple trigonometry you can see that if the number of bits in the original number is N, then we need at most 2N bits to specify this point (i.e. if the denominators are around n for the originals, then the denominators for the above points are at most n^2). Now use the algorithm in first paper which tells you in time linear in the number of bits of the problem whether or not there are any lattice points in the quadrilateral. Do a binary search by looking at the a test line half way in between the supporting lines, and testing each of the two quadrilaterals. After only about N steps (remember that N is the log of the coefficients) you'll be down to a quadrilateral with a small area and width. At that point you can quickly enumerate all lattice points in it and test them for the minimum. This algorithm probably runs in time O(N^2) where N is the number of bits in the original coefficients.

[another addition to simplify things]:

The idea in solving the problem is to do the following:

1) Find a good enough approximation to the minimum distance, and the lattice point attaining that, so that you can enclose that region in a rectangle whose sides are parallel to the $x$ and $y$ axes, with a small enough area, $A$ and perimeter $P$. It's easy to see that you can enumerate all lattice points inside of such a rectangle in time $A+P$, and then check those to see which ones give you the minimum. [Changing notation] Let $n$ be the number of bits to specify the problem and $N=2^n$. Let the angle between the two lines be denoted by $2 \theta$, and the midline by $L$.

2) As I mentioned, any disk with radius $> \sqrt{2}/2$ must contain a lattice point. We use this by placing such a disk in between the two lines, and as close to the point of intersection as possible. We see that the distance from the intersection point to the center of the disk which just fits, is $\sqrt{2}/2 \csc \theta$. If $\theta \ge \pi/6$ (say) we can enclose the whole triangle bounded by the lines and the tangent line to the disk orthogonal to $L$ in a rectangle of constant size and perimeter, so we just try all of those points. Note that in any case $\theta \ge c/N$ for some absolute constant $c$ since when $\theta$ is small enough $\cos \theta \approx 1 - \theta^2/2$, and $\cos 2 \theta$ is given by a dot product between the coefficients of the lines and so has at most $2n$ bits.

3) Otherwise we call the algorithm described in http://mpi-inf.mpg.de/~soeren/pubs/2ip.ps just once, with the objective function the dot product between $(x,y)$ and a vector parallel to $L$, and the four constraints: between the two lines, and the dot product with $L$ is $\ge 0$ and $\le$ the bound we get by placing the disk. Because $\theta$ is small enough and bounded away from 0, we can now draw a rectangle enclosing the point produced by the algorithm whose area and perimeter are bounded, independent of the number bits, which must contain the answer, and we again enumerate all points in that rectangle.

Answer 3

You can try to apply the algorithm from the paper

B. N. Delone, “An algorithm for the “divided cells” of a lattice”, Izv. Akad. Nauk SSSR Ser. Mat., 11:6 (1947), 505–538

Answer 4

Let $n$ be the input size, i.e. the number of input digits in some base, or the sum or max of the logs of the magnitudes of the input numbers-- these are all equivalent for the purposes of big-O.

Strategy: apply at most $O(n)$ shears (area preserving linear transformations leaving one coord axis fixed) to get an easy problem, solve the easy problem, and then apply the inverse sequence of shears to get the solution back in the original coordinate space.

EASY PROBLEM #1: when at least one of the two bounding rays is axis-aligned.

EASY PROBLEM #2: when the two bounding ray directions point into the interiors of different quadrants of the plane.

Details of the easy problems are omitted, since they are straightforward and not the interesting or challenging part of this problem.

So assume we have a non-easy case of the problem: that is, both bounding ray directions point into the interior of the same quadrant; w.l.o.g. both point into the first quadrant.

Choose one of the two bounding rays and call it "primary" and the other "secondary". If the primary ray has slope $< 1$, swap the coordinate axes, so that the primary ray has slope $\ge 1$.

Let $a/b \ge 1$ be the slope of the primary bounding ray, with $a \ge b \gt 0$ taken directly from the input equations. Let $q,r$ be the quotient and remainder of dividing $a$ by $b$. That is, $q=\mathrm{floor}(a/b), r=a-q b, a \ge b \gt r \ge 0$.

Apply to the problem geometry the shear that leaves the $y$ axis fixed and takes $(1,q)$ to $(1,0)$; in other words, the linear transformation that leaves $(0,1)$ fixed and takes $(1,0)$ to $(1,-q)$. Maintain the designation "primary" on the image of the primary bounding ray. This shear transformation decreases the primary ray's slope to r/b, either keeping its direction into the first quadrant (if $r>0$), or parallel to the $x$ axis (if $r=0$).

Note that the solution before this shear transformation can be described by any of the following equivalent characterizations (equivalent because both bounding ray dirs point into the first quadrant):

• (a) "the integer point in the quarter closest to the intersection point" (by definition)
• (b) "the integer point in the quarter with minimal $x+y$",
• (c) "the integer point in the quarter with minimal $x$ among all points with minimal $y$".
• (d) "the integer point in the quarter with minimal $y$ among all points with minimal $x$"

There are three cases.

Case 1: After the shear, both rays still point into the interior of the first quadrant.

In this case, the solution to the sheared problem is, again, the following equivalent characterizations:

• (a') "the integer point in the transformed quarter closest to the transformed intersection point"
• (b') "the integer point in the transformed quarter with minimal $x+y$"
• (c') "the integer point in the transformed quarter with minimal $x$ among all points with minimal $y$".
• (d') "the integer point in the transformed quarter with minimal $y$ among all points with minimal $x$"

By focusing on (d) and (d') it's clear that the solution to the transformed problem is the transformed solution to the original problem.

In the transformed problem, the primary bounding ray has slope $r/b$ with $b > r > 0$ (still in Case 1 here) so its slope is $< 1$; therefore as we start analysis of the transformed problem as before, we'll swap the coords so that its slope $b/r > 1$. Notice that what we've done so far is perform one iteration of the euclidean algorithm; that is, $a$ has been replaced by $b$, and $b$ has been replaced by the remainder of $a$ divided by $b$.

Repeat as often as we find ourselves in Case 1 (keeping the "primary" designation on the same one of the two rays even as they get transformed). Since the $a,b$ values are following the euclidean algorithm starting with two of the original input numbers, well-known analysis of the euclidean algorithm tells us Case 1 can happen at most $O(n)$ times: that is, if we get all the way to the end of the euclidean algorithm, $b$ will become 0 which will take us out of Case 1 (if we haven't already left it before that).

Case 2: The shear makes the primary bounding ray horizontal, and the sheared secondary ray still points into the first quadrant.

Exactly as in Case 1, (a')=(b')=(c')=(d') and so the solution to the transformed problem is the transform of the solution to the original problem. Furthermore the transformed problem is now a case of EASY PROBLEM #1; solve that, inverse-transform the solution back to the original space; done. (When this is encountered recursively from Case 1, it means we've made it all the way to the end of the euclidean algorithm.)

Case 3: The shear takes the secondary bounding ray direction out of the interior of the first quadrant, to either the $x$ axis or the fourth quadrant. (This is regardless of whether the primary bounding ray direction became horizontal or stayed in the first quadrant).

This case is relatively easy to solve, but we have to be careful-- in this case the transform of the point inside the quarter closest to the intersection (i.e. satisfying (a)) is not necessarily the point in the transformed quarter closest to the transformed intersection (i.e. satisfying (a')). However, it is true that the transform of the point satisfying (d) is the point satisfying (d'). Therefore the transformed problem, though easy, is not in the form of the original problem. Instead, it's a case of:

EASY PROBLEM #3: The bounding ray directions both point into the +x half plane, but do not both point into the interior of the same quadrant, and we are asked to find the point satisfying (d) rather than (a) (which are not necessarily the same in this case).

So I've described an algorithm for transforming the original problem into one of three easy problems, in $O(n)$ iterations.

Complexity analysis

The number of iterations is $O(n)$, so the overall runtime is $O(n)$ times the cost per iteration.

So how expensive is an iteration? Each iteration consists of a small constant number of arithmetic operations, each of which is at most $O(M(\mathrm{num\ digits\ in\ operands}))$ where $M$ stands for the cost of one multiplication (see the wikipedia article "Computational cost of mathematical operations"). But how many digits is that?

Well, at first glance it appears that intermediate and final values could have up to $n^2$ digits which would make the overall runtime $O(n M(n^2)))$. Furthermore one might suspect the $n$-digit witness point could help decrease this bound. But the following closer analysis reveals it's not as bad as that, and the $n$-digit witness actually makes no difference at all.

Each transform is a shear followed by a coord axis swap: $$ \begin{pmatrix} 0 & 1 \\ 1 & 0 \\ \end{pmatrix} \begin{pmatrix} 1 & 0 \\ - q_i & 1 \\ \end{pmatrix} = \begin{pmatrix} - q_i & 1 \\ 1 & 0 \\ \end{pmatrix} $$ where $q_0, q_1, ...$ is the sequence of quotients $q$ occurring in the euclidean algorithm (possibly cut short by Case 3). So the matrix norm of the cumulative matrix is $\leq$ the product of the $q_i$'s, which is at most the original value of $a$, which has $n$ digits. And the same can be said about the cumulative inverse matrix. That tells us the maximum number of digits occuring in the answer or any intermediate value is $O(n)$. Therefore the overall runtime is $O(n M(n))$.

Incidentally, the sequence of 2x2 integer matrix computations being followed is actually well known as the EEA (Extended Euclidean Algorithm) for CRT (Chinese Remainder Theorem). If taken all the way to completion, the cumulative matrix represents an invertible linear transformation with integer coeffs that takes the original primary ray direction to a coordinate axis, which can be a useful building block for many applications (notably Project Euler problems :-) ). So one might think we could simply use a prepackaged EEA implementation as a building block; but I don't see how to make that work for this problem, due to the possibility of needing to stop early due to Case 3.

Answer 5

Doesn't the following work? UPDATE: Sorry, no it doesn't. I missed the requirement that the lattice point be on the correct side of the two lines.

Find the exact coordinates, $(u,v)$ of the intersection. These are rational numbers, and can be found by solving two linear equations. Then the nearest lattice point is $(\mathrm{ROUND}(u), \mathrm{ROUND}(v))$ where $\mathrm{ROUND}$ rounds to the nearest integer.

To see that this is the closest point, translate $(\mathrm{ROUND}(u), \mathrm{ROUND}(v))$ to the origin. So we need to show that, if $|u|$ and $|v|$ are $< 1/2$, then $(u,v)$ is closer to $(0,0)$ than to any other lattice point $(a,b)$. The cases of $(a,b) = (0, \pm 1)$, $(\pm 1, 0)$ or $(\pm 1, \pm 1)$ can be checked by hand. For any other $a$, $b$, one of $|u-a|$ and $|v-b|$ is greater than $1$, so the distance from $(u,v)$ to $(a,b)$ is at least $1$; which is much greater than $\sqrt{u^2+v^2} \leq \sqrt{2}/2$.

For a formal complexity analysis, the size of your input is the log of the number of digits needed to specify the two lines. All the arithmetic operations need to solve linear equations can be done in time polynomial in this logarithm.

This question is much more interesting for lattices other than the square grid. See Voronoi decomposition for details.

Answer 6

This problem is a little bit obscure, but there are several efficient solutions which you may be interested in. I'd start by first trying to find a copy of Lee and Chang's paper, unless you're fluent in French, in which case Guigue's thesis contains an updated exposition.

Mehta, S., M. Mukherjee, and G. Nagy. "Constrained integer approximation to 2-D line intersections." Second Canadian Conf. on Computational Geometry. 1990.

Lee, H. S., and R. C. Chang. "Approximating vertices of a convex polygon with grid points in the polygon." Algorithms and Computation. Springer Berlin Heidelberg, 1992. 269-278.

Guigue, Philippe. Constructions géométriques à précision fixée. Diss. Université de Nice Sophia-Antipolis, 2003.