Naturalnie liniowe równania diofantyczne

13

Liniowy równanie diofantycznego dwóch zmiennych jest równanie postaci ax + by = C , gdzie , b oraz c są liczbami całkowitymi, stałe i x i y są liczbami całkowitymi zmiennych.

Dla wielu naturalnie występujących diofantyczne równania, x i y oznaczają ilości, które nie mogą być ujemne.

Zadanie

Napisać program lub funkcję, która przyjmuje współczynniki , b i c na wejściu i zwraca dowolną parę liczb naturalnych (0, 1, 2, ...), x i y to sprawdzenie równania ax + by = C , w razie takiej pary istnieje.

Dodatkowe zasady

  • Możesz wybrać dowolny format wejściowy i wyjściowy, który obejmuje tylko pożądane liczby całkowite i, opcjonalnie, tablicę / listę / macierz / krotkę / wektorową notację twojego języka, o ile nie osadzisz żadnego kodu na wejściu.

  • Można założyć, że współczynniki a i b są zarówno niezerową.

  • Twój kod musi działać dla każdej trojaczki liczb całkowitych od -2 60 do 2 60 ; na moim komputerze musi zakończyć się w niecałą minutę (Intel i7-3770, 16 GiB RAM).

  • Nie możesz używać żadnych wbudowanych rozwiązań, które rozwiązują równania diofantyczne, a tym samym trywializują to zadanie, takie jak Mathematica FindInstancelub FrobeniusSolve.

  • Twój kod może zachowywać się tak, jak chcesz, jeśli nie można znaleźć rozwiązania, o ile jest ono zgodne z limitem czasu, a jego wynik nie może być mylony z prawidłowym rozwiązaniem.

  • Obowiązują standardowe zasady .

Przykłady

  1. Poniższe przykłady ilustrują prawidłowe operacje we / wy dla równania 2x + 3y = 11 , które ma dokładnie dwa prawidłowe rozwiązania ( (x, y) = (4,1) i (x, y) = (1,3) ).

    Input:  2 3 11
    Output: [4 1]
    
    Input:  (11 (2,3))
    Output: [3],(1)
    
  2. Jedynym prawidłowym rozwiązaniem 2x + 3y = 2 jest para (x, y) = (1,0) .

  3. Poniższe przykłady ilustrują prawidłowe operacje we / wy dla równania 2x + 3y = 1 , które nie ma prawidłowych rozwiązań .

    Input:  (2 3 1)
    Output: []
    
    Input:  1 2 3
    Output: -1
    
    Input:  [[2], [3], [1]]
    Output: (2, -1)
    
  4. Dla (a, b, c) = (1152921504606846883, -576460752303423433, 1) wszystkie prawidłowe rozwiązania (x, y) spełniają to (x, y) = (135637824071393749 - bn, 271275648142787502 + an) dla pewnej liczby całkowitej nieujemnej n .

Dennis
źródło
Myślę, że dobrze byłoby położyć nieco większy nacisk na nieujemne liczby całkowite i że drugi przykład w rzeczywistości nie ma rozwiązania.
Sp3000,
intput 1 2 3 ma jednak prawidłowe wyjście ... [1, 1]
Jack Ammo
@JackAmmo: Wszystkie przykłady w drugim bloku kodu odpowiadają 2x + 3y = 1 .
Dennis,
W ax + bx = k wydaje mi się, że rozumiem, że rozwiązaniem musi być x> = 0 i y> = 0. Więc kim są takie x, y> = 0 rozwiązań 38 * x + 909 * y = 3?
RosLuP
W takim przypadku prawdopodobnie muszę zwrócić nieistniejące rozwiązanie ...
RosLuP

Odpowiedzi:

6

Pyth, 92 bajty

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

To całkiem potwór.

Wypróbuj online: demonstracja . Format wejściowy to, c\n[a,b]a format wyjściowy to [x,y].

W przypadku, gdy nie istnieje żadne rozwiązanie liczb całkowitych, nic nie wydrukuję, aw przypadku, gdy nie istnieje naturalne rozwiązanie liczb całkowitych, po prostu wydrukuję losowe rozwiązanie liczb całkowitych.

Wyjaśnienie (zgrubny przegląd)

  1. Najpierw znajdę całkowite rozwiązanie równania ax + by = gcd(a,b)za pomocą algorytmu Extended Euclidean.

  2. Następnie zmodyfikuję rozwiązanie (moje mnożenie ai bprzez c/gcd(a,b)), aby uzyskać całkowite rozwiązanie ax + by = c. Działa to, jeśli c/gcd(a,b)jest liczbą całkowitą. W przeciwnym razie nie ma rozwiązania.

  3. Wszystkie inne rozwiązania całkowite mają postać a(x+n*b/d) + b(y-n*a/d) = c z d = gcd(a,b)do liczby całkowitej n. Używając dwóch nierówności x+n*b/d >= 0i y-n*a/d >= 0mogę określić 6 możliwych wartości dla n. Spróbuję wszystkich 6 z nich i wydrukuję rozwiązanie o najwyższym najniższym współczynniku.

Objaśnienie (szczegółowe)

Pierwszym krokiem jest znalezienie rozwiązania liczby całkowitej w równaniu ax' + by' = gcd(a,b). Można tego dokonać za pomocą rozszerzonego algorytmu euklidesowego. Możesz dowiedzieć się, jak to działa na Wikipedii . Jedyna różnica polega na tym, że zamiast 3 kolumn ( r_i s_i t_i) użyję 6 kolumn ( r_i-1 r_i s_i-1 s_i t_i-1 t_i). W ten sposób nie muszę przechowywać ostatnich dwóch wierszy w pamięci, tylko ostatni.

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. 

