∑MathIntermediate

Linear Diophantine Equations

Key Points

•
A linear Diophantine equation ax + by = c has integer solutions if and only if gcd(a, b) divides c.
•
The Extended Euclidean Algorithm finds one particular solution by producing Bezout coefficients for a and b.
•
All solutions form an infinite arithmetic progression (a lattice line): $x = x 0$ + $\frac{b}{g}$ t and $y = y 0$ - $\frac{a}{g}$ t where $g = g c d (a$ ,b).
•
Counting solutions in ranges reduces to finding the valid interval of the parameter t using careful floor/ceil arithmetic.
•
Systems of linear Diophantine equations can be reduced to a single equation in a parameter or to congruences solved by CRT-like reasoning.
•
Careful handling of signs, zero coefficients, and 64-bit overflow is essential in competitive programming implementations.
•
Runtime is dominated by Extended GCD, which runs in O(log min(|a|, |b|)) time and constant space.
•
Modular equations a x ≡ b (mod m) and CRT are direct corollaries of linear Diophantine theory and often appear in CF problems.

Prerequisites

→Greatest Common Divisor and Euclidean Algorithm — Understanding gcd and how to compute it efficiently is necessary for solvability and for the Extended Euclidean Algorithm.
→Modular Arithmetic Basics — Interpreting congruences and residues connects Diophantine equations to modular equations and CRT.
→Integer Division and Rounding — Correctly transforming interval constraints requires sign-safe floor and ceil operations.
→Overflow and Integer Types — Products and scaling can exceed 64-bit; using wider intermediates like __int128 prevents errors.
→Linear Algebra (2×2 Determinant) — Determinant indicates dependence of equations when solving systems and simplifies coefficients.

Detailed Explanation

Tap terms for definitions

01Overview

Linear Diophantine equations are equations of the form ax + by = c, where a, b, c are given integers and we seek integer solutions (x, y). These appear everywhere: from solving modular equations to aligning periodic events and resource allocation problems. The key property is that solutions exist exactly when the greatest common divisor g = gcd(a, b) divides c. When solutions exist, the Extended Euclidean Algorithm (EEA) gives one specific solution (x0, y0). From there, all solutions can be described compactly using a single integer parameter t, producing a line of lattice points spaced regularly. This structure lets us efficiently answer more advanced questions: find any solution, find a solution with minimal |x|, count how many solutions lie inside a box [xL, xR] × [yL, yR], or solve systems of two equations by eliminating variables or translating them into congruences solved via the Chinese Remainder Theorem (CRT). In competitive programming (e.g., CF 1500–1900), mastering these transformations and implementing robust arithmetic (including safe floor/ceil for negative numbers and overflow-aware multiplications) is crucial. The algorithms run in logarithmic time with O(1) memory, so even very large inputs (up to 64-bit) are feasible with careful coding.

02Intuition & Analogies

Imagine a grid of streets forming integer coordinate points. The equation ax + by = c describes a tilted straight line on this grid. We are interested only in points where streets cross (integer coordinates). If the line passes through at least one intersection, it will pass through a whole infinite chain of them at fixed intervals; if it misses, it misses all of them. Why? Think of gcd(a, b) as the "grid’s step size" along the direction defined by (a, b). You can only move along combinations of a steps in x and b steps in y, and the net value ax + by moves in chunks of size gcd(a, b). Therefore, c must line up with that chunking or else the line never hits an intersection. Once we find a single intersection point (x0, y0), sliding along the line by adding b/g to x and subtracting a/g from y keeps you on the same line because ax + by remains c. This slide is like moving from one equally spaced lamp post to the next along a straight road. If you fence off a rectangular area [xL, xR] × [yL, yR] and want to count how many lamp posts are inside, you don’t need to scan the whole grid—you just find the first t that lands in the box and the last t that still stays in, then count how many steps fit. For multiple lines (a system), you can parameterize one line and plug it into the other to reduce the problem to a single parameter equation, or translate one into a modular constraint and combine constraints using a CRT-like merge.

03Formal Definition

A linear Diophantine equation is ax +

b y = c

with a, b, c ∈ ℤ and unknowns x, y ∈ ℤ. Let

g = g c d (a

, b). A necessary and sufficient condition for solvability is g

