Une équation diophantienne linéaire à deux variables est une équation de la forme ax + par = c , où a , b et c sont des entiers constants et x et y sont des variables entières.

Pour de nombreuses équations diophantiennes naturelles, x et y représentent des quantités qui ne peuvent pas être négatives.

Tâche

Écrivez un programme ou une fonction qui accepte les coefficients a , b et c en entrée et renvoie une paire arbitraire de nombres naturels (0, 1, 2,…) x et y qui vérifient l'équation ax + by = c , si une telle paire existe.

Règles supplémentaires

• Vous pouvez choisir n'importe quel format d'entrée et de sortie qui n'implique que les entiers souhaités et, éventuellement, la notation tableau / liste / matrice / tuple / vecteur de votre langue, tant que vous n'incorporez aucun code dans l'entrée.

• Vous pouvez supposer que les coefficients a et b sont tous deux non nuls.

• Votre code doit fonctionner pour tout triplet d'entiers compris entre -2 ⁶⁰ et 2 ⁶⁰ ; il doit se terminer en moins d'une minute sur ma machine (Intel i7-3770, 16 Go de RAM).

• Vous ne pouvez pas utiliser de modules intégrés qui résolvent les équations diophantiennes et banalisent ainsi cette tâche, comme Mathematica FindInstanceou FrobeniusSolve.

• Votre code peut se comporter comme vous le souhaitez si aucune solution ne peut être trouvée, tant qu'il respecte le délai et que sa sortie ne peut pas être confondue avec une solution valide.

• Les règles de code-golf standard s'appliquent.

Exemples

1. Les exemples ci-dessous illustrent des E / S valides pour l'équation 2x + 3y = 11 , qui a exactement deux solutions valides ( (x, y) = (4,1) et (x, y) = (1,3) ).

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

2. La seule solution valable de 2x + 3y = 2 est la paire (x, y) = (1,0) .

3. Les exemples ci-dessous illustrent des E / S valides pour l'équation 2x + 3y = 1 , qui n'a pas de solutions valides .

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

4. Pour (a, b, c) = (1152921504606846883, -576460752303423433, 1) , toutes les solutions correctes (x, y) satisfont que (x, y) = (135637824071393749 - bn, 271275648142787502 + an) pour un entier non négatif n .

Pyth, 92 octets

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

C'est tout à fait un monstre.

Essayez-le en ligne: démonstration . Le format d'entrée est c\n[a,b]et le format de sortie est [x,y].

Dans le cas où aucune solution entière n'existe, je n'imprimerai rien, et dans le cas où aucune solution entière naturelle n'existe, j'imprimerai simplement une solution entière aléatoire.

Explication (aperçu approximatif)

1. Dans un premier temps, je vais trouver une solution entière à l'équation ax + by = gcd(a,b)en utilisant l'algorithme euclidien étendu.

2. Ensuite, je vais modifier la solution (ma multiplication aet bavec c/gcd(a,b)) pour obtenir une solution entière de ax + by = c. Cela fonctionne, si c/gcd(a,b)est un entier. Sinon, il n'existe pas de solution.

3. Toutes les autres solutions entières ont la forme a(x+n*b/d) + b(y-n*a/d) = c avec d = gcd(a,b)for integer n. En utilisant les deux inégalités x+n*b/d >= 0et y-n*a/d >= 0je peux déterminer 6 valeurs possibles pour n. Je vais essayer les 6 et imprimer la solution avec le coefficient le plus bas le plus élevé.

Explication (détaillée)

La première étape consiste à trouver une solution entière à l'équation ax' + by' = gcd(a,b). Cela peut être fait en utilisant l'algorithme euclidien étendu. Vous pouvez vous faire une idée de son fonctionnement sur Wikipedia . La seule différence est qu'au lieu d'utiliser 3 colonnes ( r_i s_i t_i), j'utiliserai 6 colonnes ( r_i-1 r_i s_i-1 s_i t_i-1 t_i). De cette façon, je n'ai pas à garder les deux dernières lignes en mémoire, seulement la dernière.

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. 

