Problém řešení Diofantovy lineární rovnice: Je dána množina přirozených čísel A[1], ..., A[N] a přirozené číslo C, nalezněte takovou podmnožinu A[1], ..., A[N], aby součet jejich prvků byl roven C. Tento problém je NP-úplný. Speciální případ pro sudé C, C = (A[1] + A[2] + ... + A[N]) div 2, je tzv. problém spravedlivého rozdělení (partition problem). Uvedu velice jednoduchý exponenciální algoritmus pro řešení těchto problémů. ... ani na nejrychlejším počítači na světě?     V literatuře se můžeme dočíst (následující citáty jsou vybrány z více různých zdrojů): Jediný známý způsob řešení je prozkoumávat postupně všechny podmnožiny, ale těch je u N-prvkové množiny 2N. Například pro N = 300 je počet všech podmnožin větší než odhadovaný celkový počet elektronů, protonů a netronů v celém vesmíru. Pokud by tuto úlohu řešil počítač, který by za sekundu prozkoumal miliardu podmnožin, nedokončil by výpočet dříve, než by vyhaslo naše slunce. Už spravedlivé rozdělení věcí z trochu větší kanceláře, N = 100, nemá smysl řešit ani na NEJRYCHLEJŠÍM POČÍTAČI NA SVĚTĚ. Algortimus DIOPHANT     Uvedu jednoduchý exponenciální algoritmus pro řešení Diofantovy lineární rovnice: Je dána množina přirozených čísel A[1], ..., A[N] a přirozené číslo C, nalezněte takovou podmnožinu A[1], ..., A[N], aby součet jejich prvků byl roven C. Algoritmus Utřiď prvky pole ve vzestupném pořadí. Vypočti Ls[1] = A[1], Ls[2] = Ls[1] + A[2], ..., Ls[J] = Ls[J - 1] + A[J], ... pro J = 3, ..., N. Je-li A[1] <= C & C <= Ls[N], pak Najdi v poli Ls první index K zleva takový, že Ls[K] >= C.   Je-li Ls[K] = C, pak je řešení nalezeno a do hledané podmnožiny patří i prvky A[1], ..., A[K]; algoritmus končí.   Najdi v poli A první index R zprava takový, že A[R] <= C.   Je-li A[R] = C, pak je řešení nalezeno a do hledané podmnožiny patří také prvek A[R]; algoritmus končí.   Pro L = K, ..., R, pokud A[1] <= C - A[L], budeme předpokládat, že prvek A[L] patří do hledané podmnožiny a v podmnožině A[1], ..., A[L - 1] budeme hledat ty prvky, jejichž součet je roven C - A[L]. Tj. rekurzivní volání této prosvětlené části popisu algoritmu s hodnotami: N = L - 1 a C = C - A[L]     V tomto článku představuji nerekurzívní implementaci algoritmu DIOPHANT v jazyce Rexx. DIOPHANT call SORT 1, N Ls.1 = A.1 do I = 2 to N; Im1 = I - 1; Ls.I = A.I + Ls.Im1; end if A.1 <= C & C <= Ls.N then do   S = 1; Stack.1 = N C   do while S > 0     parse var Stack.S R C V; S = S - 1     do K = 1 while Ls.K < C; end     if Ls.K = C then call EXISTUJE V, K, 1     do while A.R > C; R = R - 1; end     if A.R = C then call EXISTUJE V, R, R     do L = K to R while A.1 <= C - A.L       D = V A.L; S = S + 1; Stack.S = (L - 1) (C - A.L) D     end   end end say "Řešení neexistuje"   Řešení je ukládáno do proměnné V, jako řetězec čísel oddělených mezerami. Do zásobníku se ukládá počet prvků pole, hodnota C pro podproblém, řetězec dosud nalezených prvků řešení. Uložení hodnot do zásobníku - push operace - je implementována instrukcemi: S = S + 1; Stack.S = (L - 1) (C - A.L) D Pokud řešení existuje, vyvolá se procedura EXISTUJE. Uvedu její jednoduchý příklad, který jen vypíše řešení na obrazovce a ukončí výpočet: EXISTUJE EXISTUJE: procedure expose A. parse arg V, B, E do J = B to E by -1; V = V A.J; end say "Řešení je:" V exit       Pokud bychom instrukci exit zaměnili za instrukce S = S - 1; return, pak by algoritmus DIOPHANT hledal všechna řešení.   Nejhorší případ - Časová složitost O(2N)      Následující příklad dokazuje, že algoritmus DIOPHANT má v nejhorším případě exponenciální složitost. Sloupec označený P(N) obsahuje počty push operací. Exponenciální složitost řešení problému spravedlivého rozdělení A[1] = 2, A[J] = 100 + A[J - 1], J = 2, ..., N N Existuje P(N) čas (sec) P(N) / P(N - 1) 2   1 0   3   1 0 1.00 4 + 3 0   5   3 0   6   5 0 1.67 7   8 0 1.60 8 + 7 0   9   27 0   10   44 0 1.63 11   77 0.05 1.75 12 + 14 0   13   265 0.11   14   429 0.16 1.60 15   772 0.28 1.80 N Existuje P(N) čas (sec) P(N) / P(N - 1) 16 + 25 0   17   2875 0.94   18   4742 1.60 1.65 19   8761 2.97 1.85 20 + 38 0.06   21   33672 11.26   22   56392 19.66 1.67 23   105581 37.57 1.87 24 + 54 0   25   414524 143.74   26   703409 250.96 1.70 27   1329633 466.22 1.89 28 + 75 0.06   29   5296942 1794.41     K problému spravedlivého rozdělení         Můžeme ušetřit zhruba polovinu push operací, když budeme hledat jen jedinou podmnožinu, která obsahuje prvek A[N] a součet jejichž prvků má být roven ((A[1] + ... + A[N]) div 2) - A[N]. To je ovšem problém řešení Diofantovy linearní rovnice pro N - 1 prvků.     V tomto článku budu problém spravedlivého rozdělení řešit jako "obyčejný" problém řešení Diofantovy lineární rovnice a to i v případě, když (A[1] + ... + A[N]) div 2 bude liché číslo.   Průměrný? případ - Příklady polynomiální složitosti     Co je průměrný případ pro problém řešení Diofantovy lineární rovnice? Toť otázka! Test #1     V prvním pokusu jsem pro dané N = 10, 20, 30, 40, 50, 100, 200, 300, 400, 500 vygeneroval pole A pomocí následujícího fragmentu programu. Poznamenejme, že v jazyce Rexx vrací volání vestavěné funkce RANDOM(Min, Max) pseudonáhodná čísla z intervalu <Min, Max>, kde Max je nejvýše rovno 100000. Generování pole A Sum = 0 do J = 1 to N    A.J = RANDOM(1, 100000); Sum = Sum + A.J end U problému spravedlivého rozdělení jsem vypočetl C pomocí celočíselného dělení. V jazyce Rexx je operátorem celočíselného dělení %. C pro problém spravedlivého rozdělení C = Sum % 2 Pro problém řešení Diofantovy lineární rovnice jsem hodnotu C určil příkazem: C pro problém řešení Diofantovy rovnice C = FORMAT(Sum * RANDI(),,0) FORMAT je vestavěná funkce jazyka, jejíž volání v tomto případě vrací zaokrouhlenou celočíselnou hodnotu Sum * RANDI(). Volání funkce RANDI přitom vrací pseudonáhodné číslo v intervalu 0 až 1. Použil jsem funkci z článku Park S.K., Miler K.W.: Random Number Generators: Good ones are hard to find. CACM 1988, Vol. 31, No. 10, p. 1192-1201. V globální proměnné Seed je uložena běžná, naposledy určená, hodnota celočíselné posloupnosti. Generátor pseudonáhodných čísel musí být inicializován tak, že proměnné Seed se přiřadí hodnota z intervalu 1 až 2147483646. Pseudonáhodná čísla z intervalu 0 až 1 pak mohou být generována opakovaným voláním funkce RANDI. Poznámka: // je operátor operace zbytek po dělení. RANDI RANDI: procedure expose Seed numeric digits 14; Seed = RANDOM(1, 100000) A = 16807; M = 2147483647; Seed = (A * Seed) // M return Seed / M     Pro dané N jsem 100krát vygeneroval pole A, vypočetl C a řešil problém spravedlivého rozdělení; 100krát jsem vygeneroval pole A a C a řešil problém řešení Diofantovy lineární rovnice. Výsledky shrnuje následující tabulka. Doba výpočtu se týká PC s procesorem Cyrix 6x86MX a 32MB RAM. Spravedlivé rozdělení C = Sum % 2 N  Existuje Prům. P(N) Prům. čas (sec) 10 0 53.3 0.0205   20 86 7314.93 2.6880   30 100 3493.63 1.3191   40 100 2478.58 0.9303   50 100 2421.02 0.9429   100 100 2035.57 0.5999   200 100 5149.09 1.1845   300 100 10868.98 2.6189   400 100 18804.96 5.5589   500 100 29206.94 12.7907     Řešení Diofantovy rovnice C = FORMAT(Sum * RANDI(), , 0) N  Existuje Prům. P(N) Prům. čas (sec) 10 0 29.42 0.0109   20 53 4673.77 1.5646   30 86 3480.79 1.1818   40 97 3022.63 1.0483   50 94 1956.00 0.6834   100 99 4050.09 0.8992   200 100 5149.09 1.1845   300 100 7814.07 1.9347   400 100 14301.56 4.3563   500 100 21542.35 62.5650       Prohlédněme si (typický) příklad řešení Řešení problému spravedlivého rozdělení pro N = 100 2501 push operací, 0.83 sec 430 1989 2414 3616 4786 4851 5081 6774 7651 9716 11362 13110 13913 16301 17058 18473 20594 20640 20886 25846 26730 28495 30352 30646 32243 32762 34159 35037 36114 36359 39237 40442 41264 42796 43803 44301 44400 44777 45728 47467 49921 51136 51454 51679 53961 54581 55479 57299 57537 58066 59220 59308 59901 60066 61181 61579 62559 62892 63906 64905 65776 67353 67776 70076 70171 71316 71796 75373 77555 80520 81273 81681 81804 81815 83093 83490 84513 84750 84931 85057 88009 88180 88419 88569 89708 90075 90418 90992 91375 92527 92877 93400 94319 94734 97172 98091 98170 98173 98337 99363 5552260  =  2776130  +  2776130   Test #2 V druhém testu jsem generoval prvky pole A pomocí následujícího fragmentu programu (jeho první instrukce byla numeric digits 16): Generování prvků pole A v testu #2 Sum = 0 do J = 1 to N   A.J = FORMAT(RANDI() * 21474836.47,,0)   Sum = Sum + A.J end Řešil jsem problém spravedlivého rozdělení věcí s maximální cenou více než 21 milionů Kč. Proměnné Seed jsem přiřadil hodnotu 214748364 a pak jsem spustil 20krát DIOPHANT pro každou z hodnot N = 100, 200, 300, 400, 500. Později jsem Seed přiřadil hodnotu 19685089 a pokračoval jsem pro N = 30, 40, 50, 20, 10. Výsledky shrnuje následující tabulka: Problém spravedlivého rozdělení #2 N Existuje Průměr P(N)   Maximum P(N)   10 0× 52   98   20 1× 22 534   40 733   30 20× 1 293 339   5 023 110   40 20× 777 078   4 592 315   50 20× 681 585   2 797 498   100 20× 117 536   435 236   200 20× 88 077   454 263   300 20× 81 816   263 396   400 20× 83 037   344 113   500 20× 64 972   151 288   Historie algortimu DIOPHANT     Rekurzívní verzi DIOPHANT v jazyce Pascal, jsem poprvé publikoval v článku Problém dvou loupežníků s podtitulem Praktická poznámka k jednomu prakticky neřešitelnému problému v časopise Bajt, říjen 1993 (36), s. 134-136. V roce 1993 jsem se zabýval pouze řešením problému spravedlivého rozdělení. Nejzajímavějším z popsaných testů praktické použitelnosti programu DIOPHANT v Bajtu bylo spravedlivé rozdělení zboží ze skladu ČKD Blansko a.s.: 10000 položek v cenách, po vynásobení stem, od 2 Kč do téměř 400 miliónů Kč. Tisíckrát jsem z nich náhodně vybral 500 a program DIOPHANT se je pokusil spravedlivě rozdělit. Jen jednou jsem musel výpočet po devíti hodinách násilně přerušit. V ostatních případech bylo řešení vždy nalezeno; nejdéle za 96 minut, v průměru za 16.5 sekundy, v 967 případech ne za déle než 5 sekund (PC AT 386/40 MHz, 8MB RAM, Turbo Pascal 6.0).   Závěr Pokud pole A bude obsahovat hodnoty skutečných věcí (okolo 500 hodnot věcí z kanceláře, garáže, hangáru, doku) můžeme najít řešení Diofantovy lineární rovnicetéměř vždy, v rozumném čase a jen s pomocí interpretru jazyka Rexx a obyčejného PC!.