# Avoiding overflow in finding a solution to $ax + by = c$

Can anyone come up with a nice way of computing a solution to the linear diophantine equation $ax + by = c$ where $a,b,c \in \mathbb{Z}$ and $\gcd(a,b) \mid c$, such that all the calculations are carried out without any of the intermediate results exceeding $\max\{|a|, |b|, |c|\}$ in absolute value?

I.e. if $a,b,c$ are 64-bit integers, all the calculations should be done using 64-bit integers only.

The standard method is of course to use Euclidean algorithm to first find a solution to $ax + by = \gcd(a,b)$ and then multiplying $x$ and $y$ by $\frac{c}{\gcd(a,b)}$. The Euclidean algorithm is fine, but this last multiplication might go out of range.

• If the last multiplication goes out of range, then the final answer is out of range, which means you're asking the impossible. What am I missing? – Jeffε Oct 30 '13 at 22:36
• This might be more suitable for Computational Science as it seems like a numerical analysis type question. – Kaveh Oct 30 '13 at 23:29
• @JeffE: I'm looking for a single solution. From the solution you get from the final multiplication you can construct smaller ones by adding multiples of $\frac{b}{\gcd(a,b)}$ to $x$ and $\frac{-a}{\gcd(a,b)}$ to $y$. If I'm not mistaken, there will always be a solution which is smaller than $\max(|a|,|b|,|c|)$. – J. J. Oct 31 '13 at 8:15
• @Kaveh: Thanks. I think there might be some interest theoretical interest in this as well. Specifically: Is it possible to somehow modify the extended Euclidean algorithm so that it stays stable for the general equation $ax + by = c$ like it does for Bezout's identity. I realize that this problem is not (usually) really a problem in practice since one can of course use multiprecision integers for the last step. – J. J. Oct 31 '13 at 8:42
• Do you mean $\gcd(a,b) \mid c$, instead of $c \mid \gcd(a,b)$? Also, WLOG we can assume $\gcd(a,b)=1$ (if not, divide $a,b,c$ by $\gcd(a,b)$). So, the problem is, given $a,b,c$ where $\gcd(a,b)=1$, compute $x,y$ such that $ax+by=c$, without overflow in any intermediate computation. – D.W. Oct 31 '13 at 23:51

One way to solve this might be through an application of the Chinese remainder theorem and some kind of Hensel lifting.

Let's fix a constant $$\alpha$$ with the property that $$2^\alpha 3^\alpha 5^\alpha > 2^{129}$$ and $$5^\alpha < 2^{64}$$. For instance, $$\alpha=27$$ would work fine.

Suppose $$ax+by=c$$ holds modulo $$2^\alpha$$, modulo $$3^\alpha$$, and modulo $$5^\alpha$$, i.e., $$ax+by \equiv c \pmod{q}$$ for $$q=2^\alpha$$, $$q=3^\alpha$$, and $$q=5^\alpha$$. Then by the Chinese remainder theorem, it follows that $$ax+by=c$$ holds in the integers. So, our goal will be to find $$x,y$$ that satisfy this equation modulo $$2^\alpha$$, modulo $$3^\alpha$$, and modulo $$5^\alpha$$.

Here are some observations that will help us find a solution that works modulo each of $$2^\alpha$$, $$3^\alpha$$, and $$5^\alpha$$:

• Given a solution that is valid modulo $$2^i$$, $$3^i$$, and $$5^i$$, it is easy to find a solution that is valid modulo $$2^{i+1}$$, $$3^{i+1}$$, and $$5^{i+1}$$. Here is how. Given $$(x,y)$$ such that $$ax+by \equiv c \pmod{q}$$ (for $$q\in \{2^i,3^i,5^i\}$$), we consider the 3481 pairs $$(x',y')$$ such that $$x' \in \{x,x+ d \cdot 2^i 3^i 5^i : |d| < 30, d \in \mathbb{Z}\}$$ and $$y' \in \{y,y+ d \cdot 2^i 3^i 5^i : |d| < 30, d \in \mathbb{Z}\}$$.

One of these pairs will satisfy $$ax+by \equiv c \pmod{q}$$ (for $$q\in \{2^{i+1}, 3^{i+1}, 5^{i+1}\}$$) and will avoid overflow in $$x',y'$$ (i.e., we will have $$0 \le x',y' < 2^{64}$$). So, you can check all 3481 candidates and find the one that you wanted. Identifying the one you wanted is easy. Given a candidate, it is easy to check for overflow and immediately rule out any $$x',y'$$ that trigger overflow or underflow in the expressions $$x+ d \cdot 2^i 3^i 5^i,y+ d \cdot 2^i 3^ i5^i$$. Also, you can check whether $$ax+by=c$$ holds modulo $$2^{i+1}$$, $$3^{i+1}$$, and $$5^{i+1}$$ separately in 64-bit arithmetic without overflow, as long as each of $$2^{i+1},3^{i+1},5^{i+1}$$ is at most $$2^{64}$$.

