Efficient modular solution to $ax - by \equiv 0\pmod{p}$?

Question

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$?

Answer 1

It is not always possible for both of $a,b$ to satisfy $a,b<\sqrt{p}$.

A small counterexample is $p=5$, $x=1$, $y=4$. The only solutions are $(a,b) \in \{ (1,4), (2,3), (3,2), (4,1) \}$.

In fact, for every $p$ the values $(x,y)=(1,p-1)$ will lead to a counterexample because it forces $a+b=0 \mod p$, and $a+b < \sqrt{p}+\sqrt{p}<p$ then gives the trivial solution $a+b=0$ only.

For any $a$, you can set $b=y^{-1}ax \mod p$, where the $y^{-1}$ is calculated with the Extended Euclidean algorithm applied to the coprime numbers $y$ and $p$. You now have a solution to $ax-by=0 \mod p$. You just cannot guarantee that if $a<\sqrt{p}$ that you will get a $b$ that satisfies $b<\sqrt{p}$.

EDIT:

To find a solution that satisfies $a,b<\sqrt{p}$ (assuming there is one) you could just start at $a=1$ and try successive values.

For example $p=101$, $x=32$, $y=6$. Starting with $a=1$, you get the solution: $$b=6^{−1}\cdot32=17\cdot32=544=39 \mod 101$$ This $b$-value isn't in the correct range. To get the other solutions for $a=2,3,4,...$ you take successive multiples of $b=39$, which are $39,78,16,55,94,32,71,9,...$. You stop then because you found $b\le10$, which gives the valid solution $(a,b)=(8,9)$.

For large $p$, trying all values in succession may be slow or even infeasible. You can however be a little cleverer.

To illustrate this, here's a larger example: $p=1000003$, $x=454463$, $y=109818$.

The first solution has $a=1$ and $b=454463*109818^{-1} = 454463*554765=409838$.

$(a,b)=(1, 409838)$

Get the first multiple of this that doesn't yet overflow, i.e. multiply it by $\lfloor{p/b}\rfloor$.
$2\cdot(1, 409838) =(2, 2\cdot409838) = (2,819676)$
In this method we keep track of the solution with the smallest value of $b$ and of the solution with the largest value of $b$. So at the moment the $a=2$ solution is the current largest solution, and the $a=1$ solution is the current smallest solution.

Now add the smallest solution to the largest solution so that $b$ overflows modulo $p$.
$(2,819676)+(1, 409838) = (3,229511)$
This gives a new and improved current smallest solution.

At this point we would normally try to improve our largest solution by adding the smallest to it. However, it would overflow so that wouldn't give a larger $b$.
Let's add it anyway to get a new smallest solution:
$(2,819676)+(3,229511) = (5,49184)$

Again, let's see what we get when we add the smallest solution to the largest. This time we can it 3 times to get a new largest solution:
$(2,819676)+3\cdot(5,49184) = (17,967228)$
And then add it one more time to make it overflow and produce a new smallest solution:
$(17,967228)+(5,49184) = (22,16409)$

Again add the smallest solution as often as possible to the largest without making it overflow:
$(17,967228)+1\cdot(22,16409) = (39,983637)$
And then add it once more to make it overflow:
$(39,983637)+(22,16409) = (61,43)$

The algorithm ends now, because at this point $b$ has become smaller than $\sqrt{p}$ as required. You would also stop if $a$ became greater than $\sqrt{p}$, in which case there is no solution satisfying the condition $a,b < \sqrt{p}$.

Answer 2

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.