Eratosthenovo síto v C - rychlá implementace

Před pár dny mě inspiroval program na výpis prvočísel od Satíka

Satíkova aplikace generuje poměrně rychle seznam všech prvočísel v řádu přibližně 10⁶. U vyšších hodnot už rychle ztrácí dech, což je celkem pochopitelné a očekávané.

Začal jsem se zajímat, jak rychle lze vygenerovat prvočísla v rozsahu 10⁹?

Napsal jsem tedy základní algoritmus, který využívá naivní implementaci Eratosthenova síta a začal dělat testy. Během pokusů jsem vyzkoušel následující implementace:

1. naive Eratosthenes sieve
2. segmented Eratosthenes sieve
3. wheel factorization
4. optimized Eratosthenes sieve

Naivní verze síta je podle očekávání rychlá pouze v rozsahu 10⁶.

Na druhou stranu rozsah 10⁹ je stále ještě směšně malý, takže optimalizace síta za využití segmentace je pouze "kanón na vrabce". I když segmentace přinesla poměrně slušné zrychlení, z hlediska celkového času nepovažuji zisk za odpovídající vůči ceně komplikovaného kódu.

Faktorizace přes "kolo" při velikosti kola 30 také zrychlila běh, ale ani ta nestačila optimalizované verzi algoritmu, pro kterou jsem se nakonec rozhodl zejména kvůli jednoduchosti implementace.

Pokud hledáme prvočísla v rozsahu 3 .. 2*n, připravíme si pole o velikosti n, kde prvek na pozici i reprezentuje hodnotu 2*n+1. Pole inicializujeme na true. Následující kód označí všechna složená čísla v daném rozsahu, zatímco prvočíselné hodnoty zůstanou nezměněny.

// iprime je 64-bitové celé číslo bez znaménka
for (iprime i = 1, j = 3; i < n; i++, j += 2)
    for (iprime k = j * j >> 1; k < n && primes[i]; k += j)
        primes[k] = 0;

Poté platí, že liché číslo x je prvočíslo, pokud primes[x >> 1] != 0.

Na procesoru Core i5 byly časy běhu algoritmu přibližně následující (nezahrnuje generování výstupního souboru):

• 10⁶ - 0.002s
• 10⁷ - 0.03s
• 10⁸ - 0.51s
• 10⁹ - 5.9s
• 10⁹ - 1.5s při segmentaci síta

Jak rychle lze tedy najít prvočísla až do miliardy? Během několika sekund.

Využitím řady dalších optimalizací je možné získat prvočísla v rozsahu 10⁹ v čase okolo jedné sekundy. Tyto úpravy zahrnují zejména segmentaci síta kombinovanou s úsporným paměťovým modelem spadajícím do L1 a L2 CPU cache. To se ovšem vyplatí až u vyšších hodnot okolo řádu 10¹². Na mém kódu se mi líbí jednoduchost spojená se dostačujícím výkonem.

Zdrojový kód v jazyce C je v příloze.

Galerie

Program byl vytvořen v roce 2014.