∣ c . I f so l v ab l e, t h eree x i s t in t e g erss, t (B ezo u t coe ff i c i e n t s) s u c h t ha t a s + b t = g . T h e n (x 0, y 0) = (s \cdot \frac{c}{g}, t \cdot \frac{c}{g}) i s a p a r t i c u l a rso l u t i o n . T h eco m pl e t ese t o f so l u t i o n s i s g i v e nb y t h eo n e - p a r am e t er f ami l y : x = x 0 + \frac{b}{g} \cdot kan d y = y 0 - \frac{a}{g} \cdot k f or an y k \in Z . T hi s f or m s a cose t o f t h e l a tt i ce g e n er a t e d b y v ec t ors (b / g, - a / g) . F orco n s t r ain t s x \in [xL, x R], y \in [y L, y R], t h e p a r am e t er km u s tl i e in t h e in t ersec t i o n o f tw o in t e g er in t er v a l s d er i v e df ro m l in e a r in e q u a l i t i es . A sys t e m o f tw oe q u a t i o n s a 1 x + b 1 y = c 1 an d a 2 x + b 2 y = c 2 c anb e han d l e d b y f i rs tp a r am e t er i z in g so l u t i o n so f t h e f i rs t, t h e n s u b s t i t u t in g in t o t h eseco n d t o g e t a l in e a re q u a t i o nink; e q u i v a l e n tl y, d er i v e a co n g r u e n ce f or x or kan d so l v e i t u s in g CRTw h e na pp ro p r ia t e . Sp ec ia l c a ses in c l u d e m o d u l a r l in e a re q u a t i o n s a x \equiv r (m o d m), w hi c h corres p o n d t o a x + m y = r, w i t h so l u t i o n se x i s t in g i ff g c d (a, m) ∣

04When to Use

Use linear Diophantine methods when: (1) You must solve ax + by = c over integers, not reals—common in number theory and scheduling problems. (2) You need to solve modular equations ax ≡ b (mod m); convert to ax + m y = b and apply EEA. (3) You want to count or enumerate solutions within bounds; parameterize all solutions and intersect with range constraints. (4) You face two linear equations in two integer variables; parameterize one equation and insert into the other to find a common solution or detect inconsistency. (5) You need CRT merges of congruences, especially with non-coprime moduli; extended GCD and Diophantine reasoning unify the approach. (6) You are minimizing a linear cost (e.g., minimize |x|, minimize ax + by under constraints) over integer solutions; shift along the solution family to the best candidate. In CF problems (1500–1900), expect tasks like finding any solution, counting solutions in rectangles, merging congruences, or solving mixed constraints that translate to Diophantine equations. The tools here run in O(log n), so they are ideal for large 64-bit inputs when implemented carefully.

⚠️Common Mistakes

Ignoring the gcd condition: attempting to solve when gcd(a, b) ∤ c leads to bogus arithmetic. Always check divisibility first. - Mishandling zero coefficients: if a = 0 or b = 0, the equation reduces to a single-variable constraint; handle these early to avoid division by zero. - Wrong parameter shifts: using x = x0 + b·t (instead of b/g) or mixing signs in y’s update; always divide steps by g and remember y decreases by (a/g) when x increases by (b/g). - Incorrect floor/ceil with negative numbers: using naive integer division for ranges can shrink or expand intervals incorrectly. Implement robust floor_div and ceil_div that handle signs. - Overflow in 64-bit: products like a2b1 - b2a1 or scaling x0, y0 by c/g can overflow. Use 128-bit intermediates (__int128) and cast back after range checks. - Forgetting normalization in modular equations: solving ax ≡ b (mod m) requires reducing by g = gcd(a, m) then working modulo m/g. - CRT with non-coprime moduli: assuming a solution always exists; in reality, r1 ≡ r2 (mod gcd(m1, m2)) is required. - Not considering negative moduli or steps: ensure step sizes and divisions work regardless of sign; write sign-safe helpers. - Returning any solution when the problem asks for one in a specific range or minimal absolute value; remember to shift to meet constraints. - Skipping degeneracy checks in systems: if the determinant (a1 b2 − a2 b1) vanishes, equations may be dependent (infinite solutions) or inconsistent (no solution).

Key Formulas

GCD

g = g cd (a, b)

