Persamaan Diophantine linear alami

13

Persamaan Diophantine linier dalam dua variabel adalah persamaan bentuk ax + by = c , di mana a , b dan c adalah bilangan bulat konstan dan x dan y adalah variabel bilangan bulat.

Untuk banyak persamaan Diophantine yang terjadi secara alami, x dan y mewakili jumlah yang tidak boleh negatif.

Tugas

Tulis sebuah program atau fungsi yang menerima koefisien a , b dan c sebagai input dan mengembalikan sepasang bilangan asli (0, 1, 2, ...) x dan y yang memverifikasi persamaan kapak + dengan = c , jika pasangan tersebut ada

Aturan tambahan

  • Anda dapat memilih format apa pun untuk input dan output yang hanya melibatkan bilangan bulat yang diinginkan dan, secara opsional, larik / daftar / matriks / notasi / vektor notasi bahasa Anda, selama Anda tidak menanamkan kode apa pun dalam input.

  • Anda dapat mengasumsikan bahwa koefisien a dan b keduanya bukan nol.

  • Kode Anda harus bekerja untuk setiap triplet bilangan bulat antara -2 60 dan 2 60 ; itu harus selesai di bawah satu menit di komputer saya (Intel i7-3770, 16 GiB RAM).

  • Anda tidak boleh menggunakan built-in yang memecahkan persamaan Diophantine dan karenanya meremehkan tugas ini, seperti Mathematica FindInstanceatau FrobeniusSolve.

  • Kode Anda mungkin berperilaku namun Anda inginkan jika tidak ada solusi yang ditemukan, asalkan sesuai dengan batas waktu dan hasilnya tidak dapat dikacaukan dengan solusi yang valid.

  • Aturan standar berlaku.

Contohnya

  1. Contoh di bawah ini menggambarkan I / O yang valid untuk persamaan 2x + 3y = 11 , yang memiliki tepat dua solusi yang valid ( (x, y) = (4,1) dan (x, y) = (1,3) ).

    Input:  2 3 11
    Output: [4 1]
    
    Input:  (11 (2,3))
    Output: [3],(1)
    
  2. Satu-satunya solusi yang valid dari 2x + 3y = 2 adalah pasangan (x, y) = (1,0) .

  3. Contoh di bawah ini menggambarkan I / O yang valid untuk persamaan 2x + 3y = 1 , yang tidak memiliki solusi yang valid .

    Input:  (2 3 1)
    Output: []
    
    Input:  1 2 3
    Output: -1
    
    Input:  [[2], [3], [1]]
    Output: (2, -1)
    
  4. Untuk (a, b, c) = (1152921504606846883, -576460752303423433, 1) , semua solusi yang benar (x, y) memenuhi bahwa (x, y) = (135637824071393749 - bn, 271275648142787502 + an) untuk beberapa non-negatif bilangan bulat n .

Dennis
sumber
Saya pikir mungkin lebih baik untuk menempatkan sedikit lebih banyak penekanan pada bilangan bulat negatif, dan bahwa contoh kedua sebenarnya tidak memiliki solusi.
Sp3000
intput 1 2 3 memiliki output yang valid ... [1, 1]
Jack Ammo
@JackAmmo: Semua contoh di blok kode kedua sesuai dengan 2x + 3y = 1 .
Dennis
Dalam ax + bx = k menurut saya untuk memahami bahwa solusinya harus x> = 0 dan y> = 0. Jadi, siapakah x, y> = 0 solusi dari 38 * x + 909 * y = 3?
RosLuP
Dalam kasus seperti itu mungkin saya harus mengembalikan solusi yang tidak ada ...
RosLuP

Jawaban:

6

Pyth, 92 byte

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

Itu monster.

Cobalah online: Peragaan . Format input adalahc\n[a,b] dan format output [x,y].

Dalam hal tidak ada solusi integer, saya tidak akan mencetak apa pun, dan dalam hal tidak ada solusi integer alami, saya hanya akan mencetak solusi integer acak.

Penjelasan (Ikhtisar Kasar)

  1. Pada awalnya saya akan menemukan solusi integer untuk persamaan ax + by = gcd(a,b)dengan menggunakan algoritma Extended Euclidean.

  2. Kemudian saya akan memodifikasi solusi (pengali adan bdengan saya c/gcd(a,b)) untuk mendapatkan solusi integer ax + by = c. Ini berfungsi, jika c/gcd(a,b)bilangan bulat. Kalau tidak, tidak ada solusi.

  3. Semua solusi integer lainnya memiliki formulir a(x+n*b/d) + b(y-n*a/d) = c dengan d = gcd(a,b)untuk integer n. Menggunakan dua ketidaksetaraan x+n*b/d >= 0dan y-n*a/d >= 0saya bisa menentukan 6 nilai yang mungkin untuk n. Saya akan mencoba semuanya dan mencetak solusinya dengan koefisien terendah tertinggi.

Penjelasan (Detail)

Langkah pertama adalah menemukan solusi integer untuk persamaan tersebut ax' + by' = gcd(a,b). Ini dapat dilakukan dengan menggunakan algoritma Euclidean yang diperluas. Anda bisa mendapatkan ide tentang cara kerjanya di Wikipedia . Satu-satunya perbedaan adalah, bahwa alih-alih menggunakan 3 kolom ( r_i s_i t_i) saya akan menggunakan 6 kolom ( r_i-1 r_i s_i-1 s_i t_i-1 t_i). Dengan cara ini saya tidak harus menyimpan dua baris terakhir dalam memori, hanya yang terakhir.

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. 

