Eine lineare diophantinische Gleichung in zwei Variablen ist eine Gleichung der Form ax + by = c , wobei a , b und c konstante ganze Zahlen und x und y ganzzahlige Variablen sind.

Für viele natürlich vorkommende diophantinische Gleichungen stehen x und y für Größen, die nicht negativ sein können.

Aufgabe

Schreiben Sie ein Programm oder eine Funktion, die die Koeffizienten a , b und c als Eingabe akzeptiert und ein beliebiges Paar natürlicher Zahlen (0, 1, 2,…) x und y zurückgibt , die die Gleichung ax + by = c verifizieren , falls ein solches Paar vorliegt existiert.

Zusätzliche Regeln

• Sie können ein beliebiges Format für die Eingabe und Ausgabe auswählen, das nur die gewünschten Ganzzahlen und optional Array / Liste / Matrix / Tupel / Vektor Ihrer Sprache enthält, sofern Sie keinen Code in die Eingabe einbetten.

• Sie können annehmen, dass die Koeffizienten a und b beide ungleich Null sind.

• Ihr Code muss für jedes Triplett von ganzen Zahlen zwischen -2 ⁶⁰ und 2 ⁶⁰ funktionieren . Auf meinem Computer (Intel i7-3770, 16 GiB RAM) muss es in weniger als einer Minute fertig sein.

• Sie dürfen keine eingebauten Elemente verwenden, die diophantische Gleichungen lösen und diese Aufgabe damit trivialisieren, wie z. B. die von Mathema- tica FindInstanceoder FrobeniusSolve.

• Ihr Code kann sich beliebig verhalten, wenn keine Lösung gefunden werden kann, solange er das Zeitlimit einhält und seine Ausgabe nicht mit einer gültigen Lösung verwechselt werden kann.

• Es gelten die Standardregeln für Code-Golf .

Beispiele

1. Die folgenden Beispiele veranschaulichen gültige E / A für die Gleichung 2x + 3y = 11 , die genau zwei gültige Lösungen enthält ( (x, y) = (4,1) und (x, y) = (1,3) ).

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

2. Die einzig gültige Lösung von 2x + 3y = 2 ist das Paar (x, y) = (1,0) .

3. Die folgenden Beispiele veranschaulichen gültige E / A für die Gleichung 2x + 3y = 1 , für die keine gültigen Lösungen vorliegen .

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

4. Für (a, b, c) = (1152921504606846883, -576460752303423433, 1) erfüllen alle korrekten Lösungen (x, y) (x, y) = (135637824071393749 - bn, 271275648142787502 + an) für einige nicht negative ganze Zahlen n .

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

Es ist ein ziemliches Monster.

Probieren Sie es online aus: Demonstration . Das Eingabeformat ist c\n[a,b]und das Ausgabeformat ist [x,y].

Für den Fall, dass keine ganzzahlige Lösung existiert, gebe ich nichts aus und für den Fall, dass keine natürliche ganzzahlige Lösung existiert, gebe ich einfach eine zufällige ganzzahlige Lösung aus.

Erklärung (grobe Übersicht)

1. Zuerst werde ich eine ganzzahlige Lösung für die Gleichung finden, ax + by = gcd(a,b)indem ich den erweiterten euklidischen Algorithmus verwende.

2. Dann ändere ich die Lösung (mein Multiplizieren aund bmit c/gcd(a,b)), um eine ganzzahlige Lösung von zu erhalten ax + by = c. Dies funktioniert, wenn c/gcd(a,b)es sich um eine Ganzzahl handelt. Sonst gibt es keine Lösung.

3. Alle anderen Ganzzahllösungen haben die Form a(x+n*b/d) + b(y-n*a/d) = c mit d = gcd(a,b)für Ganzzahl n. Mit den beiden Ungleichungen x+n*b/d >= 0und y-n*a/d >= 0kann ich 6 mögliche Werte für bestimmen n. Ich probiere alle 6 aus und drucke die Lösung mit dem höchsten und niedrigsten Koeffizienten.

Erklärung (detailliert)

Der erste Schritt besteht darin, eine ganzzahlige Lösung für die Gleichung zu finden ax' + by' = gcd(a,b). Dies kann mithilfe des erweiterten euklidischen Algorithmus erfolgen. Bei Wikipedia können Sie sich ein Bild davon machen, wie es funktioniert . Der einzige Unterschied ist, dass r_i s_i t_iich anstelle von 3 Spalten ( ) 6 Spalten ( r_i-1 r_i s_i-1 s_i t_i-1 t_i) verwende. Auf diese Weise muss ich nicht die letzten beiden Zeilen im Speicher behalten, sondern nur die letzte.

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. 