Explanation: The greatest common divisor g is the largest integer dividing both a and b. It sets the granularity of values that ax + by can achieve.

Solvability Criterion

a x + b y = c has solutions ⟺ g ∣ c

Explanation: The linear combination ax + by changes in steps of size g. Therefore c must be a multiple of g for integer solutions to exist.

Bezout Identity

a s + b t = g

Explanation: Extended GCD returns integers s, t such that this holds. Scaling by c/g yields a particular solution to ax + by = c.

General Solution

x = x_{0} + \frac{b}{g} t, y = y_{0} - \frac{a}{g} t, t \in Z

Explanation: All integer solutions lie on a line with step vector (b/g, −a/g). Starting from any one solution (x0, y0), varying t enumerates all.

t-interval for Box Constraints

t \in [⌈ \frac{x _{L} - x _{0}}{b / g} ⌉, ⌊ \frac{x _{R} - x _{0}}{b / g} ⌋] \cap [⌈ \frac{y _{0} - y _{R}}{a / g} ⌉, ⌊ \frac{y _{0} - y _{L}}{a / g} ⌋]

Explanation: Converting x and y range constraints into bounds on t yields the intersection interval. The number of integer t in this interval is the count of solutions in the box.

System Reduction Coefficient

K = \frac{a _{2} b _{1} - b _{2} a _{1}}{g _{1}}, g_{1} = g cd (a_{1}, b_{1})

Explanation: After parameterizing solutions of a1x + b1y = c1, substitution into a2x + b2y = c2 yields K t = RHS. This K equals the determinant divided by g1.

Parameter from Second Equation

If K \neq = 0 : t_{0} = \frac{c _{2} - ( a _{2} x _{0} + b _{2} y _{0} )}{K}

Explanation: A second equation pins down the parameter t uniquely if $K \neq = 0;$ solvability requires K | RHS. Then x and y follow from the general solution.

Modular-to-Diophantine

a x \equiv r (mod m) ⟺ a x + m y = r

Explanation: A modular equation is equivalent to a linear Diophantine equation in two variables. Existence requires gcd(a, m) | r.

CRT Consistency

CRT condition: r_{1} \equiv r_{2} (mod g cd (m_{1}, m_{2}))

Explanation: Two congruences can be merged if their residues agree modulo the gcd of their moduli. The combined modulus is lcm(m1, m2).

LCM-GCD Relation

lcm (m_{1}, m_{2}) = \frac{m _{1} m _{2}}{g cd ( m _{1} , m _{2} )}

Explanation: The least common multiple of two positive integers relates to their product and gcd. It is the modulus of the merged congruence in CRT.

Complexity Analysis

The core computational tool is the Extended Euclidean Algorithm (EEA), which runs in O(log min(

∣ a ∣

∣ b ∣

)) time and O(1) space. This bound comes from the fact that each recursive or iterative step reduces the remainder size similarly to the Euclidean Algorithm; the number of steps is proportional to the number of digits (in base 2) of the smaller input. Finding a particular solution to ax +

b y = c

involves a single EEA call to compute (g, s, t) such that as +

b t = g

, followed by constant-time arithmetic to scale by c/g. Shifting to the general solution

x = x 0

\frac{b}{g}

y = y 0

−

\frac{a}{g}

t is algebraic and O(1). Counting solutions within ranges uses only a few floor/ceil operations and comparisons, also O(1). Solving a system of two linear Diophantine equations via parameterization reduces to one EEA call plus constant-time operations; in degenerate cases, it remains O(1) after the EEA. Modular equation solving (ax ≡ b mod m) requires one EEA; CRT merging of two congruences likewise needs one EEA and O(1) arithmetic. In all cases, space complexity is O(1), as only a fixed number of 64-bit (and sometimes 128-bit intermediates) variables are used. Practical performance considerations include avoiding 64-bit overflow during intermediate products (use __int128), implementing sign-correct floor/ceil division, and normalizing by gcd early to keep numbers small. These safeguards do not change asymptotic complexity but are essential for correctness and robustness on large inputs (up to about 10^18).

Code Examples

