 <p><a href="http://en.wikipedia.org/wiki/Diophantine_equation#Linear_Diophantine_equations" rel="nofollow noreferrer">Linear Diophantine</a> equations take the form <code>ax + by = c</code>. If <code>c</code> is the greatest common divisor of <code>a</code> and <code>b</code> this means <code>a=z'c</code> and <code>b=z''c</code> then this is <a href="http://en.wikipedia.org/wiki/B%C3%A9zout%27s_identity" rel="nofollow noreferrer">Bézout's identity</a> of the form</p> <p><img src="https://i.stack.imgur.com/87Ulj.png" alt="enter image description here"></p> <p>with <code>a=z'</code> and <code>b=z''</code> and the equation has an infinite number of solutions. So instead of trial searching method you can check <strong>if <code>c</code> is the <em>greatest common divisor</em> (GCD) of <code>a</code> and <code>b</code></strong> (in your case this translates into <code>bx - dy = c - a</code>)</p> <p>If indeed <code>a</code> and <code>b</code> are multiples of <code>c</code> then <code>x</code> and <code>y</code> can be computed using <a href="http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm" rel="nofollow noreferrer">extended Euclidean algorithm</a> which finds integers <code>x</code> and <code>y</code> (one of which is typically negative) that satisfy Bézout's identity</p> <p><img src="https://i.stack.imgur.com/ULthf.png" alt="enter image description here"></p> <p>and your answer is:</p> <p><code>a = k*x</code>, <code>b = k*y</code>, <code>c - a = k * gcd(a,b)</code> for any integer <strong>k</strong>.</p> <p>(as a side note: this holds also for any other <a href="http://en.wikipedia.org/wiki/Euclidean_domain" rel="nofollow noreferrer">Euclidean domain</a>, i.e. polynomial ring &amp; every Euclidean domain is <a href="http://en.wikipedia.org/wiki/Unique_factorization_domain" rel="nofollow noreferrer">unique factorization domain</a>). You can use <strong>Iterative Method</strong> to find these solutions:</p> <hr> <h2>Iterative method</h2> <p>By routine algebra of expanding and grouping like terms (refer to last section of <a href="http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm" rel="nofollow noreferrer">wikipedia article</a> mentioned before), the following algorithm for iterative method is obtained:</p> <ul> <li>1 . Apply Euclidean algorithm, and let qn (n starts from 1) be a finite list of quotients in the division.</li> <li>2 . Initialize x0, x1 as 1, 0, and y0, y1 as 0,1 respectively. <ul> <li>2.1 Then for each i so long as qi is defined,</li> <li>2.2 Compute xi+1 = xi−1 − qixi</li> <li>2.3 Compute yi+1 = yi−1 − qiyi</li> <li>2.4 Repeat the above after incrementing i by 1.</li> </ul></li> <li>3 . The answers are the second-to-last of xn and yn.</li> </ul> <p>pseudocode:</p> <pre><code>function extended_gcd(a, b) x := 0 lastx := 1 y := 1 lasty := 0 while b ≠ 0 quotient := a div b (a, b) := (b, a mod b) (x, lastx) := (lastx - quotient*x, x) (y, lasty) := (lasty - quotient*y, y) return (lastx, lasty) </code></pre> <p>So I have written example algorithm which calculates <em>greatest common divisor</em> using <strong>Euclidean Algorithm</strong> iterative method for non-negative <code>a</code> and <code>b</code> (for negative - <a href="http://math.stackexchange.com/questions/37806/extended-euclidean-algorithm-with-negative-numbers">these</a> extra steps are needed), it returns GCD and stores solutions for <code>x</code> and <code>y</code> in variables passed to it by reference:</p> <pre><code>int gcd_iterative(int a, int b, int&amp; x, int&amp; y) { int c; std::vector&lt;int&gt; r, q, x_coeff, y_coeff; x_coeff.push_back(1); y_coeff.push_back(0); x_coeff.push_back(0); y_coeff.push_back(1); if ( b == 0 ) return a; while ( b != 0 ) { c = b; q.push_back(a/b); r.push_back(b = a % b); a = c; x_coeff.push_back( *(x_coeff.end()-2) -(q.back())*x_coeff.back()); y_coeff.push_back( *(y_coeff.end()-2) -(q.back())*y_coeff.back()); } if(r.size()==1) { x = x_coeff.back(); y = y_coeff.back(); } else { x = *(x_coeff.end()-2); y = *(y_coeff.end()-2); } std::vector&lt;int&gt;::iterator it; std::cout &lt;&lt; "r: "; for(it = r.begin(); it != r.end(); it++) { std::cout &lt;&lt; *it &lt;&lt; "," ; } std::cout &lt;&lt; "\nq: "; for(it = q.begin(); it != q.end(); it++) { std::cout &lt;&lt; *it &lt;&lt; "," ; } std::cout &lt;&lt; "\nx: "; for(it = x_coeff.begin(); it != x_coeff.end(); it++){ std::cout &lt;&lt; *it&lt;&lt;",";} std::cout &lt;&lt; "\ny: "; for(it = y_coeff.begin(); it != y_coeff.end(); it++){ std::cout &lt;&lt; *it&lt;&lt;",";} return a; } </code></pre> <p>by passing to it an <a href="http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm" rel="nofollow noreferrer">example</a> from wikipedia for <code>a = 120</code> and <code>b = 23</code> we obtain:</p> <pre><code>int main(int argc, char** argv) { // 120x + 23y = gcd(120,23) int x_solution, y_solution; int greatestCommonDivisor = gcd_iterative(120, 23, x_solution, y_solution); return 0; } </code></pre> <p>r: 5,3,2,1,0,</p> <p>q: 5,4,1,1,2,</p> <p>x: 1,0,1,-4,5,-9,23,</p> <p>y: 0,1,-5,21,-26,47,-120,</p> <p>what is in accordance with the given table for this example:</p> <p><img src="https://i.stack.imgur.com/HVqy0.png" alt="enter image description here"></p>