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.

  • 3
    $\begingroup$ 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? $\endgroup$
    – Jeffε
    Oct 30 '13 at 22:36
  • $\begingroup$ This might be more suitable for Computational Science as it seems like a numerical analysis type question. $\endgroup$
    – Kaveh
    Oct 30 '13 at 23:29
  • 1
    $\begingroup$ @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|)$. $\endgroup$
    – J. J.
    Oct 31 '13 at 8:15
  • $\begingroup$ @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. $\endgroup$
    – J. J.
    Oct 31 '13 at 8:42
  • $\begingroup$ 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. $\endgroup$
    – 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.

  • $\begingroup$ 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$? $\endgroup$
    – J. J.
    Nov 6 '13 at 19:05
  • $\begingroup$ @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}$. $\endgroup$
    – D.W.
    Nov 6 '13 at 19:19
  • $\begingroup$ 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}$. $\endgroup$
    – D.W.
    Nov 6 '13 at 19:22
  • $\begingroup$ 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$? $\endgroup$
    – Mr.
    Nov 7 '13 at 5:04
  • 1
    $\begingroup$ @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. $\endgroup$
    – 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.

  • $\begingroup$ I appreaciate your nice work. Do you think this could somehow be extended to the case $c > b$? $\endgroup$
    – J. J.
    Nov 11 '13 at 17:17
  • $\begingroup$ @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. $\endgroup$
    – 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.