Find any solution to ax + by = c and enumerate all solutions (Extended GCD)

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 struct EGResult {
5     long long g; // gcd(a, b)
6     long long x; // coefficient for a
7     long long y; // coefficient for b
8 };
9 
10 // Extended Euclidean Algorithm: returns g = gcd(a, b) and x, y with ax + by = g
11 EGResult extended_gcd(long long a, long long b) {
12     if (b == 0) {
13         return { llabs(a), (a >= 0 ? 1 : -1), 0 };
14     }
15     EGResult t = extended_gcd(b, a % b);
16     long long g = t.g;
17     long long x1 = t.y;
18     long long y1 = t.x - (a / b) * t.y;
19     return { g, x1, y1 };
20 }
21 
22 // Try to find any solution to ax + by = c. Returns false if none.
23 bool find_any_solution(long long a, long long b, long long c,
24                        long long &x0, long long &y0, long long &g) {
25     if (a == 0 && b == 0) {
26         if (c == 0) { x0 = y0 = 0; g = 0; return true; } // all (x,y) are solutions
27         return false; // 0x + 0y = c has no solution if c != 0
28     }
29     EGResult eg = extended_gcd(a, b);
30     g = eg.g;
31     if (c % g != 0) return false;
32     // Scale the Bezout coefficients by c/g. Use 128-bit temporarily to avoid overflow.
33     __int128 scale = c / g;
34     __int128 xi = (__int128)eg.x * scale;
35     __int128 yi = (__int128)eg.y * scale;
36     x0 = (long long)xi;
37     y0 = (long long)yi;
38     return true;
39 }
40 
41 // Shift a known solution by cnt steps along the solution line.
42 void shift_solution(long long &x, long long &y, long long a, long long b, long long g, long long cnt) {
43     // x += cnt * (b/g), y -= cnt * (a/g)
44     __int128 nx = (__int128)x + (__int128)cnt * (b / g);
45     __int128 ny = (__int128)y - (__int128)cnt * (a / g);
46     x = (long long)nx;
47     y = (long long)ny;
48 }
49 
50 int main() {
51     ios::sync_with_stdio(false);
52     cin.tie(nullptr);
53 
54     long long a, b, c; 
55     // Example: solve 30x + 18y = 6
56     a = 30; b = 18; c = 6;
57     long long x0, y0, g;
58     if (!find_any_solution(a, b, c, x0, y0, g)) {
59         cout << "No solution\n";
60         return 0;
61     }
62     cout << "One solution: x0=" << x0 << ", y0=" << y0 << " (gcd=" << g << ")\n";
63     // Enumerate a few solutions by shifting t = -2..2
64     for (long long t = -2; t <= 2; ++t) {
65         long long x = x0, y = y0;
66         shift_solution(x, y, a, b, g, t);
67         cout << "t=" << t << ": x=" << x << ", y=" << y << "; check a*x+b*y=" << (a*x + b*y) << "\n";
68     }
69     return 0;
70 }
71

Uses Extended GCD to compute a particular solution to ax + by = c and shows how to generate all solutions by stepping along the lattice line using the parameter t. The shift uses the canonical step vector (b/g, −a/g).

Time: O(log min(|a|, |b|))Space: O(1)