Jetzt möchte ich eine Lösung finden ax + by = c. Dies ist nur möglich, wenn c mod gcd(a,b) == 0. Wenn diese Gleichung erfüllt ist, multipliziere ich einfach x',y'mit 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]

Wir haben eine ganzzahlige Lösung für ax + by = c. Beachten Sie , dass x, yoder beide negativ sein kann. Unser Ziel ist es daher, diese in nicht negative umzuwandeln.

Das Schöne an Diophantine-Gleichungen ist, dass wir alle Lösungen mit nur einer Anfangslösung beschreiben können. Wenn (x,y)es sich um eine Lösung handelt, haben alle anderen Lösungen die Form (x-n*b/gcd(a,b),y+n*a/gcd(a,b))einer nGanzzahl.

Deshalb wollen wir ein n, wo x-n*b/gcd(a,b) >= 0und y+n*a/gcd(a,b >= 0. Nach einer gewissen Transformation kommen wir zu den beiden Ungleichungen n >= -x*gcd(a,b)/bund n >= y*gcd(a,b)/a. Beachten Sie, dass das Ungleichungssymbol aufgrund der Division mit einem potenziellen negativen aoder in die andere Richtung zeigen kann b. Es ist mir egal, ich sage einfach, dass eine Zahl -x*gcd(a,b)/b - 1, -x*gcd(a,b)/b, -x*gcd(a,b)/b + 1definitiv die Ungleichung 1 und eine Zahl y*gcd(a,b)/a - 1, y*gcd(a,b)/a, y*gcd(a,b)/a + 1die Ungleichung 2 erfüllt. Es gibt eine n, die beide Ungleichungen erfüllt, und eine der 6 Zahlen tut dies auch.

Dann berechne ich die neuen Lösungen (x-n*b/gcd(a,b),y+n*a/gcd(a,b))für alle 6 möglichen Werte von n. Und ich drucke die Lösung mit dem höchsten niedrigsten Wert.

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

Das Sortieren nach ihrer sortierten Reihenfolge funktioniert folgendermaßen. Ich benutze das Beispiel2x + 3y = 11

Ich sortiere jede der 6 Lösungen (diese werden als Schlüssel bezeichnet) und sortiere die ursprünglichen Lösungen nach ihren Schlüsseln:

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]

Dies sortiert eine vollständige nicht negative Lösung bis zum Ende (falls vorhanden).

• Nach Dennis 'Bemerkungen, die meine vorherige Idee auf den Kopf stellten, musste ich den Code von Grund auf ändern, und es dauerte lange, bis ich das Debugging durchführte.

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

• Nun, ich weiß, es ist kein Golfsport, da diese Art von Sprachen nicht für die Reduzierung der Codelänge geeignet ist, aber ich kann sicherstellen, dass die zeitliche Komplexität optimal ist.

Erläuterung:

• Der Code verwendet drei Invarianten a, b, c als Eingabe. Diese letzten werden einigen Bedingungen unterworfen, bevor mit der Berechnung fortgefahren wird:

  1- wenn (a + b> c) und (a, b, c> 0) keine Lösung!

  2- wenn (a + b <c), (a, b, c <0) keine Lösung!

  3- wenn (a, b) gemeinsame entgegengesetzte Vorzeichen von c haben: keine Lösung!

  4- Wenn GCD (a, b) c nicht teilt, dann keine Lösung mehr! Ansonsten teilen Sie alle Varianten durch GCD.

• Danach müssen wir eine andere Bedingung überprüfen, sie sollte den Weg zur gewünschten Lösung erleichtern und verkürzen.

  5- Wenn c a oder b dividiert, ist die Lösung s = (x oder y) = (c- [ax, yb]) / [b, a] = C / [b, a] + [ax, yb] / [b a] = S + [ax, yb] / [b, a], wobei S natürlich ist, so dass ax / b oder durch / a fortan nicht negative direkte Lösungen haben würden, die jeweils x = b oder y = a sind. (Beachten Sie, dass Lösungen nur Nullwerte sein können, wenn bei vorherigen beliebigen Lösungen Negative festgestellt werden.)

• Wenn das Programm diese Stufe erreicht, wird ein engerer Bereich von Lösungen für x = (c-yb) / a gewobbelt, anstatt, dank der Kongruenz, größere Bereiche von Zahlen zu wobbeln, was durch regelmäßige Zyklen wiederholt wird. Das größte Suchfeld ist [xa, x + a], wobei a der Divisor ist.

VERSUCH ES

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 und ein 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

Bei den anderen möglichen 'Lösungen' gab es einen Fehler, weil versucht wurde, die unendlichen Lösungen in einer Liste zu speichern. Jetzt ist die Grenze von 80 Lösungen max auferlegt