Efficient modular solution to $ax - by \equiv 0\pmod{p}$? Given prime $p$, integers $x$ and $y$ where both $x, y < p$ and $x \neq y$, is there an efficient way to find nontrivial coefficients $a, b$ where $a, b < \sqrt p$ such that 
$$ax - by \equiv 0\bmod p$$
Further, assume that we are told that such a pair $a, b$ exists; the question is, what's the best way to find them?
If there is no efficient (non-brute force) method, then assume we have many such pairs $a_ix_i - a_jy_j \equiv 0\bmod p$ where the $a_i, a_j$ are known to exist for their respective $x_i, y_j$ pairs, and are bounded by $\sqrt p$ as above; is there an efficient way to find at least one satisfying pair $a_i, a_j$?
 A: I'm somewhat new to this board, but if Python is permitted, here is an annotated function that implements Jaap's efficient algorithm (with some helpers up front):
(Again, sorry if posting code is inappropriate for this board - please let me know, and I will delete it.  If appropriate, feel free to delete this comment! ;-)
IMPLEMENTATION
def xgcd(b, n):
    '''
    Returns GCD of b and n - specifically, returns the GCD(b,n), 
    and Bezout's coefficients x and y.  
    From https://en.wikibooks.org/wiki/Algorithm_Implementation/Mathematics/Extended_Euclidean_algorithm#Iterative_algorithm_3
    '''
    x0, x1, y0, y1 = 1, 0, 0, 1
    while n != 0:
        q, b, n = b // n, n, b % n
        x0, x1 = x1, x0 - q * x1
        y0, y1 = y1, y0 - q * y1
    return  b, x0, y0

def mulinv(b,n):
    '''
    Returns multiplicative inverse of b mod n, or None if no inverse exists.
    '''
    g, x, _ = xgcd(b, n)
    if g == 1:
        return x % n
    else:
        return None

def jaapfastsolve(x, y, p):
    '''
    Returns solution of a, b, (with number of steps) to the form: ax - by = 0 mod p
    '''
    isqrt = lambda x: int(pow(x,0.5)) # convenience def
    mins = (1, x*mulinv(y,p)%p)
    maxs = mins
    n = 0
    while mins[1]>isqrt(p):
        n = n + 1
        # increase max (by max) to largest under p
        tm = p//maxs[1]
        news = (tm*maxs[0], tm*maxs[1]%p)
        maxs = max([maxs, news], key=lambda t:t[1])
        # increase again by min to largest under p to get new max
        tn = (p-maxs[1])//mins[1]
        news = ((tn*mins[0])+maxs[0], (tn*mins[1] + maxs[1])%p)
        maxs = max([maxs, news], key=lambda t:t[1])
        # add one more min to overflow p to get new min
        news = (mins[0]+maxs[0], (mins[1] + maxs[1])%p)
        mins = min([mins, news], key=lambda t:t[1])
        # print mins
    return (n, mins)

EXAMPLE
p=1000003; x=454463; y=109818
n, mins = jaapfastsolve(x, y, p)
print "fast solve answers"
print n, mins, (mins[0]*x-mins[1]*y)%p==0

Returns 4 (61, 43) True.