Count solutions of ax + by = c inside [xL, xR] × [yL, yR]

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 struct EGResult { long long g, x, y; };
5 
6 EGResult extended_gcd(long long a, long long b) {
7     if (b == 0) return { llabs(a), (a >= 0 ? 1 : -1), 0 };
8     EGResult t = extended_gcd(b, a % b);
9     return EGResult{ t.g, t.y, t.x - (a / b) * t.y };
10 }
11 
12 bool find_any_solution(long long a, long long b, long long c,
13                        long long &x0, long long &y0, long long &g) {
14     if (a == 0 && b == 0) { if (c == 0) { x0 = y0 = 0; g = 0; return true; } else return false; }
15     EGResult eg = extended_gcd(a, b); g = eg.g;
16     if (c % g != 0) return false;
17     __int128 scale = c / g;
18     x0 = (long long)((__int128)eg.x * scale);
19     y0 = (long long)((__int128)eg.y * scale);
20     return true;
21 }
22 
23 // Floor division a // b towards -infinity, handles negatives correctly.
24 long long floor_div(long long a, long long b) {
25     assert(b != 0);
26     long long q = a / b; long long r = a % b;
27     if (r != 0 && ((r > 0) != (b > 0))) --q; // adjust when trunc towards 0 differs from floor
28     return q;
29 }
30 
31 // Ceil division towards +infinity, handles negatives correctly.
32 long long ceil_div(long long a, long long b) {
33     assert(b != 0);
34     long long q = a / b; long long r = a % b;
35     if (r != 0 && ((r > 0) == (b > 0))) ++q;
36     return q;
37 }
38 
39 // Count integer solutions (x, y) to ax + by = c with x in [xL,xR], y in [yL,yR]
40 long long count_in_box(long long a, long long b, long long c,
41                        long long xL, long long xR, long long yL, long long yR) {
42     if (xL > xR || yL > yR) return 0;
43     long long x0, y0, g;
44     if (!find_any_solution(a, b, c, x0, y0, g)) return 0;
45     if (g == 0) { // a = b = c = 0: whole plane; count is (xR-xL+1)*(yR-yL+1)
46         __int128 cnt = (__int128)(xR - xL + 1) * (yR - yL + 1);
47         if (cnt > LLONG_MAX) return LLONG_MAX; // saturated for safety
48         return (long long)cnt;
49     }
50     long long stepX = b / g; // delta x per t
51     long long stepY = a / g; // delta y per t (note: y = y0 - stepY * t)
52 
53     // Derive t interval from x bounds: x = x0 + stepX * t ∈ [xL, xR]
54     long long tLx = ceil_div(xL - x0, stepX);
55     long long tRx = floor_div(xR - x0, stepX);
56     if (tLx > tRx) return 0; // no t satisfies x-bounds
57 
58     // Derive t interval from y bounds: y = y0 - stepY * t ∈ [yL, yR]
59     // Equivalent to y0 - yR ≤ stepY * t ≤ y0 - yL
60     long long tLy = ceil_div(y0 - yR, stepY);
61     long long tRy = floor_div(y0 - yL, stepY);
62 
63     long long tL = max(tLx, tLy);
64     long long tR = min(tRx, tRy);
65     if (tL > tR) return 0;
66     // number of integers in [tL, tR]
67     __int128 cnt = (__int128)(tR - tL + 1);
68     if (cnt > LLONG_MAX) return LLONG_MAX;
69     return (long long)cnt;
70 }
71 
72 int main() {
73     ios::sync_with_stdio(false);
74     cin.tie(nullptr);
75 
76     // Example: Count solutions to 30x + 18y = 6 with 0 ≤ x ≤ 10, 0 ≤ y ≤ 10
77     long long a = 30, b = 18, c = 6;
78     long long xL = 0, xR = 10, yL = 0, yR = 10;
79     cout << count_in_box(a, b, c, xL, xR, yL, yR) << "\n";
80 
81     return 0;
82 }
83

Counts solutions by transforming the box constraints into an interval for the parameter t. It uses sign-safe floor/ceil division so that negative steps and values are handled correctly. Saturation is shown for the fully degenerate a = b = c = 0 case.

Time: O(log min(|a|, |b|)) for the EEA; O(1) arithmetic otherwiseSpace: O(1)