Teraz chcę znaleźć rozwiązanie ax + by = c. Jest to możliwe tylko wtedy, gdy c mod gcd(a,b) == 0. Jeżeli równanie to jest spełnione, po prostu mnożąc x',y'się 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]

Mamy rozwiązanie dla liczb całkowitych ax + by = c. Zauważmy, że x, yalbo oba mogą być ujemne. Naszym celem jest więc przekształcenie ich w nieujemne.

Zaletą równań diofantycznych jest to, że możemy opisać wszystkie rozwiązania za pomocą tylko jednego początkowego rozwiązania. Jeśli (x,y)to rozwiązanie, że wszelkie inne rozwiązania są formy (x-n*b/gcd(a,b),y+n*a/gcd(a,b))do nliczby całkowitej.

Dlatego chcemy znaleźć n, gdzie x-n*b/gcd(a,b) >= 0i y+n*a/gcd(a,b >= 0. Po pewnej transformacji powstają dwie nierówności n >= -x*gcd(a,b)/bi n >= y*gcd(a,b)/a. Zauważ, że symbol nierówności może patrzeć w innym kierunku z powodu podziału z potencjalnym ujemnym alub b. Nie dbam o to tak bardzo, po prostu mówię, że jedna liczba -x*gcd(a,b)/b - 1, -x*gcd(a,b)/b, -x*gcd(a,b)/b + 1zdecydowanie spełnia nierówność 1, a jedna liczba y*gcd(a,b)/a - 1, y*gcd(a,b)/a, y*gcd(a,b)/a + 1spełnia nierówność 2. Jest taka n, która zaspokaja obie nierówności, jedna z 6 liczb również.

Następnie obliczam nowe rozwiązania (x-n*b/gcd(a,b),y+n*a/gcd(a,b))dla wszystkich 6 możliwych wartości n. I drukuję rozwiązanie o najwyższej najniższej wartości.

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

Sortowanie według sortowanej kolejności działa w następujący sposób. Korzystam z tego przykładu2x + 3y = 11

Sortuję każde z 6 rozwiązań (nazywane to kluczami), a oryginalne rozwiązania sortuję według kluczy:

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]

To sortuje kompletne nieujemne rozwiązanie do końca (jeśli istnieje).

Jakube
źródło
1
  • po uwagach Dennisa, które sprawiły, że mój poprzedni pomysł został wywrócony do góry nogami, musiałem zmienić kod z jego korzeni i zajęło mi to długoterminowe debugowanie i kosztowało mnie dwa razy więcej bajtów: ”(.

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
  • Wiem, że to nie gra w golfa, ponieważ tego rodzaju języki nie są przystosowane do zmniejszania długości kodu, ale mogę zapewnić, że złożoność czasu jest najlepsza.

Wyjaśnienie:

  • kod przyjmuje trzy wejściowe niezmienniki a, b, c, te ostatnie są poddane kilku warunkom przed przystąpieniem do obliczania:

    1- jeśli (a + b> c) i (a, b, c> 0) nie ma rozwiązania!

    2- jeśli (a + b <c), (a, b, c <0) brak rozwiązania!

    3- jeśli (a, b) mają wspólne przeciwne znaki c: brak rozwiązania!

    4 - jeśli GCD (a, b) nie podzieli c, to znowu nie ma rozwiązania! , w przeciwnym razie podziel wszystkie warianty przez GCD.

  • po tym musimy sprawdzić inny warunek, powinien on ułatwić i skrócić drogę do pożądanego rozwiązania.

    5- jeśli c podziel a lub b, rozwiązanie s = (x lub y) = (c- [ax, yb]) / [b, a] = C / [b, a] + [ax, yb] / [b , a] = S + [ax, yb] / [b, a], gdzie S jest naturalne, więc ax / b lub przez / a będą odtąd mieć nieujemne bezpośrednie rozwiązania, które odpowiednio wynoszą x = b lub y = a. (zauważ, że rozwiązania mogą być zerowymi wartościami w przypadku, gdy poprzednie arbitralne rozwiązania zostaną ujawnione jako negatywne)

  • gdy program osiąga ten etap, zamiana węższych rozwiązań dla x = (c-yb) / a jest zamiatana, dzięki zgodności, zamiatania większych zakresów liczb, które są powtarzane cyklicznie. największym polem wyszukiwania jest [xa, x + a], gdzie a jest dzielnikiem.

SPRÓBUJ

Abr001am
źródło
euuh, problem z dużymi liczbami, naprawię to (zastanawiam się, kiedy lol)
Abr001am
Myślę, że to wciąż niewielki błąd do naprawienia, dotyczący dużych liczb całkowitych, wciąż nie rozumiem, dlaczego dzielenie 1152921504606846800.000000 / 576460752303423420.000000 wychodzi z liczbą naturalną 2, chociaż ten ostatni wynik jest zaokrąglony.
Abr001am,
och zapomniałem naprawić ten błąd: p dzięki za zauważenie go @Jakube
Abr001am
0

Aksjomat, 460 bajtów

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

golf i trochę testu

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

W innych możliwych „rozwiązaniach” wystąpił błąd, ponieważ próbował zapisać nieskończone rozwiązania na jednej liście; teraz jest narzucony limit maksymalnie 80 rozwiązań

RosLuP
źródło