<p>Here is a modified version of the Euclidean algorithm that is stable.</p>

<p>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)$.</p>

<p>We call a solution $(x,y)$ <strong>good</strong> if $0 \le x &lt; \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})$.</p>

<p>If $(x,y)$ is a good solution, then
$$\frac{c}{b} - \frac{a}{d} &lt; 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)$.</p>

<p>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' &lt; \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} &lt; y' \le \frac{c}{a}.$$
Hence it follows that
$$\frac{c}{a} - \frac{b}{d} &lt; 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.</p>

<p>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} &lt; 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.</p>

CC BY-SA 3.0
2013-11-21T19:54:49.427Z
2013-11-21T19:54:49.427Z
J. J.
<p>One way to solve this might be through an application of the Chinese remainder theorem and some kind of Hensel lifting.</p>
<p>Let's fix a constant <span class="math-container">$\alpha$</span> with the property that <span class="math-container">$2^\alpha 3^\alpha 5^\alpha &gt; 2^{129}$</span> and <span class="math-container">$5^\alpha &lt; 2^{64}$</span>.  For instance, <span class="math-container">$\alpha=27$</span> would work fine.</p>
<p>Suppose <span class="math-container">$ax+by=c$</span> holds modulo <span class="math-container">$2^\alpha$</span>, modulo <span class="math-container">$3^\alpha$</span>, and modulo <span class="math-container">$5^\alpha$</span>, i.e., <span class="math-container">$ax+by \equiv c \pmod{q}$</span> for <span class="math-container">$q=2^\alpha$</span>, <span class="math-container">$q=3^\alpha$</span>, and <span class="math-container">$q=5^\alpha$</span>.  Then by the Chinese remainder theorem, it follows that <span class="math-container">$ax+by=c$</span> holds in the integers.  So, our goal will be to find <span class="math-container">$x,y$</span> that satisfy this equation modulo <span class="math-container">$2^\alpha$</span>, modulo <span class="math-container">$3^\alpha$</span>, and modulo <span class="math-container">$5^\alpha$</span>.</p>
<p>Here are some observations that will help us find a solution that works modulo each of <span class="math-container">$2^\alpha$</span>, <span class="math-container">$3^\alpha$</span>, and <span class="math-container">$5^\alpha$</span>:</p>
<ul>
<li><p>Given a solution that is valid modulo <span class="math-container">$2^i$</span>, <span class="math-container">$3^i$</span>, and <span class="math-container">$5^i$</span>, it is easy to find a solution that is valid modulo <span class="math-container">$2^{i+1}$</span>, <span class="math-container">$3^{i+1}$</span>, and <span class="math-container">$5^{i+1}$</span>.  Here is how.  Given <span class="math-container">$(x,y)$</span> such that <span class="math-container">$ax+by \equiv c \pmod{q}$</span> (for <span class="math-container">$q\in \{2^i,3^i,5^i\}$</span>), we consider the 3481 pairs <span class="math-container">$(x',y')$</span> such that <span class="math-container">$x' \in \{x,x+ d \cdot 2^i 3^i 5^i : |d| &lt; 30, d \in \mathbb{Z}\}$</span> and <span class="math-container">$y' \in \{y,y+ d \cdot 2^i 3^i 5^i : |d| &lt; 30, d \in \mathbb{Z}\}$</span>.</p>
<p>One of these pairs will satisfy <span class="math-container">$ax+by \equiv c \pmod{q}$</span> (for <span class="math-container">$q\in \{2^{i+1}, 3^{i+1}, 5^{i+1}\}$</span>) and will avoid overflow in <span class="math-container">$x',y'$</span> (i.e., we will have <span class="math-container">$0 \le x',y' &lt; 2^{64}$</span>).  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 <span class="math-container">$x',y'$</span> that trigger overflow or underflow in the expressions <span class="math-container">$x+ d \cdot 2^i 3^i 5^i,y+ d \cdot 2^i 3^ i5^i$</span>.  Also, you can check whether <span class="math-container">$ax+by=c$</span> holds modulo <span class="math-container">$2^{i+1}$</span>, <span class="math-container">$3^{i+1}$</span>, and <span class="math-container">$5^{i+1}$</span> separately in 64-bit arithmetic without overflow, as long as each of <span class="math-container">$2^{i+1},3^{i+1},5^{i+1}$</span> is at most <span class="math-container">$2^{64}$</span>.</p>
</li>
<li><p>It is also easy to find a solution that is valid modulo <span class="math-container">$2^0 3^0 5^0$</span>: simply take <span class="math-container">$x=0,y=0$</span>.</p>
</li>
</ul>
<p>Given all of this background, we can now use induction to find a solution.  We start by forming a solution that works for <span class="math-container">$i=0$</span>.  At each step, we increment <span class="math-container">$i$</span>, until we reach <span class="math-container">$i=\alpha$</span>.  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.</p>
<hr />
<p>I feel like there ought to be some way to work modulo <span class="math-container">$2^i$</span>, <span class="math-container">$3^j$</span>, and <span class="math-container">$5^k$</span>, and at each step separately increment either <span class="math-container">$i$</span>, <span class="math-container">$j$</span>, or <span class="math-container">$k$</span>, 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.</p>