Solve a system of two linear Diophantine equations

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 struct EGResult { long long g, x, y; };
5 
6 EGResult extended_gcd(long long a, long long b) {
7     if (b == 0) return { llabs(a), (a >= 0 ? 1 : -1), 0 };
8     EGResult t = extended_gcd(b, a % b);
9     return { t.g, t.y, t.x - (a / b) * t.y };
10 }
11 
12 bool find_any_solution(long long a, long long b, long long c,
13                        long long &x0, long long &y0, long long &g) {
14     if (a == 0 && b == 0) { if (c == 0) { x0 = y0 = 0; g = 0; return true; } else return false; }
15     EGResult eg = extended_gcd(a, b); g = eg.g;
16     if (c % g != 0) return false;
17     __int128 scale = c / g;
18     x0 = (long long)((__int128)eg.x * scale);
19     y0 = (long long)((__int128)eg.y * scale);
20     return true;
21 }
22 
23 // Solve the system:
24 //   a1 x + b1 y = c1
25 //   a2 x + b2 y = c2
26 // Returns (true, x, y) if a solution exists; handles degenerate cases.
27 bool solve_two_diophantine(long long a1, long long b1, long long c1,
28                            long long a2, long long b2, long long c2,
29                            long long &x, long long &y) {
30     // Step 1: Solve the first equation
31     long long x0, y0, g1;
32     if (!find_any_solution(a1, b1, c1, x0, y0, g1)) return false;
33 
34     // General solution of the first: x = x0 + (b1/g1)t, y = y0 - (a1/g1)t
35     // Substitute into the second equation to solve for t:
36     // a2*(x0 + (b1/g1)t) + b2*(y0 - (a1/g1)t) = c2
37     // ((a2*b1 - b2*a1)/g1) * t + (a2*x0 + b2*y0) = c2
38 
39     __int128 det = (__int128)a2 * b1 - (__int128)b2 * a1; // a2*b1 - b2*a1
40     long long K; // K = det / g1 (must be divisible)
41     if (det % g1 != 0) return false; // no way to get integer t
42     K = (long long)(det / g1);
43 
44     __int128 rhs = (__int128)c2 - ((__int128)a2 * x0 + (__int128)b2 * y0);
45 
46     if (K == 0) {
47         // Second equation becomes (a2*x0 + b2*y0) = c2; rhs must be 0
48         if (rhs != 0) return false; // inconsistent parallel lines
49         // Infinitely many solutions along t; return one (t = 0)
50         x = x0; y = y0; return true;
51     }
52 
53     if (rhs % K != 0) return false; // t must be integer
54     long long t = (long long)(rhs / K);
55 
56     // Compute final x, y
57     __int128 xi = (__int128)x0 + (__int128)(b1 / g1) * t;
58     __int128 yi = (__int128)y0 - (__int128)(a1 / g1) * t;
59     x = (long long)xi; y = (long long)yi;
60     return true;
61 }
62 
63 int main() {
64     ios::sync_with_stdio(false);
65     cin.tie(nullptr);
66 
67     // Example system:
68     //  6x + 9y = 3
69     // 10x + 7y = 1
70     long long a1=6,b1=9,c1=3; long long a2=10,b2=7,c2=1;
71     long long x,y;
72     if (solve_two_diophantine(a1,b1,c1,a2,b2,c2,x,y)) {
73         cout << "x=" << x << ", y=" << y << "\n";
74         cout << "Check1: " << a1*x + b1*y << " == " << c1 << "\n";
75         cout << "Check2: " << a2*x + b2*y << " == " << c2 << "\n";
76     } else {
77         cout << "No solution\n";
78     }
79     return 0;
80 }
81

Parameterizes all solutions of the first equation, substitutes into the second, and solves for the parameter t. This mirrors a CRT-like reduction where the second equation pins down the residue class of the parameter. Degenerate cases with zero determinant are handled explicitly.

Time: O(log min(|a1|, |b1|))Space: O(1)

Solve modular linear equations and merge two congruences (CRT-like)

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 struct EGResult { long long g, x, y; };
5 
6 EGResult extended_gcd(long long a, long long b) {
7     if (b == 0) return { llabs(a), (a >= 0 ? 1 : -1), 0 };
8     EGResult t = extended_gcd(b, a % b);
9     return EGResult{ t.g, t.y, t.x - (a / b) * t.y };
10 }
11 
12 // Solve a x ≡ b (mod m). Returns (ok, x0, mod), where all solutions are x ≡ x0 (mod mod).
13 tuple<bool, long long, long long> solve_linear_mod(long long a, long long b, long long m) {
14     EGResult eg = extended_gcd(a, m);
15     long long g = eg.g;
16     if (b % g != 0) return {false, 0, 0};
17     __int128 scale = b / g;
18     // Minimal non-negative solution
19     __int128 x0 = ((__int128)eg.x * scale) % (m / g);
20     if (x0 < 0) x0 += (m / g);
21     return {true, (long long)x0, m / g};
22 }
23 
24 // Merge x ≡ r1 (mod m1) and x ≡ r2 (mod m2). Returns (ok, r, lcm)
25 // Works with non-coprime moduli.
26 tuple<bool, long long, long long> crt_merge(long long r1, long long m1, long long r2, long long m2) {
27     // Normalize residues
28     r1 %= m1; if (r1 < 0) r1 += m1;
29     r2 %= m2; if (r2 < 0) r2 += m2;
30 
31     EGResult eg = extended_gcd(m1, m2);
32     long long g = eg.g;
33     if ((r2 - r1) % g != 0) return {false, 0, 0}; // inconsistency
34 
35     // Solve m1 * k ≡ (r2 - r1) (mod m2)
36     __int128 k = (__int128)(r2 - r1) / g;
37     __int128 mod2_g = m2 / g;
38     __int128 inv = eg.x; // inverse of m1/g modulo m2/g
39     __int128 t = (k * inv) % mod2_g;
40     if (t < 0) t += mod2_g;
41 
42     __int128 lcm = (__int128)m1 / g * m2;
43     __int128 r = (__int128)r1 + t * m1;
44     r %= lcm; if (r < 0) r += lcm;
45     return {true, (long long)r, (long long)lcm};
46 }
47 
48 int main() {
49     ios::sync_with_stdio(false);
50     cin.tie(nullptr);
51 
52     // Example 1: Solve 12 x ≡ 18 (mod 30)
53     {
54         auto [ok, x0, mod] = solve_linear_mod(12, 18, 30);
55         if (ok) cout << "x ≡ " << x0 << " (mod " << mod << ")\n";
56         else cout << "No solution\n";
57     }
58 
59     // Example 2: Merge x ≡ 2 (mod 6) and x ≡ 3 (mod 15)
60     {
61         auto [ok, r, m] = crt_merge(2, 6, 3, 15);
62         if (ok) cout << "x ≡ " << r << " (mod " << m << ")\n";
63         else cout << "Inconsistent congruences\n";
64     }
65 
66     return 0;
67 }
68

