Given two lines on a plane, how to find integer points closest to their intersection?

This problem falls into the category of Integer Convex Optimization.

Presented here is a mathematical way to approach the problem. I don't expect you to actually use it - a lot of complicated techniques are required, and other algorithmic approaches (such as "searching" for the appropriate point) will likely do just fine. However, interest has been expressed in the "true" solution, so here it is.

It can be solved in three stages:

1. First, determine which side of each line the answer will be on, as illustrated by TheMachineCharmer's answer.
2. Once that is known, the problem can be rewritten as a convex optimization problem (see Wikipedia for details). The function to be optimized is minimizing (x - x0)^2 + (y - y0)^2, with x0 and y0 the coordinates of the intersection of the two lines. The two lines each become a linear inequality, e.g. "x+y >= 0", together forming the convex region the answer can be found in. I will note that the solution will be (x=x0, y=y0) - what you need from this stage a way of expressing the problem, analagous to a feasible tableau for the simplex method.
3. Third, an integer solution can be found by repeatedly adding cuts to further constrain the feasible region until the solution to the convex optimization problem is integral. This stage may take a lot of iterations in the general case, but looking at the problem presented, and in particular the 2D nature of it, I believe it will be solved with at most two cuts.

I show here how a "difficult" instance of this problem can be solved. I think this method can be generalized. I have put another simpler instance in the comments of the original post.

Consider the two lines:

10000019 * X - 10000015 * Y + 909093 >= 0    (L1)
-10000022 * X + 10000018 * Y + 1428574 >= 0  (L2)
A = 10000019, B = -10000015, C = -909093

The intersection point is H:

Hx = -5844176948071/3, Hy = -5844179285738/3

For a point M(X,Y), the squared distance HM^2 is:

HM^2 = (9*X^2+35065061688426*X
    +68308835724213587680825685
    +9*Y^2+35065075714428*Y)/9

g = gcd(A,B) = 1: the equation of L1 A*X+B*Y+909093 can take any integer value.

Bezout coefficients, U=2500004 and V=2500005 verify:

A * U + B * V = 1

We now rewrite the problem in the "dual" basis (K,T) defined by:

X = T*U - K*B
Y = T*V + K*A

After substitution, we get:

T+909093 >= 0
2*T+12*K+1428574 >= 0
minimize 112500405000369*T^2
   +900003150002790*T*K
   +1800006120005274*K^2
   +175325659092760325844*T
   +701302566240903900522*K
   +Constant

After further translating (first on T, then on K to minimize the constant in the second equation), T=T1-909093, K=K1+32468:

T1 >= 0
2*T1+4+12*K1 >= 0
minimize 112500405000369*T1^2
    +300001050000930*T1
    +900003150002790*T1*K1
    +1200004080003516*K1
    +1800006120005274*K1^2
    +Constant

The algorithm I proposed is to loop on T1. Actually, we don't need to loop here, since the best result is given by T1=K1=0, corresponding to:

X = -1948055649352, Y = -1948056428573

My initial post below.

Here is another idea of algorithm. It may work, but I did not implement it...

With appropriate change of signs to match the posit

The problem of checking whether a point is part of a mathematical cone is fairly simple. Given 2 vectors, v, w, any point in the cone defined by (v, w) will be on the form: z = a***v** + b***w**, where a,b >= 0. Note, for this to work, you will have to move Origo to the intersection of the 2 lines. Since we cannot assume finite precision of the intersection, you will have to do floating point math and decide whether something is close enough to what you want.

1. Find vectors defining the 4 cones (there's infinitely many of them, we just need 2 for each cone), that are defined by the 2 lines.
2. Find out which cone contains our point, call that cone for C.
3. Take the 2 vectors defining C, and find the median vector (the vector that would split C in 2 identical cones), call it m.
4. Now is time to initiate the loop. For simplicity sake I'm going to assume that we limit ourself to n-bits on the x and y axis. Note you'll need an integer larger than n-bits for the length of m. Now do a binary search along the length of m, and check the 2 rings around every time (I suspect 1 ring around will be enough). When you've found the smallest length that do contain points C, check which of those points are the closest.

The worst case growth would be O(log(sqrt(2*n^2)), where n is the length we use to represent the x and y axis.

It is possible to do a "reverse binary search" so to speak, if you don't know the length of *m. Just keep doubling the the length you go out until you find points in C. Then you know 2 points on m and you can do a binary search between them.

The main problem with all this is precision, so keep this in mind. Alternative ways