Let $n$ be the input size, i.e. the log of the max magnitude of any of the input numbers.

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-qb, 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.

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$ 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.  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 leaves the primary bounding ray direction in the interior of the first quadrant,
but it takes the secondary bounding ray dir out of the interior of the first quadrant: that is, to either the $x$ axis or the fourth 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](http://en.wikipedia.org/wiki/Computational_complexity_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.