• It is also easy to find a solution that is valid modulo $$2^0 3^0 5^0$$: simply take $$x=0,y=0$$.

Given all of this background, we can now use induction to find a solution. We start by forming a solution that works for $$i=0$$. At each step, we increment $$i$$, until we reach $$i=\alpha$$. Each step can be done in 64-bit arithmetic semi-efficiently (with a few thousand basic arithmetical operations), so the whole computation should be semi-efficient.

I feel like there ought to be some way to work modulo $$2^i$$, $$3^j$$, and $$5^k$$, and at each step separately increment either $$i$$, $$j$$, or $$k$$, but I haven't worked out the exact details of how to do that yet. If you can make that work, the resulting scheme would be significantly more efficient.

• An interesting way of approach. However I'm not seeing why in the step where you are finding solutions for $i+1$ one of those pairs will work. Also is there a reason we are working in $2^i 3^j 5^k$ and not just $2^i$? – J. J. Nov 6 '13 at 19:05
• @J.J., that's the Hensel lifting part. Are you familiar with why Hensel lifting works? It can be proven in the same way. If you're not familiar with it, it's a nice exercise to prove it (it's by induction on $i$, where the inductive hypothesis is that the partial solution you've currently got mod $2^i 3^i 5^i$ is congruent, modulo $2^i 3^i 5^i$, to the desired solution in the integers). That's enough to show that there exists a way to extend the current partial solution mod $2^i 3^i 5^i$ to one mod $2^{i+1} 3^{i+1} 5^{i+1}$. – D.W. Nov 6 '13 at 19:19
• As far as why I work mod $2^i 3^i 5^i$ instead of mod $2^i$: it's because you said you want to work in 64-bit arithmetic only. What we can say is that if $ax+by \equiv c \pmod q$ holds (modulo $q$), and if $q > 2^{129}$, then it follows that $ax+by=c$ (in the integers). If you wanted to just work mod $2^i$, you'd need to work mod $2^{129}$, which can't be done in 64-bit arithmetic. Where does $2^{129}$ come from? It's a consequence of the fact that $|a|,|x|,|b|,|y|,|c| < 2^{64}$; consequently, for any solution in the integers (where $ax+by=c$), we must have $|ax+by| < 2^{129}$. – D.W. Nov 6 '13 at 19:22
• why cant he just work with $2^\alpha 3^\alpha >2^{129}$? $\alpha$ can be $50$ now? Or just $5^\alpha > 2^{129}$? Here $\alpha$ would be $56$? – T.... Nov 7 '13 at 5:04
• @JAS, I already answered that question. Again, $5^{55}$ is greater than $2^{64}$, so you cannot do computations modulo $5^{56}$ within 64-bit arithmetic. For example, to go from a solution that's valid modulo $5^{54}$ to a solution that's valid modulo $5^{55}$ requires checking which candidate solution is valid modulo $5^{55}$, which requires computing with numbers that are potentially as large as $5^{55}-1$, and those numbers cannot be expressed within 64-bit arithmetic. – D.W. Nov 7 '13 at 7:16

Here is a better solution that should be extremely efficient (within a $2\times$ factor as fast as the extended Euclidean algorithm) and never overflows 64-bit arithmetic. Without loss of generality, let's assume that $a>b>0$ and $\gcd(a,b)=1$.

Definition. Suppose $0\le c\le b$ and $x,y \in \mathbb{N}$. We'll say that $(x,y)$ is a representation of $c$ if (i) $ax-by=c$ and (ii) $0 \le x < b$ and $0 \le y < a$.

The basic idea. Our goal is to find a representation of $c$, using 64-bit arithmetic. Here's the overall approach I suggest. First, I'll show below how to obtain a representation of $1$. Next, I'll show that if we have a representation of $c_1$ and a representation of $c_2$, we can obtain a representation of $c_1+c_2$ (subject to a few minor conditions). This will make it easy to obtain a representation of $c$: we form an addition chain that ends with $c$, and then use the preceding observations to calculate a representation of $c$. Details follow below.

Fact 1. If $0 \le c \le b$, then there exists a representation of $c$, and this representation is unique.

Proof sketch. To show existence, we can take $y$ be $b^{-1} c \bmod a$ and $x = (by+c)/a$. To show uniqueness, if there are two representations of $c$, subtract them; we obtain $x,y$ such that $ax-by=0$, $0\le x < b$, $|y| < a$, but this is only possible when $x=y=0$.

Fact 2. We can find a representation of $1$, within 64-bit arithmetic.

Proof sketch. Use the extended Euclidean algorithm. None of the intermediate values exceed $\max(a,b,c)$.

For what comes next, define an operator $\oplus$ on representation as follows:

$$(x_1,y_1) \oplus (x_2,y_2) = \begin{cases} (x_1+x_2,y_1+y_2) &\text{if y_1+y_2<a}\\ (x_1+x_2-b,y_1+y_2-a) &\text{otherwise.} \end{cases}$$