Maintenant, je veux trouver une solution ax + by = c. Cela n'est possible que lorsque c mod gcd(a,b) == 0. Si cette équation est satisfaite, je multiplie simplement x',y'avec 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]

Nous avons une solution entière pour ax + by = c. Notez que x, you les deux peuvent être négatifs. Notre objectif est donc de les transformer en non-négatifs.

La bonne chose à propos des équations diophantiennes est que nous pouvons décrire toutes les solutions en utilisant une seule solution initiale. Si (x,y)est une solution, que toutes les autres solutions sont de la forme (x-n*b/gcd(a,b),y+n*a/gcd(a,b))pour nentier.

Par conséquent, nous voulons trouver un n, où x-n*b/gcd(a,b) >= 0et y+n*a/gcd(a,b >= 0. Après une certaine transformation, nous nous retrouvons avec les deux inégalités n >= -x*gcd(a,b)/bet n >= y*gcd(a,b)/a. Notez que le symbole d'inégalité peut regarder dans l'autre sens en raison de la division avec un potentiel négatif aou b. Je m'en fiche pas mal, je dis simplement qu'un nombre de -x*gcd(a,b)/b - 1, -x*gcd(a,b)/b, -x*gcd(a,b)/b + 1satisfait définitivement l'inégalité 1, et un nombre de y*gcd(a,b)/a - 1, y*gcd(a,b)/a, y*gcd(a,b)/a + 1satisfait l'inégalité 2. Il y a unn , qui satisfait les deux inégalités, l'un des 6 chiffres aussi.

Ensuite, je calcule les nouvelles solutions (x-n*b/gcd(a,b),y+n*a/gcd(a,b))pour les 6 valeurs possibles de n. Et j'imprime la solution avec la valeur la plus basse la plus élevée.

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

Le tri par leur ordre de tri fonctionne de la manière suivante. J'utilise l'exemple2x + 3y = 11

Je trie chacune des 6 solutions (ce sont les clés) et trie les solutions originales par leurs clés:

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]

Cela trie une solution non négative complète à la fin (le cas échéant).

• après les remarques de Dennis, qui ont bouleversé mon idée précédente, j'ai dû changer le code de ses racines et cela m'a pris un débogage à long terme, et m'a coûté deux fois le n ° d'octets: '(.

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

• Eh bien, je sais que ce n'est pas un golf, car ce type de langages n'est pas adapté à la réduction de la longueur du code, mais je peux m'assurer que la complexité temporelle est à son meilleur.

Explication:

• le code prend en entrée trois invariants a, b, c, ces derniers sont soumis à quelques conditions avant de procéder au calcul:

  1- si (a + b> c) et (a, b, c> 0) pas de solution!

  2- si (a + b <c), (a, b, c <0) pas de solution!

  3- si (a, b) ont des signes opposés communs de c: pas de solution!

  4- Si GCD (a, b) ne divise pas c, alors plus de solution! sinon, divisez toutes les variantes par GCD.

• après cela, nous devons vérifier une autre condition, elle devrait faciliter et raccourcir le chemin vers la solution souhaitée.

  5- si c divise a ou b, solution s = (x ou y) = (c- [ax, yb]) / [b, a] = C / [b, a] + [ax, yb] / [b , a] = S + [ax, yb] / [b, a] où S est naturel donc ax / b ou by / a aurait désormais des solutions directes non négatives qui sont respectivement x = b ou y = a. (notez que les solutions peuvent être des valeurs nulles au cas où des solutions arbitraires précédentes seraient révélées négatives)

• lorsque le programme atteint ce stade, une gamme plus étroite de solutions pour x = (c-yb) / a est balayée à la place, grâce à la congruence, de balayer de plus grandes gammes de nombres, qui se répètent de manière répétitive par des cycles réguliers. le plus grand champ de recherche est [xa, x + a] où a est le diviseur.

ESSAYEZ-LE

Axiome, 460 octets

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 et un 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

Dans les autres «solutions» possibles, il y avait un bogue car il tentait de sauvegarder les solutions infinies dans une liste; maintenant, il est imposé la limite de 80 solutions max