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, …) 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 rules apply.

Examples

1. 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)

2. The only valid solution of 2x + 3y = 2 is the pair (x,y) = (1,0).

3. 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)

4. 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.

• I think it might be good to place a bit more emphasis on nonnegative integers, and that the second example in fact has no solution. Jun 20, 2015 at 7:31
• intput 1 2 3 has a valid output though... [ 1, 1 ] Jun 20, 2015 at 16:46
• @JackAmmo: All examples in the second code block correspond to 2x + 3y = 1. Jun 20, 2015 at 16:48
• In ax+bx=k it seems to me to understand that the solution has to be x>=0 and y>=0. So who are such x,y>=0 solutions of 38*x+909*y=3?
– user58988
Aug 27, 2017 at 14:06
• In such case probably I have to return that not exist solution...
– user58988
Aug 27, 2017 at 14:12

Pyth, 92 bytes

I!%vzhK%2u?sm,ed-hd*ed/F<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)

1. At first I'll find an integer solution to the equation ax + by = gcd(a,b) by using the Extended Euclidean algorithm.

2. 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.

3. 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 >= 0 and y-n*a/d >= 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<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 => [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<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) >= 0 and y+n*a/gcd(a,b >= 0. After some transformation we end up with the two inequalities n >= -x*gcd(a,b)/b and n >= 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 => [y,x]
*LhK      multiply each value with K[0] => [y*gcd,x*gcd]
/V      Q   vectorized division => [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 => [a*n,-b*n]
_                          reverse => [-b*n,a*n]
/RhK                           divide by K[0] => [-b*n/gcd,a*n/gcd]
-VJ                               vectorized subtraction with J
=> [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).

• after Dennis' remarks, that made my previous idea upside-down, i had to change the code from its roots and it took me long term debugging, and cost me twice n° of bytes :'(.

Matlab (660)

a=input('');b=input('');c=input('');if((min(a*c,b*c)>c*c)&&a*c>0&&b*c>0)||(a*c<0&&b*c<0),-1,return,end,g=abs(gcd(a,b));c=c/g;a=a/g;b=b/g;if(c~=floor(c)),-1,return,end,if(c/a==floor(c/a)&&c/a>0),e=c/a-b;if(e>0),e,a,return,else,c/a,0,return,end,end,if(c/b==floor(c/b)&&c/b>0),e=c/b-a;if(e>0),b,e,return,else,0,c/b,return,end,end,f=max(abs(a),abs(b));if f==abs(a),f=b;b=a;a=f;g=0.5;end,e=(c-b)/a;f=(c-2*b)/a;if(e<0&&f<e),-1,elseif(e<0&&f>e),for(i=abs(c*a):abs((c+1)*a)),e=(c-i*b);if(mod(e,a)==0)if(g==0.5),i,e/a;else,e/a,i,end,return,end,end,else for(i=1:abs(a)),e=(c-i*b);if(e/a<0),-1,elseif(mod(e,a)==0),if(g==0.5),i,e/a,else,e/a,i,end,return,end,end,end,-1

• Well , i know its not golfed, since that type of languages isnt adapted for code length reduction, but, i can ensure that time-complexity is at its best.

Explanation:

• the code takes three invariants a,b,c as input, these last ones are subdued to couple of conditions before proceeding to calculate:

1- if (a+b>c) and (a,b,c>0) no solution!

2- if (a+b < c) ,(a,b,c<0) no solution!

3- if (a, b) have common opposite signs of c : no solution!

4- if GCD(a,b) dosnt divide c, then no solution again! , otherwise, divide all variants by GCD.

• after this , we have to check another condition out, it should ease and shorteb the way to desired solution.

5- if c divide a or b , solution s= (x or y)=(c-[ax,yb])/[b,a]=C/[b,a]+[ax,yb]/[b,a]=S+[ax,yb]/[b,a] where S is natural so ax/b or by/a would have henceforth non-negative direct solutions which are respectively x=b or y=a . (notice that solutions can be just nil values in case previous arbitrary solutions are revealed negatives)

• when the program reaches this stage, a narrower range of solutions for x=(c-yb)/a is swept instead, thanks to congruence, of sweeping larger ranges of numbers ,which is come accross repetitively by regular cycles. the largest search field is [x-a,x+a] where a is the divisor.

TRY IT

• euuh , large numbers issue , gonna fix that (wonder when lol) Jun 21, 2015 at 4:06
• I think its still minor bug to fix, about large integers, i still dont get why dividing 1152921504606846800.000000/576460752303423420.000000 comes out with natural number 2 ,though this last result is rounded . Jun 21, 2015 at 9:53
• ohh . i forgot to fix that bugging bug :p thanks for noticing it @Jakube Jun 21, 2015 at 17:26

Axiom, 460 bytes

w(a,b,x,u)==(a=0=>[b,x];w(b rem a,a,u,x-u*(b quo a)))
d(a,b,k)==(o:List List INT:=[];a=0 and b=0=>(k=0=>[1,1];[]);a=0=>(k=0=>[[1,0]];k rem b=0=>[1,k quo b];[]);b=0=>(k=0=>[[0,1]];k rem a=0=>[k quo a,1];[]);r:=w(a,b,0,1);q:=k quo r.1;(y,x,u,v):=(q*(r.1-r.2*a)quo b,q*r.2,b quo r.1,a quo r.1);m:=min(80,4+abs(k)quo min(abs(a),abs(b)));l:=y quo v;x:=x+l*u;y:=y-l*v;for n in -m..m repeat(t:=x+n*u;z:=y-n*v;t>=0 and z>=0 and t*a+z*b=k=>(o:=cons([t,z],o)));sort(o))


ungolf and some test

-- input a b and k for equation a*x+b*y=k
-- result one List of List of elments [x,y] of solution of
-- that equation with x and y NNI (not negative integers)
-- or Void list [] for no solution
diopanto(a,b,k)==
o:List List INT:=[]
a=0 and b=0=>(k=0=>[1,1];[])
a=0=>(k=0=>[[1,0]];k rem b=0=>[1,k quo b];[])
b=0=>(k=0=>[[0,1]];k rem a=0=>[k quo a,1];[])
r:=w(a,b,0,1)
q:=k quo r.1
(y,x,u,v):=(q*(r.1-r.2*a)quo b,q*r.2,b quo r.1,a quo r.1)
m:=min(80,4+abs(k)quo min(abs(a),abs(b)))
l:=y quo v           -- center the interval
x:=x+l*u; y:=y-l*v
for n in -m..m repeat
t:=x+n*u;z:=y-n*v
t>=0 and z>=0 and t*a+z*b=k=>(o:=cons([t,z],o))
sort(o)

------------------------------------------------------
(4) -> d(0,-9,0)
(4)  [[1,0]]
Type: List List Integer
(5) -> d(2,3,11)
(5)  [[4,1],[1,3]]
Type: List List Integer
(6) -> d(2,3,2)
(6)  [[1,0]]
Type: List List Integer
(7) -> d(2,3,1)
(7)  []
Type: List List Integer
(8) -> d(1152921504606846883,-576460752303423433,1)
(8)
[[135637824071393749,271275648142787502],
[712098576374817182,1424197152749634385],
[1288559328678240615,2577118657356481268],
[1865020080981664048,3730040161963328151],
[2441480833285087481,4882961666570175034]]
Type: List List Integer


In the other 'solutions' possible there was a bug because it tried to save the infinite solutions in one List; now it is imposed the limit of 80 solutions max