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.


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.


  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.

  • \$\begingroup\$ 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. \$\endgroup\$
    – Sp3000
    Jun 20, 2015 at 7:31
  • \$\begingroup\$ intput 1 2 3 has a valid output though... [ 1, 1 ] \$\endgroup\$
    – Jack Ammo
    Jun 20, 2015 at 16:46
  • \$\begingroup\$ @JackAmmo: All examples in the second code block correspond to 2x + 3y = 1. \$\endgroup\$
    – Dennis
    Jun 20, 2015 at 16:48
  • \$\begingroup\$ 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? \$\endgroup\$
    – user58988
    Aug 27, 2017 at 14:06
  • \$\begingroup\$ In such case probably I have to return that not exist solution... \$\endgroup\$
    – user58988
    Aug 27, 2017 at 14:12

3 Answers 3


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

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


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


  • \$\begingroup\$ euuh , large numbers issue , gonna fix that (wonder when lol) \$\endgroup\$
    – Abr001am
    Jun 21, 2015 at 4:06
  • \$\begingroup\$ 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 . \$\endgroup\$
    – Abr001am
    Jun 21, 2015 at 9:53
  • \$\begingroup\$ ohh . i forgot to fix that bugging bug :p thanks for noticing it @Jakube \$\endgroup\$
    – Abr001am
    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
  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];[])
  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>=0 and z>=0 and t*a+z*b=k=>(o:=cons([t,z],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)
                                                  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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.