Naturally linear Diophantine equations A linear Diophantine equation in two variables is an equation of the form ax + by = c, where a, b and c are constant integers and x and y are integer variables.

For many naturally occurring Diophantine equations, x and y represent quantities that cannot be negative.

Task

Write a program or function that accepts the coefficients a, b and c as input and returns an arbitrary pair of natural numbers (0, 1, 2, &hellip;) x and y that verify the equation ax + by = c, if such a pair exists.

Additional rules


You can choose any format for input and output that involves only the desired integers and, optionally, array/list/matrix/tuple/vector notation of your language, as long as you don't embed any code in the input.
You may assume that the coefficients a and b are both non-zero.
Your code must work for any triplet of integers between -260 and 260; it must finish in under a minute on my machine (Intel i7-3770, 16 GiB RAM).
You may not use any built-ins that solve Diophantine equations and thus trivialize this task, such as Mathematica's FindInstance or FrobeniusSolve.
Your code may behave however you want if no solution can be found, as long as it complies with the time limit and its output cannot be confused with a valid solution.
Standard code-golf rules apply.


Examples


The examples below illustrate valid I/O for the equation 2x + 3y = 11, which has exactly two valid solutions ( (x,y) = (4,1) and (x,y) = (1,3) ).

Input:  2 3 11
Output: [4 1]

Input:  (11 (2,3))
Output: [3],(1)

The only valid solution of 2x + 3y = 2 is the pair (x,y) = (1,0).
The examples below illustrate valid I/O for the equation 2x + 3y = 1, which has no valid solutions.

Input:  (2 3 1)
Output: []

Input:  1 2 3
Output: -1

Input:  [[2], [3], [1]]
Output: (2, -1)

For (a, b, c) = (1152921504606846883, -576460752303423433, 1), all correct solutions (x,y) satisfy that (x,y) = (135637824071393749 - bn, 271275648142787502 + an) for some non-negative integer n.
Pyth, 92 bytes

I!%vzhK%2u?sm,ed-hd*ed/F&lt;G2cG2@G1G+~Q,hQ_eQj9 2)J*L/vzhKtKeoSNm-VJ/RhK_*LdQsm+LdtM3/V*LhK_JQ


It's quite a monster. 

Try it online: Demonstration. The input format is c\n[a,b] and the output format is [x,y]. 

In the case that no integer solution exists, I'll print nothing, and in the case that no natural integer solution exists, I'll simply print a random integer solution. 

Explanation (Rough Overview)


At first I'll find an integer solution to the equation ax + by = gcd(a,b) by using the Extended Euclidean algorithm. 
Then I'll modify the solution (my multiplying a and b with c/gcd(a,b)) to get an integer solution of ax + by = c. This works, if c/gcd(a,b) is an integer. Otherwise there doesn't exist a solution. 
All the other integer solutions have the form a(x+n*b/d) + b(y-n*a/d) = c  with d = gcd(a,b) for integer n. Using the two inequalities x+n*b/d &gt;= 0 and y-n*a/d &gt;= 0 I can determine 6 possible values for n. I'll try all 6 of them and print the solution with the highest lowest coefficient. 


Explanation (Detailed)

The first step is to find an integer solution to the equation ax' + by' = gcd(a,b). This can be done by using the extended Euclidean algorithm. You can get an idea on how it works at Wikipedia. The only difference is, that instead of using 3 columns (r_i s_i t_i) I'll use 6 columns (r_i-1 r_i s_i-1 s_i t_i-1 t_i). This way I don't have to keep the last two rows in memory, only the last one. 