CC BY-SA 3.0
2013-11-06T19:05:56.637Z
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.
CC BY-SA 3.0
2013-11-06T19:19:12.397Z
@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.
CC BY-SA 3.0
2013-11-06T19:22:12.953Z
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.
CC BY-SA 3.0
2013-11-07T05:04:43.183Z
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$?
Turbo
CC BY-SA 3.0
2013-11-07T05:36:15.360Z
@JAS, we can't work modulo $3^{50}$, because $3^{50} > 2^{64}$, and we're trying to work with just 64-bit arithmetic.  Same problem with for $5^{56}$.  My parameter choices were chosen so that $2^{\alpha}$, $3^{\alpha}$, and $5^{\alpha}$ are all smaller than $2^{64}$, so you can work modulo each of them using only 64-bit arithmetic.  I'll edit my answer to make this clearer -- I just noticed a typo in my answer that made this unnecessarily opaque.
D.W.
CC BY-SA 3.0
2013-11-07T07:16:57.100Z
@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.
CC BY-SA 3.0
2013-11-01T00:14:22.390Z
2013-11-07T05:36:31.690Z
2020-06-17T09:38:08.363Z
D.W.
<p>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&gt;b&gt;0$ and $\gcd(a,b)=1$.</p>

<p><strong>Definition.</strong> 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 &lt; b$ and $0 \le y &lt; a$.</p>

<p><strong>The basic idea.</strong>  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 <a href="https://en.wikipedia.org/wiki/Addition_chain" rel="nofollow">addition chain</a> that ends with $c$, and then use the preceding observations to calculate a representation of $c$.  Details follow below.</p>

<hr>

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

<p><strong>Proof sketch.</strong> 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 &lt; b$, $|y| &lt; a$, but this is only possible when $x=y=0$.</p>

<p><strong>Fact 2.</strong> We can find a representation of $1$, within 64-bit arithmetic.</p>

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

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

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

<p><strong>Fact 3.</strong> 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$.</p>

<p><strong>Proof sketch.</strong> 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 &lt; a$, then $0 \le x_3 &lt; b$.  Now the fact that $(x_3,y_3)$ is a representation of $c_1+c_2$ is immediate, by linearity.</p>

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

<hr>

<p>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 &lt; k$).  There are standard ways to choose such a sequence of length $\le 2 \lg c$.</p>

<p>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.</p>

<p>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.</p>

CC BY-SA 3.0
2013-11-11T17:17:06.323Z
I appreaciate your nice work. Do you think this could somehow be extended to the case $c > b$?
J. J.
CC BY-SA 3.0
2013-11-11T19:04:19.177Z
@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.
CC BY-SA 3.0
2013-11-08T05:07:08.620Z
2013-11-08T05:26:13.800Z
2013-11-08T05:26:13.800Z
D.W.
<p>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?</p>

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

<p>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.</p>

CC BY-SA 3.0
2013-10-30T22:36:14.577Z
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ε
CC BY-SA 3.0
2013-10-30T23:29:47.380Z
This might be more suitable for [scicomp.se] as it seems like a numerical analysis type question.
Kaveh
CC BY-SA 3.0
2013-10-31T08:15:07.257Z
@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.
CC BY-SA 3.0
2013-10-31T08:42:39.673Z
@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.
CC BY-SA 3.0
2013-10-31T23:51:56.023Z
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.
CC BY-SA 3.0
2013-11-07T04:42:29.080Z
@J.J. Is there a motivation for your question?
Turbo
CC BY-SA 3.0
2013-11-11T17:16:38.573Z
@JAS: Our team happened to run into an overflow problem with our linear diophantine solver when we did some practicing for programming competitions. The question sparked from there.
J. J.
<ds.algorithms>
Avoiding overflow in finding a solution to $ax + by = c$
J. J.
# Avoiding overflow in finding a solution to $ax + by = c$

(Asked by: username_0 on 2013-10-30)

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.

**Tags:** ds.algorithms

> **Comments:**
> **username_1**: 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?  
> **username_2**: This might be more suitable for [scicomp.se] as it seems like a numerical analysis type question.  
> **username_0**: @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|)$.  
> **username_0**: @username_2: 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.  
> **username_3**: 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.  
> **username_4**: @J.J. Is there a motivation for your question?  
> **username_0**: @JAS: Our team happened to run into an overflow problem with our linear diophantine solver when we did some practicing for programming competitions. The question sparked from there.  

---

## 3 answers

### Answer by username_0 on 2013-11-21. Score: 1 | Accepted answer

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.


---

### Answer by username_3 on 2013-11-01. Score: 3

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.

> **Comments:**
> **username_0**: 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$?  
> **username_3**: @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}$.  
> **username_3**: 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}$.  
> **username_4**: 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$?  
> **username_3**: @JAS, we can't work modulo $3^{50}$, because $3^{50} > 2^{64}$, and we're trying to work with just 64-bit arithmetic.  Same problem with for $5^{56}$.  My parameter choices were chosen so that $2^{\alpha}$, $3^{\alpha}$, and $5^{\alpha}$ are all smaller than $2^{64}$, so you can work modulo each of them using only 64-bit arithmetic.  I'll edit my answer to make this clearer -- I just noticed a typo in my answer that made this unnecessarily opaque.  
> **username_3**: @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.  

---

### Answer by username_3 on 2013-11-08. Score: 2

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](https://en.wikipedia.org/wiki/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

**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.

> **Comments:**
> **username_0**: I appreaciate your nice work. Do you think this could somehow be extended to the case $c > b$?  
> **username_3**: @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.  

---

