Niin, siis aloin tuossa eräs päivä opiskelemaan säikeitä ja päätin tehdä pienen ns. benchmark-ohjelman joka laskee alkulukuja. Brute-force -menetelmä kävi hitaaksi joten päätin tutkia josko sattuisi olemaan nopeampaa algoritmia tarjolla ja sellainenhan löytyi, nimittäin Eratostheneen seula.
Wikipediasta:
Seulan toimintaperiaate (algoritmi) on seuraava:
1. Kirjoitetaan lista kaikista luonnollisista luvuista alkaen kakkosesta ja päättyen johonkin valittuun suurimpaan lukuun n.
2. Poistetaan listasta kaikki luvun 2 monikerrat (4, 6, 8 jne.).
3. Listan seuraava jäljellä oleva luku on alkuluku.
4. Poistetaan listasta kaikki ne luvut, jotka ovat sekä edellisessä vaiheessa löydettyä alkulukua suurempia että sen monikertoja.
5. Toistetaan vaiheita 3 ja 4, kunnes listan seuraava jäljellä oleva luku on suurempi kuin listan suurimman luvun n neliöjuuri.
6. Nyt listassa on jäljellä vain alkulukuja.
Ongelma onkin siinä, että tämä listan koko on vakio. Se ei sinänsä tuota ongelmaa, mutta ideaalinen algoritmi tässä tapauksessa toimisi siten, että laskettavien alkulukujen määrä olisi loputon(eli selvillä olevien alkulukujen määrä kasvaisi sitä mukaa kun ohjelma rullaa), tai, että itse työ voitaisiin jakaa haluttuun määrään osia. Esim. siten, että säie #0 seulaa luvut 0-1000, säie #1 seulaa luvut 1000-2000, säie #2 seulaa luvut 3000-4000, säie #n seulaa luvut (1000*n)-(1000*n+1000) jne. Onko tämä jollain tapaa mahdollista tällaista seula-algoritmia käyttäen?
Tässä vielä implementaatio algoritmista:
/* 1. Create a list of consecutive integers from two to n: (2, 3, 4, ..., n). * 2. Initially, let p equal 2, the first prime number. * 3. Strike from the list all multiples of p greater than p. * 4. Find the first number remaining on the list greater than p (this number is the next prime); let p equal this number. * 5. Repeat steps 3 and 4 until p2 is greater than n. * 6. All the remaining numbers on the list are prime. */ #include <iostream> #define MAX 3580 int main() { std::cout << "Finding all the primes less than " << MAX << "..." << std::endl; // 1. unsigned int arr[MAX]; unsigned int n=2; while(n < MAX) { arr[n] = n; n++; } // 2. unsigned long long int p; p = 2; // 3. while(p*p < MAX) { n=p+1; while(n<MAX) { if(n%p==0) arr[n] = 0; n++; } // 4. n=p+1; while(1) { if(arr[n] > 0) { p = arr[n]; break; } n++; } } n=0; while(n<MAX) { if(arr[n] > 0) std::cout << arr[n] << std::endl; n++; } std::cout << "Done!" << std::endl; }
Saako tämän algoritmin räätälöityä jotenkin siten, että algoritmi olisi ns. "online", eli etsittävien alkulukujen määrä kasvaisi ohjelman rullatessa, vai jääkö brute-force -menetelmän optimointi ainoaksi vaihtoehdoksi?
Itse säikeistäminen ei ole niinkään ongelma, kun voi suorittaa tuon lukujen kerronnaisten etsimisen ja alkioiden poistamisen eri säikeissä.
Claw kirjoitti:
Ongelma onkin siinä, että tämä listan koko on vakio. Se ei sinänsä tuota ongelmaa, mutta ideaalinen algoritmi tässä tapauksessa toimisi siten, että laskettavien alkulukujen määrä olisi loputon(eli selvillä olevien alkulukujen määrä kasvaisi sitä mukaa kun ohjelma rullaa)
Loputon alkulukujen määrä tuolla algoritmilla käyttäisi rajattomasti muistia ja rajattomasti prosessointiaikaa. Toki voit laskea ensin vaikka lukuihin 1000 asti, sitten kun ne on valmiit kasvattaa määrää 10000:een ja suorittaa vaiheet luvusta 2 alkaen uudelleen luvuille 1001-10000.
Claw kirjoitti:
tai, että itse työ voitaisiin jakaa haluttuun määrään osia.
Tokihan työ voidaan hyvin helposti jakaa haluttuun määrään osia, kuvaamallasi tavalla. Eli käsket 1. säikeen poistaa kakkoset ensimmäisestä 1/n lukujoukosta, 2. säikeen poistaa kakkoset 2. 1/n lukujoukosta jne.
Ohjelman vaihe 3 (ei-alkulukujen poisto listasta) on nyt toteutettu niin epätehokkaasti, että algoritmi ei tuossa muodossa ole paljoa brute forcea tehokkaampi.
Listalta halutaan poistaa luvun p moninkerrat, eli luvut 2p, 3p, 4p... Turha siis kokeilla joka luvulla onko luku p:n moninkerta, kun kyseiset luvut voidaan löytää suoraankin:
n = p+p; while(n<MAX) { arr[n] = 0; n += p; }
Koodia voidaan edelleen nopeuttaa hieman alustamalla n aluksi arvoon p*p, koska kaikki sitä pienemmät ei-alkuluvut on jo poistettu, sillä niillä on varmasti jokin p:tä pienempi alkutekijä.
Sitten varsinaiseen kysymykseesi: seula on mahdollista saada tuottamaan äärettömästi (tai niin pitkälle kuin lukutyypin koko riittää) alkulukuja viivyttämällä ei-alkulukujen karsintaa niin, että kukin ei-alkuluku poistetaan vasta, kun tarkastellaan, onko kyseinen luku alkuluku.
Kun luku x todetaan alkuluvuksi, Eratostheneen seula sulkee pois luvut 2x, 3x, 4x... Halutaan kyetä tuottamaan äärettömästi alkulukuja, joten ei voida käydä näitä kaikkia lukuja läpi. Sen sijaan lukujono voidaan tallentaa myöhempää tarkastelua varten lukuparina (n,p), jossa n tarkoittaa lukujonon ensimmäistä arvoa ja p arvoa, jolla n:ää kasvatetaan seuraavan alkion tuottamiseksi, eli tässä kohtaa n=2x ja p=x. Kun lukuja käydään läpi järjestyksessä, voidaan lukupari (n,p) korvata lukuparilla (n+p,p), kun kaikki lukua n pienemmät arvot on saatu käsiteltyä.
Kun ollaan käytä läpi kaikki lukua y pienemmät luvut ja halutaan tietää, onko y alkuluku, valitaan kaikista tallennetuista lukujonoista (n,p) se, jonka alkuarvo n on pienin. Tarkastellaan arvoja n ja y:
- Jos n>y, y on varmasti alkuluku, koska mikään aiemmin löydetty alkuluku ei ole sulkenut lukua y pois. Tällöin lukupari (2y,y) lisätään lukujonojen joukkoon jotta myöhemmin vastaan tulevat y:n tekijät osataan sulkea pois.
- Jos n=y, y ei ole alkuluku, sillä p on sen tekijä.
- Jos n<y, on kaikki n:ää pienemmät luvut käsitelty, joten voidaan siirtyä p:n tuottamassa lukujonossa eteenpäin korvaamalla pari (n,p) parilla (n+p,p). Tällöin ei vielä osata sanoa, onko y alkuluku, vaan on etsittävä uudelleen lukujono, jonka aloituskohta on pienin, ja testattava sillä.
Meni varmaan aika sekavaksi tämä selitys, mutta toivottavasti saat koodista jotain selvää (toivon mukaan std::priority_queue on tuttu rakenne):
void generoi() { typedef std::pair<long long,long long> Lukujono; // Säilytetään jonoja prioriteettijonossa, josta pienimmän alkion löytäminen onnistuu tehokkaasti. std::priority_queue<Lukujono,std::vector<Lukujono>,std::greater<Lukujono> > jonot; // hoidetaan luku 2 erikoistapauksena std::cout << 2 << std::endl; jonot.push(Lukujono(4,2)); long long i=3; while(true) { Lukujono pienin = jonot.top(); if (pienin.first > i) { // i on alkuluku std::cout << i << std::endl; jonot.push(Lukujono(2*i,i)); ++i; } else if (pienin.first==i) { // i ei ole alkuluku ++i; } else { // ei tiedetä vielä, onko i alkuluku jonot.pop(); jonot.push(Lukujono(pienin.first+pienin.second, pienin.second)); } } }
Communicating Sequential Processes (CSP) - tyylinen ratkaisu olisi varmaan yksinkertaisin. Samalla se toimii hyvänä esimerkkinä, miten säikeiden tai prosessien avulla voidaan ongelmaa yksinkertaistaa ja saadaan niistä hyötyä, vaikka ne eivät edes toimisi yhtäaikaa.
En tiedä, miten alla oleva toteutetaan helpoiten C:llä säikeiden kanssa ilman kanavia, joita olen tottunut käyttämään Infernon Limbolla.
Elikkä, ensimmäinen prosessi eliminoi 2:lla jaolliset, seuraava 3:lla jaolliset, sitä seurava 5:llä jaolliset ja niin edespäin.
--------- --------- --------- ---------- - 2 ->| print 2 | | | | | | | - 3 ->| |- 3 ->| print 3 | | | | | - 4 ->| drop | | | | | | | - 5 ->| |- 5 ->| |- 5 ->| print 5 | | | - 6 ->| drop | | | | | | | - 7 ->| |- 7 ->| |- 7 ->| |- 7 ->| | - 8 ->| drop | | | | | | | - 9 ->| |- 9 ->| drop | | | | | - 10 ->| drop | | | | | | | - 11 ->| |- 11 ->| |- 11 ->| |- 11 ->| print 11 | ... | ... | ... | ... | ... | ... | ... | ... |
jalski: Ehdotuksesi on varmaan mahtavan tehokas, kun muistetaan, että esimerkiksi jo alle miljoonan suuruisia alkulukuja on 78498. Ajattelitko tehdä jokaiselle näistä oman prosessin? Paljonko muistia kuluu? Paljonko laskenta-aikaa menee hukkaan tarpeettoman raskaiden rakenteiden vuoksi?
Sisuaskin ratkaisu toteuttaa käsittääkseni suunnilleen saman algoritmin – järkevillä resursseilla.
Kiitokset kaikille vastauksista.
Sain nyt rutistettua tuota algoritmia vähän siedettävämpään muotoon;
#include <vector> #define MAX 1024*1024*1024 int main() { std::vector<bool> arr((size_t)MAX,1); unsigned int p = 2; unsigned int n = 0; for(p=2;p*p<MAX;p++) { if(arr[p]) { for(n=p*p;n<MAX;n+=p) { arr[n] = 0; } } } }
Tuo lienee jo aika pitkälle loppuunsa veivattu jo? Unix-työkalu time:n mukaan 1.6 GHz Celeron M seuloo nuo alkuluvut reilun miljardin sisältä n. 59 sekunnissa, kun käyttää g++:n -O3 lippua.
Seuraavana lienee tuo Sisuaskin ehdotuksen räätälöinti. Tässä nyt vielä 28 tuntia aikaa ennen asepalveluksen alkamista. Hikeä pukkaa, meinaa kiire tulla.
Metabolix kirjoitti:
jalski: Ehdotuksesi on varmaan mahtavan tehokas, kun muistetaan, että esimerkiksi jo alle miljoonan suuruisia alkulukuja on 78498. Ajattelitko tehdä jokaiselle näistä oman prosessin? Paljonko muistia kuluu? Paljonko laskenta-aikaa menee hukkaan tarpeettoman raskaiden rakenteiden vuoksi?
Sisuaskin ratkaisu toteuttaa käsittääkseni suunnilleen saman algoritmin – järkevillä resursseilla.
En väittänyt, että ehdottamani toteutus olisi nopein tai tehokkain laskettaessa ääretöntä määrää alkulukuja. Infernon prosessi on kevyt ja laskettaessa esim. muutamaa tuhatta ensimmäistä alkulukua, niin suorituskyky on varmasti riittävä vanhalla ja hitaallakin koneella.
Tarkoitus oli vain näyttää eri lähestymistapa ongelmaan. Et voi kieltää, etteikö CSP-toteutus olisi yksinkertainen ja helppo hahmottaa. Toteutus on lyhyt ja sen kirjoittaa nopeasti, testaus ja debuggaus onnistuu muutamassa minuutissa.
Ratkaisu, jossa säikeiden määrä kasvaa (rajatta) ongelman koon mukaan, ei ole järkevää rinnakkaisohjelmointia. Myöskään ohjelman rakenteen selkeyttäminen säikeiden avulla ei tarkoita mielivaltaisten osaongelmien ripottelua omiin säikeisiin riippumatta siitä, sisältävätkö ne järkevässä määrin riippumatonta laskentaa (kuten tässä tapauksessa eivät).
Käyttöjärjestelmälle aiheutuva kuorma on kohtuuton, ohjelman suorituskyky erittäin heikko ja sen arviointi vaikeaa. Tarpeettomasta säikeiden rohmuamisesta voi aiheutua haittaa myös muille prosesseille.
Laskettaessa esim. muutamaa tuhatta ensimmäistä alkulukua pääsee kaikkein helpommalla kaikkein triviaaleimmallakin algoritmilla eli kokeilemalla jakoa kaikilla aiemmilla luvuilla. Siis esimerkiksi Pythonilla hyvinkin lyhyesti:
alkuluvut = [2, 3] for ehdokas in range(5, 10000, 2): ok = True for jakaja in alkuluvut: if ehdokas % jakaja == 0: ok = False if not ok or jakaja * jakaja > ehdokas: break if ok: alkuluvut.append(ehdokas)
Kirjoittaminen kesti alle minuutin, testaus kesti kymmenisen sekuntia kaikkineen, algoritmi toimi hetkessä. Ideakin on luultavasti paljon helpompi hahmottaa tästä kuin CSP-toteutuksesta.
Grez kirjoitti:
Toki voit laskea ensin vaikka lukuihin 1000 asti, sitten kun ne on valmiit kasvattaa määrää 10000:een ja suorittaa vaiheet luvusta 2 alkaen uudelleen luvuille 1001-10000.
Tämä vaihtoehto vaikutti järkevimmältä tähän hätään.
Nykyisellä viritelmällä työn hoitaminen pienissä osissa kyllä onnistuu:
#define MAX 100 std::vector<bool> arr(MAX, 1); int p=0; int n=0; int list = 0; for(int loops = 0; list < MAX; loops++) { list = 25+loops*25; for(p=2;p*p<list;p++) { if(arr[p]) { for(n = p + p;n < list; n+=p) { arr[n] = 0; } } } }
mutta en saa päähäni mitään järkevää tapaa siihen, että samoja lukuja ei "tiputettaisi"(eli asetettaisi nollaksi) aina uudelleen. Nimittäin jos yrittää tuohon lisätä jonkinlaista offset(loops*25)-kikkailua siten, että jokaisella loopilla muokataan vaan niitä lisättyjä lukuja, ilmenee ainakin yksi ongelma:
arr[n+offset] = 0; -> n on jo korotettu for-lauseessa yhden p:n verran, eli jokaisen "lisätyn" lukumäärän ensimäisen pudotettavan kerroin jää ykköseksi.
Koitin sitten vääntää purkkaa, joka sekään ei tunnu oikein tuottavan tulosta;
#define MAX 100 std::vector<bool> arr(MAX, 1); int p=0; int n=0; int list = 0; int sequence = 20; for(int loops = 0; list < MAX; loops++) { list = sequence+loops*sequence; for(p=2;p*p<list;p++) { if(arr[p]) { if(loops==0) { // Tarkistetaan for(n=p+p;n<list;n+=p) { arr[n] = 0; } } else { /* Jos hypätään "uusille" luvuille, merkataan myös ensimäinen luku, ja lisätään mukaan offsetti */ for(n = p + loops*sequence;n < list; n+=p) { arr[n] = 0; } } } } }
koodi sattumoisin jättää joitain lukuja pudottamatta, tai pudottelee vääriä lukiuja. Esim. osa vöitetyistä alkuluvuista on jaollisia kahdella.
Tekee jo mieli luovuttaa kuuden tunnin pähkäilyn jälkeen. Luultavimmin vika on kuitenkin jossain liiankin yksinkertaisessa paikassa.
os kirjoitti:
Myöskään ohjelman rakenteen selkeyttäminen säikeiden avulla ei tarkoita mielivaltaisten osaongelmien ripottelua omiin säikeisiin riippumatta siitä, sisältävätkö ne järkevässä määrin riippumatonta laskentaa (kuten tässä tapauksessa eivät).
Mielivaltaisten osaongelmien ripottelua... ? Tajusitkohan nyt ollenkaan laittamaani kuvaa ja sen ajatusta? Säikeisiinhän ei tässä ole ripoteltuna mitään mielivaltaisia osaongelmia. Kyseessähän on täysin lineaarinen suoritus putki, jossa säie toimii käytännössä suodattimena. Kaikki tapahtuu järjestyksessä, suorituksessa ei ole mitään satunnaista tai mielivaltaista.
Metabolix: Tarkoituksena oli yksinkertaistaa Eratostheneen seulan algoritmia. Tarkoituksenani ei ollut kiistellä siitä, mikä on nopein tai paras tapa.
Tiedän kyllä, että CSP-toteutus käy raskaaksi, jos ratkaistavana on miljoonia alkulukuja.
Jos jotain oikeasti kiinnostaa tutustua CSP-ohjelmointimalliin, niin voi aloittaa tutustumisen vaikka seuraavasta linkistä: http://swtch.com/~rsc/thread/squint.pdf
Kokeilin CSP-toteutusta alkulukujen etsimisestä miljoonaan asti alla olevalla koodilla:
implement Primes; include "draw.m"; Primes: module { init: fn(nil: ref Draw->Context, argl: list of string); }; include "sys.m"; sys: Sys; MAX: con 1000000; BUFSZ: con 256; init(nil: ref Draw->Context, nil: list of string) { sys = load Sys Sys->PATH; c := chan[BUFSZ] of int; spawn prime(c); for(n := 2; n <= MAX; n++) c <-= n; c <-= 1; } prime(c: chan of int) { p := <-c; if(p == 1) exit; sys->print("%d\n", p); nc := chan[BUFSZ] of int; spawn prime(nc); for(;;){ n := <-c; if(n%p) nc <-= n; if(n == 1) exit; } }
Claw kirjoitti:
Tuo lienee jo aika pitkälle loppuunsa veivattu jo?
Yksi helppo optimointi on käydä läpi pelkästään parittomat luvut, jolloin taulukon koko puolittuu. Pienellä lisävaivalla voi jättää pois myös kolmella jaolliset luvut, joka tiputtaa muistinkäytön kolmasosaan alkuperäisestä. Tässä tapauksessa käydään siis läpi vain luvut 5, 7, 11, 13, 17, 19, ... eli ne jotka ovat 1 tai 5 mod 6. Luvut 2 ja 3 voi tulostaa alussa erikoistapauksena. Tätä ideaa voi periaatteessa jatkaa pidemmällekin, mutta homma menee nopeasti aika hankalaksi ja hyöty pienenee.
Kun nyt triviaaleista ja hitaista algoritmeistä puhutaan, 1000 ensimmäistä alkulukua löytyy alle sekunnissa tällä:
primes = go [2..] where go (p:xs) = p : go (filter ((> 0) . (`mod` p)) xs) main = mapM_ print $ take 1000 primes
Ylärajattoman Eratostheneen seulan voi toteuttaa Haskellilla aika elegantisti: http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf (PDF). Se on tosin huomattavasti hitaampi kuin optimoitu bittitaulukkoa käyttävä toteutus.
jalski kirjoitti:
os kirjoitti:
Myöskään ohjelman rakenteen selkeyttäminen säikeiden avulla ei tarkoita mielivaltaisten osaongelmien ripottelua omiin säikeisiin riippumatta siitä, sisältävätkö ne järkevässä määrin riippumatonta laskentaa (kuten tässä tapauksessa eivät).
Mielivaltaisten osaongelmien ripottelua... ? Tajusitkohan nyt ollenkaan laittamaani kuvaa ja sen ajatusta? Säikeisiinhän ei tässä ole ripoteltuna mitään mielivaltaisia osaongelmia. Kyseessähän on täysin lineaarinen suoritus putki, jossa säie toimii käytännössä suodattimena. Kaikki tapahtuu järjestyksessä, suorituksessa ei ole mitään satunnaista tai mielivaltaista.
"Mielivaltaisten osaongelmien ripottelulla" tarkoitin, että se, että ongelman saa jaettua jonkinlaisiin osaongelmiin, ei tarkoita sitä, että jokainen osaongelma kannattaisi lähtökohtaisesti ratkaista omassa säikeessä. Ja en tosiaan kuvastasi hahmottanut, miten säikeet ja kanavat CSP-algoritmissa toimivat, enkä yksinkertaisimman menetelmän titteliä antaisi tälle algorimille myöskään koodiesimerkin perusteella.
Mutta täytyy myöntää, että tässä tapauksessa tuo CSP-toteutus puskuroitujen kanavien kanssa toimii Infernon virtuaalikoneella yllättävän hyvin ja lähestymistapa on kieltämättä mielenkiintoinen. Toteutus tietenkin vaatii ehdottomasti tuollaisia kevyitä säikeitä käyttävän virtuaalikoneen. Esimerkiksi Windowsin tai Linuxin säikeillä vastaavaa ei kannata yrittää.
os kirjoitti:
Mutta täytyy myöntää, että tässä tapauksessa tuo CSP-toteutus puskuroitujen kanavien kanssa toimii Infernon virtuaalikoneella yllättävän hyvin ja lähestymistapa on kieltämättä mielenkiintoinen. Toteutus tietenkin vaatii ehdottomasti tuollaisia kevyitä säikeitä käyttävän virtuaalikoneen. Esimerkiksi Windowsin tai Linuxin säikeillä vastaavaa ei kannata yrittää.
Tulipas vaan tässä työn ohessa mieleen... Missään ei muuten sanota, että CSP-toteutuksen prosessien (säikeiden) pitäisi kaikkien toimia samalla koneella. Tämä kyseinen algoritmihan on luonnostaan sellainen, että se skaalautuu helposti useammalle koneelle (eihän tässä tarvita muuta, kuin tapa yhdistää eri koneilla toimivat prosessi putket).
Voisi olla mielenkiintoinen harjoitustehtävä...
Aihe on jo aika vanha, joten et voi enää vastata siihen.