A lineal ecuación Diophantine en dos variables es una ecuación de la forma ax + by = c , donde un , b y c son números enteros constantes y x y y son enteros variables.

Para muchas ecuaciones naturales de Diophantine, x e y representan cantidades que no pueden ser negativas.

Tarea

Escribir un programa o función que acepta los coeficientes de un , b y c como entrada y devuelve un par arbitrario de números naturales (0, 1, 2, ...) x y y que verifican la ecuación ax + by = c , si un par tal existe

Reglas adicionales

• Puede elegir cualquier formato para entrada y salida que involucre solo los enteros deseados y, opcionalmente, la notación de matriz / lista / matriz / tupla / vector de su idioma, siempre que no incruste ningún código en la entrada.

• Puede suponer que los coeficientes a y b son ambos distintos de cero.

• Su código debe funcionar para cualquier triplete de enteros entre -2 ⁶⁰ y 2 ⁶⁰ ; debe terminar en menos de un minuto en mi máquina (Intel i7-3770, 16 GiB RAM).

• No puede utilizar ningún elemento integrado que resuelva las ecuaciones de Diophantine y, por lo tanto, trivialice esta tarea, como las de Mathematica FindInstanceo FrobeniusSolve.

• Su código puede comportarse como quiera si no se encuentra una solución, siempre que cumpla con el límite de tiempo y su salida no se pueda confundir con una solución válida.

• Se aplican reglas estándar de código de golf .

Ejemplos

1. Los siguientes ejemplos ilustran E / S válidas para la ecuación 2x + 3y = 11 , que tiene exactamente dos soluciones válidas ( (x, y) = (4,1) y (x, y) = (1,3) ).

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

2. La única solución válida de 2x + 3y = 2 es el par (x, y) = (1,0) .

3. Los siguientes ejemplos ilustran E / S válidas para la ecuación 2x + 3y = 1 , que no tiene soluciones válidas .

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

4. Para (a, b, c) = (1152921504606846883, -576460752303423433, 1) , todas las soluciones correctas (x, y) satisfacen que (x, y) = (135637824071393749 - bn, 271275648142787502 + an) para algún número entero no negativo 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 todo un monstruo.

Pruébelo en línea: demostración . El formato de entrada es c\n[a,b]y el formato de salida es [x,y].

En el caso de que no exista una solución entera, no imprimiré nada, y en el caso de que no exista una solución entera natural, simplemente imprimiré una solución entera aleatoria.

Explicación (descripción general)

1. Al principio, encontraré una solución entera para la ecuación ax + by = gcd(a,b)usando el algoritmo Euclidiano Extendido.

2. Luego modificaré la solución (mi multiplicación ay bcon c/gcd(a,b)) para obtener una solución entera de ax + by = c. Esto funciona, si c/gcd(a,b)es un número entero. De lo contrario, no existe una solución.

3. Todas las otras soluciones de enteros tienen la forma a(x+n*b/d) + b(y-n*a/d) = c con d = gcd(a,b)for integer n. Usando las dos desigualdades x+n*b/d >= 0y y-n*a/d >= 0puedo determinar 6 valores posibles para n. Probaré los 6 e imprimiré la solución con el coeficiente más bajo más alto.

Explicación (detallada)

El primer paso es encontrar una solución entera a la ecuación ax' + by' = gcd(a,b). Esto se puede hacer usando el algoritmo euclidiano extendido. Puede hacerse una idea de cómo funciona en Wikipedia . La única diferencia es que, en lugar de usar 3 columnas ( r_i s_i t_i), usaré 6 columnas ( r_i-1 r_i s_i-1 s_i t_i-1 t_i). De esta manera no tengo que mantener las dos últimas filas en la memoria, solo la última.

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. 

Ahora quiero encontrar una solución para ax + by = c. Esto solo es posible cuando c mod gcd(a,b) == 0. Si se satisface esta ecuación, simplemente multiplico x',y'con 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]

Tenemos una solución entera para ax + by = c. Aviso, que x, yo ambos puede ser negativo. Por lo tanto, nuestro objetivo es transformarlos en no negativos.

Lo bueno de las ecuaciones de diofantina es que podemos describir todas las soluciones usando solo una solución inicial. Si (x,y)es una solución, todas las demás soluciones tienen la forma (x-n*b/gcd(a,b),y+n*a/gcd(a,b))de nentero.

Por lo tanto, queremos encontrar un n, donde x-n*b/gcd(a,b) >= 0y y+n*a/gcd(a,b >= 0. Después de alguna transformación, terminamos con las dos desigualdades n >= -x*gcd(a,b)/by n >= y*gcd(a,b)/a. Observe que el símbolo de desigualdad podría mirar en la otra dirección debido a la división con un potencial negativo ao b. No me importa mucho, simplemente digo que un número de -x*gcd(a,b)/b - 1, -x*gcd(a,b)/b, -x*gcd(a,b)/b + 1satisface definitivamente la desigualdad 1, y un número de y*gcd(a,b)/a - 1, y*gcd(a,b)/a, y*gcd(a,b)/a + 1satisface la desigualdad 2. Si hay un n, que satisface ambas desigualdades, uno de los 6 números también lo hace.

Luego calculo las nuevas soluciones (x-n*b/gcd(a,b),y+n*a/gcd(a,b))para los 6 valores posibles de n. E imprimo la solución con el valor más bajo más alto.

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

La ordenación por su orden ordenado funciona de la siguiente manera. Estoy usando el ejemplo2x + 3y = 11

Ordeno cada una de las 6 soluciones (esto se llama claves), y clasifico las soluciones originales por sus claves:

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]

Esto ordena una solución no negativa completa hasta el final (si hay alguna).

• Después de los comentarios de Dennis, que hicieron que mi idea anterior se volcara, tuve que cambiar el código desde sus raíces y me llevó una depuración a largo plazo, y me costó dos veces n ° de 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

• Bueno, sé que no es golf, ya que ese tipo de idiomas no está adaptado para la reducción de la longitud del código, pero puedo asegurar que la complejidad del tiempo sea óptima.

Explicación:

• el código toma tres invariantes a, b, c como entrada, estos últimos se someten a un par de condiciones antes de proceder a calcular:

  1- si (a + b> c) y (a, b, c> 0) no hay solución!

  2- si (a + b <c), (a, b, c <0) no hay solución!

  3- si (a, b) tienen signos opuestos comunes de c: ¡no hay solución!

  4- si el dosificador GCD (a, b) divide c, ¡entonces no hay solución otra vez! , de lo contrario, divida todas las variantes por GCD.

• después de esto, tenemos que verificar otra condición, debería facilitar y acortar el camino hacia la solución deseada.

  5- si c divide a o b, solución s = (x o y) = (c- [ax, yb]) / [b, a] = C / [b, a] + [ax, yb] / [b , a] = S + [ax, yb] / [b, a] donde S es natural, por lo que ax / b o by / a tendrían soluciones directas no negativas que son respectivamente x = b o y = a. (observe que las soluciones pueden ser solo valores nulos en caso de que las soluciones arbitrarias anteriores se revelen negativas)

• cuando el programa alcanza esta etapa, se barre un rango más estrecho de soluciones para x = (c-yb) / a, gracias a la congruencia, de barrer rangos más grandes de números, que se encuentran repetidamente en ciclos regulares. el campo de búsqueda más grande es [xa, x + a] donde a es el divisor.

INTENTALO

Axioma, 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 y alguna prueba

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

En las otras 'soluciones' posibles había un error porque intentaba guardar las soluciones infinitas en una Lista; ahora se impone el límite de 80 soluciones máximo