Um linear equação Diofantina em duas variáveis é uma equação da forma ax + by = C , em que um , b e c são números inteiros constantes e x e y são inteiros variáveis.

Para muitas equações diofantinas de ocorrência natural, x e y representam quantidades que não podem ser negativas.

Tarefa

Escrever um programa ou função que aceita os coeficientes de um , b e c como entrada e devolve um par arbitrário de números naturais (0, 1, 2, ...) x e y que verificar a equação ax + by = c , se um tal par existe.

Regras adicionais

• Você pode escolher qualquer formato de entrada e saída que envolva apenas os números inteiros desejados e, opcionalmente, notação de matriz / lista / matriz / tupla / vetor do seu idioma, desde que não incorpore nenhum código à entrada.

• Você pode assumir que os coeficientes de um e b são ambos não-zero.

• Seu código deve funcionar para qualquer trio de números inteiros entre -2 ⁶⁰ e 2 ⁶⁰ ; ele deve terminar em menos de um minuto na minha máquina (Intel i7-3770, 16 GiB RAM).

• Você não pode usar nenhum built-in que resolva equações diofantinas e, portanto, banalize esta tarefa, como as do Mathematica FindInstanceou FrobeniusSolve.

• Seu código pode se comportar como quiser, se nenhuma solução puder ser encontrada, desde que cumpra o prazo e sua saída não possa ser confundida com uma solução válida.

• Aplicam-se regras padrão de código de golfe .

Exemplos

1. Os exemplos abaixo ilustram E / S válidas para a equação 2x + 3y = 11 , que possui exatamente duas soluções válidas ( (x, y) = (4,1) e (x, y) = (1,3) ).

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

2. A única solução válida de 2x + 3y = 2 é o par (x, y) = (1,0) .

3. Os exemplos abaixo ilustram E / S válidas para a equação 2x + 3y = 1 , que não possui soluções 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 as soluções corretas (x, y) satisfazem que (x, y) = (135637824071393749 - bn, 271275648142787502 + an) para algum número inteiro não negativo n .

Pitão, 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

É um monstro e tanto.

Experimente online: Demonstração . O formato de entrada é c\n[a,b]e o formato de saída é [x,y].

Caso não exista uma solução inteira, não imprimirei nada e, caso não exista uma solução natural natural, simplesmente imprimirei uma solução inteira aleatória.

Explicação (Visão Geral Bruta)

1. Inicialmente, encontrarei uma solução inteira para a equação ax + by = gcd(a,b)usando o algoritmo Euclidiano Estendido.

2. Então eu modificarei a solução (minha multiplicação ae bcom c/gcd(a,b)) para obter uma solução inteira de ax + by = c. Isso funciona, se c/gcd(a,b)é um número inteiro. Caso contrário, não existe uma solução.

3. Todas as outras soluções inteiras têm o formato a(x+n*b/d) + b(y-n*a/d) = c com d = gcd(a,b)for inteiro n. Usando as duas desigualdades x+n*b/d >= 0e y-n*a/d >= 0eu posso determinar 6 valores possíveis para n. Vou tentar todos os 6 deles e imprimir a solução com o maior coeficiente mais baixo.

Explicação (detalhada)

O primeiro passo é encontrar uma solução inteira para a equação ax' + by' = gcd(a,b). Isso pode ser feito usando o algoritmo euclidiano estendido. Você pode ter uma idéia de como isso funciona na Wikipedia . A única diferença é que, em vez de usar 3 colunas ( r_i s_i t_i), usarei 6 colunas ( r_i-1 r_i s_i-1 s_i t_i-1 t_i). Dessa forma, não preciso manter as duas últimas linhas na memória, apenas a ú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. 

Agora eu quero encontrar uma solução para ax + by = c. Isso é possível apenas quando c mod gcd(a,b) == 0. Se esta equação for satisfeita, eu simplesmente multiplico x',y'com 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]

Temos uma solução inteira para ax + by = c. Aviso, que x, you ambos, podem ser negativo. Portanto, nosso objetivo é transformá-los em não negativos.

O bom das equações diofantinas é que podemos descrever todas as soluções usando apenas uma solução inicial. Se (x,y)é uma solução, todas as outras soluções têm a forma (x-n*b/gcd(a,b),y+n*a/gcd(a,b))de nnúmero inteiro.

Portanto, queremos encontrar a n, onde x-n*b/gcd(a,b) >= 0e y+n*a/gcd(a,b >= 0. Após alguma transformação, terminamos com as duas desigualdades n >= -x*gcd(a,b)/be n >= y*gcd(a,b)/a. Observe que o símbolo de desigualdade pode olhar na outra direção devido à divisão com um potencial negativo aou b. Não me importo muito com isso, simplesmente digo que um número de -x*gcd(a,b)/b - 1, -x*gcd(a,b)/b, -x*gcd(a,b)/b + 1satisfaz definitivamente a desigualdade 1 e um número de y*gcd(a,b)/a - 1, y*gcd(a,b)/a, y*gcd(a,b)/a + 1satisfaz a desigualdade 2. Existe um nque satisfaz ambas as desigualdades, um dos 6 números também satisfaz.

Então eu calculo as novas soluções (x-n*b/gcd(a,b),y+n*a/gcd(a,b))para todos os 6 valores possíveis de n. E imprimo a solução com o menor valor mais 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

A classificação por ordem de classificação funciona da seguinte maneira. Estou usando o exemplo2x + 3y = 11

Classifico cada uma das 6 soluções (chamadas de chaves) e classifico as soluções originais por suas chaves:

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]

Isso classifica uma solução completa não negativa até o fim (se houver).

• Após as observações de Dennis, que fizeram minha idéia anterior de cabeça para baixo, eu tive que mudar o código de suas raízes e isso levou a depuração de longo prazo e me custou duas vezes 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

• Bem, eu sei que não é um jogador de golfe, pois esse tipo de linguagem não é adaptado para reduzir o tamanho do código, mas posso garantir que a complexidade do tempo seja a melhor possível.

Explicação:

• o código usa três invariantes a, b, c como entrada, esses últimos são subdivididos em duas condições antes de proceder ao cálculo:

  1- se (a + b> c) e (a, b, c> 0) não há solução!

  2- se (a + b <c), (a, b, c <0) não há solução!

  3- se (a, b) têm sinais opostos comuns de c: sem solução!

  4- se o GCD (a, b) não dividir c, não haverá solução novamente! caso contrário, divida todas as variantes por GCD.

• depois disso, temos que verificar outra condição, deve facilitar e facilitar o caminho para a solução desejada.

  5- se c divide a ou b, solução s = (x ou y) = (c- [ax, yb]) / [b, a] = C / [b, a] + [ax, yb] / [b , a] = S + [ax, yb] / [b, a] onde S é natural, então ax / b ou by / a teriam doravante soluções diretas não-negativas que são respectivamente x = b ou y = a. (observe que as soluções podem ser apenas valores nulos, caso soluções arbitrárias anteriores sejam reveladas negativas)

• quando o programa chega a esse estágio, uma gama mais estreita de soluções para x = (c-yb) / a é varrida, graças à congruência de varrer faixas maiores de números, que são repetidas repetidamente por ciclos regulares. o maior campo de pesquisa é [xa, x + a] onde a é o divisor.

TENTE

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 e algum teste

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

Nas outras 'soluções' possíveis, houve um erro porque tentou salvar as infinitas soluções em uma lista; agora é imposto o limite de 80 soluções no máximo