Shows two core corollaries: solving a modular linear equation via Diophantine methods, and merging two congruences using a generalized CRT that works for non-coprime moduli. This is the CRT-like technique often needed when reducing systems.

Time: O(log max(m1, m2)) for CRT merge; O(log m) for a single modular equationSpace: O(1)

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	struct EGResult {
5	long long g; // gcd(a, b)
6	long long x; // coefficient for a
7	long long y; // coefficient for b
8	};
9
10	// Extended Euclidean Algorithm: returns g = gcd(a, b) and x, y with ax + by = g
11	EGResult extended_gcd(long long a, long long b) {
12	if (b == 0) {
13	return { llabs(a), (a >= 0 ? 1 : -1), 0 };
14	}
15	EGResult t = extended_gcd(b, a % b);
16	long long g = t.g;
17	long long x1 = t.y;
18	long long y1 = t.x - (a / b) * t.y;
19	return { g, x1, y1 };
20	}
21
22	// Try to find any solution to ax + by = c. Returns false if none.
23	bool find_any_solution(long long a, long long b, long long c,
24	long long &x0, long long &y0, long long &g) {
25	if (a == 0 && b == 0) {
26	if (c == 0) { x0 = y0 = 0; g = 0; return true; } // all (x,y) are solutions
27	return false; // 0x + 0y = c has no solution if c != 0
28	}
29	EGResult eg = extended_gcd(a, b);
30	g = eg.g;
31	if (c % g != 0) return false;
32	// Scale the Bezout coefficients by c/g. Use 128-bit temporarily to avoid overflow.
33	__int128 scale = c / g;
34	__int128 xi = (__int128)eg.x * scale;
35	__int128 yi = (__int128)eg.y * scale;
36	x0 = (long long)xi;
37	y0 = (long long)yi;
38	return true;
39	}
40
41	// Shift a known solution by cnt steps along the solution line.
42	void shift_solution(long long &x, long long &y, long long a, long long b, long long g, long long cnt) {
43	// x += cnt * (b/g), y -= cnt * (a/g)
44	__int128 nx = (__int128)x + (__int128)cnt * (b / g);
45	__int128 ny = (__int128)y - (__int128)cnt * (a / g);
46	x = (long long)nx;
47	y = (long long)ny;
48	}
49
50	int main() {
51	ios::sync_with_stdio(false);
52	cin.tie(nullptr);
53
54	long long a, b, c;
55	// Example: solve 30x + 18y = 6
56	a = 30; b = 18; c = 6;
57	long long x0, y0, g;
58	if (!find_any_solution(a, b, c, x0, y0, g)) {
59	cout << "No solution\n";
60	return 0;
61	}
62	cout << "One solution: x0=" << x0 << ", y0=" << y0 << " (gcd=" << g << ")\n";
63	// Enumerate a few solutions by shifting t = -2..2
64	for (long long t = -2; t <= 2; ++t) {
65	long long x = x0, y = y0;
66	shift_solution(x, y, a, b, g, t);
67	cout << "t=" << t << ": x=" << x << ", y=" << y << "; check ax+by=" << (ax + by) << "\n";
68	}
69	return 0;
70	}
71

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	struct EGResult { long long g, x, y; };
5
6	EGResult extended_gcd(long long a, long long b) {
7	if (b == 0) return { llabs(a), (a >= 0 ? 1 : -1), 0 };
8	EGResult t = extended_gcd(b, a % b);
9	return EGResult{ t.g, t.y, t.x - (a / b) * t.y };
10	}
11
12	bool find_any_solution(long long a, long long b, long long c,
13	long long &x0, long long &y0, long long &g) {
14	if (a == 0 && b == 0) { if (c == 0) { x0 = y0 = 0; g = 0; return true; } else return false; }
15	EGResult eg = extended_gcd(a, b); g = eg.g;
16	if (c % g != 0) return false;
17	__int128 scale = c / g;
18	x0 = (long long)((__int128)eg.x * scale);
19	y0 = (long long)((__int128)eg.y * scale);
20	return true;
21	}
22
23	// Floor division a // b towards -infinity, handles negatives correctly.
24	long long floor_div(long long a, long long b) {
25	assert(b != 0);
26	long long q = a / b; long long r = a % b;
27	if (r != 0 && ((r > 0) != (b > 0))) --q; // adjust when trunc towards 0 differs from floor
28	return q;
29	}
30
31	// Ceil division towards +infinity, handles negatives correctly.
32	long long ceil_div(long long a, long long b) {
33	assert(b != 0);
34	long long q = a / b; long long r = a % b;
35	if (r != 0 && ((r > 0) == (b > 0))) ++q;
36	return q;
37	}
38
39	// Count integer solutions (x, y) to ax + by = c with x in [xL,xR], y in [yL,yR]
40	long long count_in_box(long long a, long long b, long long c,
41	long long xL, long long xR, long long yL, long long yR) {
42	if (xL > xR \|\| yL > yR) return 0;
43	long long x0, y0, g;
44	if (!find_any_solution(a, b, c, x0, y0, g)) return 0;
45	if (g == 0) { // a = b = c = 0: whole plane; count is (xR-xL+1)*(yR-yL+1)
46	__int128 cnt = (__int128)(xR - xL + 1) * (yR - yL + 1);
47	if (cnt > LLONG_MAX) return LLONG_MAX; // saturated for safety
48	return (long long)cnt;
49	}
50	long long stepX = b / g; // delta x per t
51	long long stepY = a / g; // delta y per t (note: y = y0 - stepY * t)
52
53	// Derive t interval from x bounds: x = x0 + stepX * t ∈ [xL, xR]
54	long long tLx = ceil_div(xL - x0, stepX);
55	long long tRx = floor_div(xR - x0, stepX);
56	if (tLx > tRx) return 0; // no t satisfies x-bounds
57
58	// Derive t interval from y bounds: y = y0 - stepY * t ∈ [yL, yR]
59	// Equivalent to y0 - yR ≤ stepY * t ≤ y0 - yL
60	long long tLy = ceil_div(y0 - yR, stepY);
61	long long tRy = floor_div(y0 - yL, stepY);
62
63	long long tL = max(tLx, tLy);
64	long long tR = min(tRx, tRy);
65	if (tL > tR) return 0;
66	// number of integers in [tL, tR]
67	__int128 cnt = (__int128)(tR - tL + 1);
68	if (cnt > LLONG_MAX) return LLONG_MAX;
69	return (long long)cnt;
70	}
71
72	int main() {
73	ios::sync_with_stdio(false);
74	cin.tie(nullptr);
75
76	// Example: Count solutions to 30x + 18y = 6 with 0 ≤ x ≤ 10, 0 ≤ y ≤ 10
77	long long a = 30, b = 18, c = 6;
78	long long xL = 0, xR = 10, yL = 0, yR = 10;
79	cout << count_in_box(a, b, c, xL, xR, yL, yR) << "\n";
80
81	return 0;
82	}
83