Sekarang saya ingin mencari solusinya ax + by = c. Ini hanya mungkin, ketika c mod gcd(a,b) == 0. Jika persamaan ini puas, saya hanya mengalikan x',y'dengan 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]

Kami memiliki solusi integer untuk ax + by = c. Perhatikan, bahwa x, yatau keduanya mungkin negatif. Jadi tujuan kami adalah mengubahnya menjadi non-negatif.

Hal yang menyenangkan tentang persamaan Diophantine adalah, kita dapat menggambarkan semua solusi hanya dengan menggunakan satu solusi awal. Jika (x,y)solusi, bahwa semua solusi lainnya adalah dari bentuk (x-n*b/gcd(a,b),y+n*a/gcd(a,b))untuk ninteger.

Karena itu kami ingin mencari n, di mana x-n*b/gcd(a,b) >= 0dan y+n*a/gcd(a,b >= 0. Setelah beberapa transformasi kita berakhir dengan dua ketidaksetaraan n >= -x*gcd(a,b)/bdan n >= y*gcd(a,b)/a. Perhatikan bahwa simbol ketimpangan mungkin melihat ke arah lain karena pembagian dengan potensi negatif aatau b. Saya tidak terlalu peduli tentang hal itu, saya hanya mengatakan bahwa satu angka -x*gcd(a,b)/b - 1, -x*gcd(a,b)/b, -x*gcd(a,b)/b + 1pasti memuaskan ketidaksetaraan 1, dan satu angka y*gcd(a,b)/a - 1, y*gcd(a,b)/a, y*gcd(a,b)/a + 1memuaskan ketidaksetaraan 2. Ada n, yang memenuhi kedua ketidaksetaraan, salah satu dari 6 angka juga.

Lalu saya menghitung solusi baru (x-n*b/gcd(a,b),y+n*a/gcd(a,b))untuk semua 6 nilai yang mungkin dari n. Dan saya mencetak solusinya dengan nilai terendah tertinggi.

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

Urutkan berdasarkan urutan pesanan mereka berfungsi dengan cara berikut. Saya menggunakan contohnya2x + 3y = 11

Saya mengurutkan masing-masing dari 6 solusi (ini disebut kunci), dan mengurutkan solusi asli dengan kunci mereka:

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]

Ini semacam solusi non-negatif lengkap sampai akhir (jika ada).

Jakube
sumber
1
  • setelah komentar Dennis, yang membuat ide saya sebelumnya terbalik, saya harus mengubah kode dari akarnya dan butuh saya debugging jangka panjang, dan biaya saya dua kali byte byte: '(.

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
  • Yah, saya tahu itu tidak golf, karena jenis bahasa tidak disesuaikan untuk pengurangan panjang kode, tetapi, saya dapat memastikan bahwa kompleksitas waktu adalah yang terbaik.

Penjelasan:

  • kode mengambil tiga invarian a, b, c sebagai input, yang terakhir ini ditundukkan ke beberapa kondisi sebelum melanjutkan untuk menghitung:

    1- jika (a + b> c) dan (a, b, c> 0) tidak ada solusi!

    2- jika (a + b <c), (a, b, c <0) tidak ada solusi!

    3 - jika (a, b) memiliki tanda kebalikan yang sama dari c: tidak ada solusi!

    4 - jika GCD (a, b) dosnt membagi c, maka tidak ada solusi lagi! , jika tidak, bagi semua varian dengan GCD.

  • setelah ini, kita harus memeriksa kondisi lain, itu harus memudahkan dan mempersingkat jalan menuju solusi yang diinginkan.

    5- jika c membagi a atau b, solusi s = (x atau y) = (c- [kapak, yb]) / [b, a] = C / [b, a] + [kapak, yb] / [b , a] = S + [kapak, yb] / [b, a] di mana S adalah alami sehingga kapak / b atau oleh / a akan memiliki solusi langsung selanjutnya non-negatif yang masing-masing x = b atau y = a. (perhatikan bahwa solusi dapat hanya bernilai nihil jika solusi arbitrer sebelumnya dinyatakan negatif)

  • ketika program mencapai tahap ini, serangkaian solusi yang lebih sempit untuk x = (c-yb) / a malah disapu, berkat kongruensi, menyapu rentang angka yang lebih besar, yang datang secara berulang-ulang dengan siklus reguler. bidang pencarian terbesar adalah [xa, x + a] dengan a adalah pembagi.

COBALAH

Abr001am
sumber
euuh, masalah jumlah besar, akan memperbaikinya (heran ketika lol)
Abr001am
Saya pikir bug masih kecil untuk memperbaiki, tentang bilangan bulat besar, saya masih tidak mengerti mengapa membagi 1152921504606846800.000000 / 576460752303423420.000000 keluar dengan bilangan alami 2, meskipun hasil terakhir ini dibulatkan.
Abr001am
ohh saya lupa untuk memperbaiki bug yang mengganggu itu: p terima kasih telah memperhatikannya @ Jakube
Abr001am
0

Aksioma, 460 byte

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 dan beberapa tes

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

Di 'solusi' lain yang mungkin ada bug karena mencoba menyimpan solusi tak terbatas dalam satu Daftar; sekarang diberlakukan batas 80 solusi maks

RosLuP
sumber