DIOFANTOVA
ROVNICE
 
 

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
NExistujeP(N)čas (sec)P(N) / P(N - 1)
2 1 0
3 1 01.00
4 + 3 0
5 3 0
6 5 01.67
7 8 01.60
8 + 7 0
9 270
10 4401.63
11 770.051.75
12 +14 0
13 265 0.11
14 429 0.161.60
15 772 0.281.80
 
NExistujeP(N)čas (sec)P(N) / P(N - 1)
16 + 25 0
17 2875 0.94
18 4742 1.601.65
19 87612.971.85
20 + 380.06
21 3367211.26
22 56392 19.661.67
23 105581 37.571.87
24 + 54 0
25 414524 143.74
26 703409 250.961.70
27 1329633 466.221.89
28 + 75 0.06
29 52969421794.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 ExistujePrů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 ExistujePrů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.071.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
1136213110 13913 16301 17058 18473 20594 20640 20886 25846
26730 2849530352 30646 32243 32762 34159 35037 36114 36359
39237 40442 4126442796 43803 44301 44400 44777 4572847467
49921 51136 51454 5167953961 54581 55479 57299 5753758066
59220 59308 59901 60066 6118161579 62559 62892 6390664905
65776 67353 67776 70076 70171 7131671796 75373 7755580520
81273 81681 81804 81815 83093 83490 8451384750 8493185057
88009 88180 88419 88569 89708 90075 90418 909929137592527
92877 93400 94319 94734 97172 98091 98170 98173 9833799363

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
NExistujePrůměr P(N) Maximum P(N)
10 52 98
20 22 534 40 733
30 20× 1 293 339 5 023 110
4020×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!.


hlavní stránka english
změněno 19. července 2001
Copyright © 1998-2001 Vladimír Zábrodský