Linear Diophantine Equations, in depth: existence, every solution from one, counting in a range, non-negative solutions and overflow

The geometry in one picture

When a solution exists

Let g = gcd(a, b). Every integer of the form a x + b y is a multiple of g, because g divides both a and b. So if g does not divide c, there is no solution, and that is the end. 6x + 9y = 20 fails immediately: gcd is 3, and 3 does not divide 20.

The converse is the Bezout identity: there always exist integers x and y with a x + b y = g, and the extended Euclid algorithm finds them. Multiply both sides by c / g and you have a solution of the original equation. So the complete existence test is one gcd and one remainder: a x + b y = c has integer solutions if and only if gcd(a, b) divides c (with a and b not both zero).

Extended Euclid in code

The Euclid algorithm replaces (a, b) by (b, a mod b) until the remainder is zero; the last non-zero remainder is the gcd. The extended version also carries coefficients that express each remainder as a combination of the original a and b. The iterative form avoids recursion and makes the invariant easy to see: at every step old_r == a*old_x + b*old_y and r == a*x + b*y.

def ext_gcd(a, b):
    '''Return (g, x, y) with a*x + b*y == g == gcd(a, b) >= 0.'''
    old_r, r = a, b
    old_x, x = 1, 0
    old_y, y = 0, 1
    while r != 0:
        q = old_r // r
        old_r, r = r, old_r - q * r
        old_x, x = x, old_x - q * x
        old_y, y = y, old_y - q * y
    if old_r < 0:                       # negative inputs can leave a negative gcd
        old_r, old_x, old_y = -old_r, -old_x, -old_y
    return old_r, old_x, old_y

For a = 35 and b = 55 the remainders run 35, 55, 35, 20, 15, 5, 0, so gcd is 5. Back-substituting: 5 = 20 - 15, 15 = 35 - 20, 20 = 55 - 35, so 5 = 2 times 20 - 35 = 2 times 55 - 3 times 35. The function returns (5, -3, 2): 35 times -3 plus 55 times 2 equals 5. The number of iterations is logarithmic in the smaller input, bounded by about five times its decimal digit count, so the solver costs nothing even for 64-bit values. The number theory algorithms overview proves the bound and uses the same routine for modular inverses.

From one solution to all of them

One solution gives all of them. Suppose (x0, y0) solves a x + b y = c. For any other solution (x, y), subtract: a (x - x0) = -b (y - y0). Divide by g, writing a' = a / g and b' = b / g, which are coprime: a' (x - x0) = -b' (y - y0). Since b' divides the right side and shares no factor with a', it must divide x - x0. So x = x0 + k b' for some integer k, and substituting back gives y = y0 - k a'. Every such pair is a solution, and there are no others:

x = x0 + k * (b // g)
y = y0 - k * (a // g)        # for every integer k

Geometrically, the real solutions form a line, and the integer solutions are lattice points on it, evenly spaced by the vector (b/g, -a/g). That step is the smallest possible; dividing by g is what makes it so. Forgetting the division, and stepping by (b, -a), skips all but every g-th solution, which is a common bug that only shows up when g is greater than 1.

Solving and counting in a range

def solve(a, b, c):
    '''Solve a*x + b*y = c for nonzero a, b.
    Returns (x0, y0, dx, dy) meaning all solutions are (x0 + k*dx, y0 - k*dy), or None.'''
    g, x, y = ext_gcd(a, b)
    if c % g != 0:
        return None
    m = c // g
    return x * m, y * m, b // g, a // g


def ceil_div(p, q):
    return -((-p) // q)


def count_in_box(a, b, c, x_lo, x_hi, y_lo, y_hi):
    '''How many solutions have x_lo <= x <= x_hi and y_lo <= y <= y_hi.'''
    s = solve(a, b, c)
    if s is None:
        return 0
    x0, y0, dx, dy = s

    def k_range(base, step, lo, hi):      # all k with lo <= base + k*step <= hi
        if step > 0:
            return ceil_div(lo - base, step), (hi - base) // step
        return ceil_div(hi - base, step), (lo - base) // step

    k1, k2 = k_range(x0, dx, x_lo, x_hi)
    k3, k4 = k_range(y0, -dy, y_lo, y_hi)
    return max(0, min(k2, k4) - max(k1, k3) + 1)

Each bound on x or y becomes an interval of k, because x and y are linear in k. The answer is the size of the intersection of two intervals. The direction of each inequality depends on the sign of the step, which is why k_range branches: with negative a or b, dx or dy is negative, and dividing an inequality by a negative number flips it. Python's // is floor division for negatives, and ceil_div is built from it; in C, Java, Go or Rust, / truncates toward zero, so you must write explicit floor and ceiling helpers or the counts are off by one for negative intermediates.

Worked example: filling a truck

A warehouse ships orders in crates of 35 kg and 55 kg and must fill a truck to exactly 1,000 kg. How many crates of each?

Existence. gcd(35, 55) = 5, and 5 divides 1000, so there are solutions.

Particular solution. From above, 35 times -3 plus 55 times 2 = 5. Multiply by 1000 / 5 = 200: x0 = -600, y0 = 400. Check: 35 times -600 = -21,000 and 55 times 400 = 22,000, total 1,000.

General solution. b / g = 11 and a / g = 7, so x = -600 + 11k and y = 400 - 7k.

Constraints. Crates cannot be negative. x at least 0 needs 11k at least 600, so k at least 55 (600 / 11 is about 54.5). y at least 0 needs 7k at most 400, so k at most 57 (400 / 7 is about 57.1). So k is 55, 56 or 57, giving (5, 15), (16, 8) and (27, 1): 175 + 825, 560 + 440 and 945 + 55, each 1,000 kg. count_in_box(35, 55, 1000, 0, 10**9, 0, 10**9) returns 3.

Choosing among them. Once the solutions are a single parameter k, optimisation is easy. Total crate count is x + y = -200 + 4k, which increases with k, so the fewest crates is k = 55: 5 small and 15 large. Any objective linear in x and y is linear in k, so the optimum is at one end of the k interval; no search is needed.

Non-negative solutions and the coin problem

Non-negative solutions are the coin problem: with coins of a and b, which totals c can you pay exactly? For coprime a and b, Sylvester's result gives a sharp answer. The largest total that cannot be paid is a b - a - b, and exactly (a - 1)(b - 1) / 2 non-negative totals cannot be paid. With 3 and 5, the largest impossible total is 7 and the impossible totals are 1, 2, 4 and 7, which is four of them, matching (2 times 4) / 2. Every total from 8 upward is payable. For non-coprime a and b, divide by g first: only multiples of g are reachable, and the same result applies to a / g and b / g.

The count of non-negative solutions for a given c also follows from the k interval and grows roughly like c g / (a b). With three or more denominations no closed form exists for the largest impossible total in general, and you move to dynamic programming; see the two coin change variants and knapsack.

Where the equation shows up

Modular inverses. Solving a x ≡ 1 (mod m) is the Diophantine equation a x + m y = 1, which is solvable exactly when gcd(a, m) = 1. The x from ext_gcd, reduced mod m, is the inverse. This is how RSA private exponents are computed, alongside modular exponentiation.

Aligning periodic events. Two jobs start at offsets r1 and r2 with periods p1 and p2. They coincide when r1 + p1 i = r2 + p2 j, that is p1 i - p2 j = r2 - r1: a Diophantine equation that is solvable only if gcd(p1, p2) divides the offset difference, and whose solutions then repeat every lcm(p1, p2). This is the two-congruence Chinese remainder theorem in another form.

Strided memory access. If thread t reads address t s and memory has n banks, threads t and t' collide when s (t - t') is a multiple of n. The colliding pairs are solutions of a linear Diophantine equation, and n consecutive threads touch only n / gcd(s, n) distinct banks, gcd(s, n) threads on each. This is why GPU kernels pad shared-memory rows to an odd stride: making gcd(s, n) = 1 removes the conflicts.

Failure modes

• Not dividing by g. Stepping by (b, -a) instead of (b/g, -a/g) misses solutions whenever gcd is greater than 1.
• Truncating division. C-family / and % round toward zero; ceil and floor of negative quotients come out wrong, and so do counts. Python's // floors, so code ported from Python to C breaks silently.
• Overflow. x0 = x times c / g can be far larger than a, b or c: the Bezout coefficient is up to about b / g and c / g can be up to 2^63. Reduce first: compute x0 as x times ((c / g) mod (b / g)), taken mod b / g, so it stays below b / g, then derive y0 = (c - a x0) / b, using a mulmod or 128-bit intermediates for the products, which can still approach (b / g) squared. Python integers do not overflow, which hides this in prototypes.
• Zero coefficients. If a = 0 the equation is b y = c: y is fixed if b divides c, and x is free. If both are zero, either everything or nothing is a solution. Handle these before counting: the general solver still returns a solution, but one step component is zero and the range counter would divide by it.
• Negative gcd. Some extended Euclid implementations return a negative g for negative inputs, flipping the sign of the step. Normalise g to be positive, as the code above does.

What to do next

1. Implement ext_gcd and solve exactly as above, in your language, with explicit floor and ceiling division helpers.
2. Write a brute-force test over small a, b, c and boxes, including negatives and gcd greater than 1, and compare counts.
3. Add the zero-coefficient and both-zero cases as explicit branches with their own tests.
4. Check overflow: test with coefficients near your integer type's limit, or reduce c / g modulo b / g before multiplying.
5. When you need non-negative or bounded solutions, convert the bounds to a k interval rather than searching.
6. For three or more variables, reduce one pair to its gcd, solve, and substitute back, or switch to dynamic programming for counting.