Fact 3. If $(x_1,x_2)$ is a representation of $c_1$ and $(x_2,y_2)$ is a representation of $c_2$ and $0 \le c_1,c_2,c_1+c_2 \le b$, then $(x_1,y_1) \oplus (x_2,y_2)$ is a representation of $c_1+c_2$.

Proof sketch. We can extend the proof of uniqueness above to demonstrate that if $(x_3,y_3) = (x_1,y_1) \oplus (x_2,y_2)$ and if $0 \le y_3 < a$, then $0 \le x_3 < b$. Now the fact that $(x_3,y_3)$ is a representation of $c_1+c_2$ is immediate, by linearity.

Notice that the definition of the operator $\oplus$ ensures that you never overflow 64-bit arithmetic.

So, the algorithm becomes straightforward. Choose an addition chain that ends in $c$. In other words, choose a sequence $c_1,c_2,\dots,c_m$ such that $c_1=1$ and $c_m=c$ and for all $k$, there exists $i,j$ such that $c_i+c_j=c_k$ (where $1 \le i,j < k$). There are standard ways to choose such a sequence of length $\le 2 \lg c$.

Next, use the extended Euclidean algorithm to find a representation of $c_1=1$. Finally, iteratively sweep forward: for $k=1,2,\dots,m$, find a representation of $c_k$ (using the representations of $c_1,\dots,c_{k-1}$ and $\oplus$, as suggested by the addition chain). We end with a representation of $c_m=c$, as desired. No step of this process overflows 64-bit arithmetic, so this satisfies all your desired requirements.

This property is quite efficient. The extended Euclidean algorithm requires at most $2 \lg c$ iterations, and the addition chain is of length at most $2 \lg c$, so you do at most $4 \lg c$ steps, where you do a handful of simple 64-bit arithmetic operations in each step. That's about as efficient as you could hope for.

• I appreaciate your nice work. Do you think this could somehow be extended to the case $c > b$? – J. J. Nov 11 '13 at 17:17
• @J.J., oh, good question! I don't know. Thank you for calling attention to that gap in my answer -- I hope someone will have some ideas that might help. – D.W. Nov 11 '13 at 19:04

Here is a modified version of the Euclidean algorithm that is stable.

When considering the equation $$ax + by = c \quad (1)$$ we may w.l.o.g. assume that $a, b, c \ge 0$, $b \ge 1$ and $b \ge a$. From now on we will only consider such equations. Moreover, we will set $d = \gcd(a,b)$.

We call a solution $(x,y)$ good if $0 \le x < \frac{b}{d}$. Such a solution always exists, since all solutions are parametrized by $(x + t \cdot \frac{b}{d}, y - t \cdot \frac{a}{d})$.

If $(x,y)$ is a good solution, then $$\frac{c}{b} - \frac{a}{d} < y \le \frac{c}{b} \quad (2)$$ This follows simply from the bounds on $x$ and the fact that $y = \frac{c}{b} - \frac{ax}{b}$. Notice that in particular $|y| \le \max(a,c)$.

We will now solve the problem by induction. Assume first that $a=0$. Then $x=0$, $y = \frac{c}{b}$ is a good solution. Let us then suppose by induction that we can obtain good solutions for all coefficients smaller than $a,b,c$. We let $(x',y')$ be a good solution to the equation $$(b - \left\lfloor \frac{b}{a} \right\rfloor a) x' + a y' = c \quad (3)$$ Then $(\tilde{x}, \tilde{y}) = (y' - \left\lfloor \frac{b}{a} \right\rfloor x', x')$ is a solution to (1). First of all we can calculate $\left\lfloor \frac{b}{a} \right\rfloor x'$ without overflow since $0 \le x' < \frac{a}{d}$. We note that (2) applied to the equation (3) gives us $$\frac{c}{a} - \frac{b - \left\lfloor \frac{b}{a}\right\rfloor a}{d} < y' \le \frac{c}{a}.$$ Hence it follows that $$\frac{c}{a} - \frac{b}{d} < y' - \left\lfloor \frac{b}{a} \right\rfloor x' \le \frac{c}{a} \quad (4)$$ Thus $|\tilde{x}| \le \max(b,c)$ and $\tilde{x}$ can be computed without overflow.

Finally it remains to normalize $(\tilde{x}, \tilde{y})$ to a good solution. This involves computing $x = \tilde{x} + t \cdot \frac{b}{d}$ and $y = \tilde{y} - t \cdot \frac{a}{d}$ where $t = \left\lceil \frac{-\tilde{x}}{b/d} \right\rceil$. Assume first that $\tilde{x} < 0$. From (4) we see that then $t=1$ and everything can be done without overflow. Assume then that $\tilde{x} \ge 0$. Then $|t| \le \frac{\tilde{x}}{b/d}$. Because $b \ge a$, both products $t \cdot \frac{b}{d}$ and $t \cdot \frac{a}{d}$ can be calculated without overflow. We are done.