K%2u?sm,ed-hd*ed/F&lt;G2cG2@G1G+~Q,hQ_eQj9 2)   implicit: Q = [a,b] (from input)
                                     j9 2    convert 9 to base 2: [1,0,0,1]
                            + Q              add to Q =&gt; [a,b,1,0,0,1]
                                             this is the initial row
   u                                     )   start with G = ^ and update G repeatedly
                                             by the following expression, until
                                             the value of G doesn't change anymore
    ?                   @G1                    if G[1] != 0:
                     cG2                         split G into parts of 2
      m                                          map the parts d to:
       ,                                           the pair 
        ed                                           d[1]
          -hd*ed/F&lt;G2                                d[0]-d[1]*G[0]/G[1]
     s                                           unfold
                                               else:
                           G                     G (don't change it, stop criterion for u)
 %2                                          take every second element
                                             we get the list [gcd(a,b),x',y']
K                                            store this list in K
                             ~Q,hQ_eQ        afterwards change Q to [Q[0],-Q[1]] = [a,-b]
                                             This will be important for the other parts. 


Now I want to find a solution to ax + by = c. This is possible only, when c mod gcd(a,b) == 0. If this equation is satisfied, I simply multiplying x',y' with c/gcd(a,b). 

I!%vzhK...J*L/vzhKtK   implicit: z = c in string format (from input)
  %vzhK                evaluated(z) mod K[0] (=gcd(a,b))
I!                     if not ^ than: 
             /vzhK        c/K[0]
           *L     tK      multipy ^ to each element in K[1:] (=[x',y'])
          J               and store the result in J, this is now [x,y]


We have an integer solution for ax + by = c. Notice, that x, y or both may be negative. So our goal is to transform these to non-negative. 

The nice thing about Diophantine equations is, that we can describe all solution using only one initial solution. If (x,y) is a solution, that all other solutions are of the form (x-n*b/gcd(a,b),y+n*a/gcd(a,b)) for n integer. 

Therefore we want to find a n, where x-n*b/gcd(a,b) &gt;= 0 and y+n*a/gcd(a,b &gt;= 0. After some transformation we end up with the two inequalities n &gt;= -x*gcd(a,b)/b and n &gt;= y*gcd(a,b)/a. Notice that the inequality symbol might look in the other direction due the division with a potential negative a or b. I don't care that much about it, I simply say that one number of -x*gcd(a,b)/b - 1, -x*gcd(a,b)/b, -x*gcd(a,b)/b + 1 definitly satisfies inequality 1, and one number of y*gcd(a,b)/a - 1, y*gcd(a,b)/a, y*gcd(a,b)/a + 1 satisfies inequality 2. It there is a n, that satisfies both inequalities, one of the 6 numbers also does. 

Then I calculate the new solutions (x-n*b/gcd(a,b),y+n*a/gcd(a,b)) for all 6 possible values of n. And I print the solution with the highest lowest value. 

eoSNm-VJ/RhK_*LdQsm+LdtM3/V*LhK_JQ
                               _J    reverse J =&gt; [y,x]
                           *LhK      multiply each value with K[0] =&gt; [y*gcd,x*gcd]
                         /V      Q   vectorized division =&gt; [y*gcd/a,-x*gcd/b]
                  m                  map each d of ^ to:
                      tM3              [-1,0,1]
                   +Ld                 add d to each ^
                 s                   unfold
                                     these are the possible values for n
    m                                map each d (actually n) of ^ to:
             *LdQ                      multiply d to Q =&gt; [a*n,-b*n]
            _                          reverse =&gt; [-b*n,a*n]
        /RhK                           divide by K[0] =&gt; [-b*n/gcd,a*n/gcd]
     -VJ                               vectorized subtraction with J
                                       =&gt; [x+b*n/gcd,y-a*n/gcd]
 oSN                                 order the solutions by their sorted order
e                                    print the last one


The sort by their sorted order thing works the following way. I'm using the example 2x + 3y = 11

I sort each of the 6 solutions (this are called keys), and sort the original solutions by their keys:

solutions: [1, 3], [4, 1], [7, -1], [-5, 7], [-2, 5], [1, 3]
keys:      [1, 3], [1, 4], [-1, 7], [-5, 7], [-2, 5], [1, 3]
sort by key:
solutions: [-5, 7], [-2, 5], [7, -1], [1, 3], [1, 3], [4, 1]
keys:      [-5, 7], [-2, 5], [-1, 7], [1, 3], [1, 3], [1, 4]


This sorts a complete non-negative solution to the end (if there is any).