Волшебное решето Эратосфена

Наверняка все, кто читает этот пост не раз использовали, или хотя бы слышали о решете Эратосфена — методе отыскания простых чисел. Сама проблема получения простых чисел занимает ключевое место в математике, на ней основаны некоторые криптографические алгоритмы, например RSA. Есть довольно много подходов к данной задаче, но в этой статье я остановлюсь на некоторых модификациях самого простого из них — решета Эратосфена.
Принцип решета прост: пускай нам нужно отыскать простые числа в промежутке от единицы до некоторого N <= 10^6. Мы заводим массив на N элементов и заполняем его true. Затем последовательно проходим по нему до корня из N, и встречая true, вычеркиваем все числа с этим шагом до N. Алгоритм выглядит компактно и просто, привожу его на языке java.

1. Arrays.fill(isPrime,true);
2. isPrime[1] = false;
3. for (int i=2; i*i < N; i++)
4.    if (isPrime[i])
5.       for (int j=i*i; j < N; j+=i)
6.          isPrime[j] = false;
7.  

Алгоритм работает за O(N*log(log(N))), поэтому для небольших чисел он вполне подходит. Рассмотрим случай, когда нам необходимо найти простые числа не от одного до N, а от n до m. Пускай мы имеем ограничения 1 <= m <= n <= 10^9; n-m <= 10^6. Здесь нам необходимо применить алгоритм, называемый двойным решетом. Он заключается в том чтобы найти простые числа до корня из n, затем сохранить их в отдельный массив и точно так же «вычеркивать» эти простые числа с определенным шагом, но уже из необходимого нам промежутка [m, n]. Кратко этот алгоритм будет выглядеть так:

1. int primes = new int[P]; // заполняем его простыми до корня из n  
2. boolean sieve = new boolean[n-m+1]; // вторичное решето  
3. Arrays.fill(sieve, true);
4. for (int i= 0; i<P; i++) {
5.    int h = m % primes[i];
6.    int j = h ==  0 ?  0 : primes[i] - h;
7.    for (; j<=n-m; j+=primes[i])
8.       sieve[j] = false;
9.  

Так можно справляться с диапазонами, достаточно удаленными от нуля. Теперь подходим к самому интересному моменту: а что, если нам необходимо вывести простые числа от 0 до N = 10^8? Если вспомнить асимптотику, то это на первый взгляд кажется реальным даже для обычного решета, но посмотрев на затраты памяти: 10^8 * 4 = 400МБ мы видим что задача не такая тривиальная, особенно, если у нас есть жесткие временные ограничения. Для решения этой задачи можно выделить два подхода, а именно:

• упаковка решета
• учитывание кэш-памяти

Рассмотрим их немного подробнее. Упаковка заключается вот в чем: в современных языках программирования тип boolean занимает в среднем от одного до четырех байт, хотя его можно закодировать всего одним битом. Поэтому мы можем создавать массив целых чисел и работать с каждым из них как с побитовой маской, тем самым экономя память в 8-30 раз. Преимущества видны сразу, но теперь остановимся чуть подробнее на второй идее, которая гарантированно ускоряет решето в несколько раз
Многим известно, что в нашем компьютере между процессорной и оперативной памятью находится кэш-память, работа с ней проводится намного быстрее, чем с оперативной, но ее размеры ограничены. Например, при работе с большим массивом, процессор загружает в кэш некоторую его часть, работает с ней, потом переносит обратно в оперативную, загружает другую и так далее. А теперь вспомним наш алгоритм решета: каждое простое число мы вычеркивали из всего массива, проходясь по нему от начала до конца. Поэтому процессор много раз будет загружать в кэш разные отрезки массива и скорость на этом будет сильно теряться. Данный подход предлагает минимизировать затраты на копирование массива из одной памяти в другую. Это сделать несложно если весь наш промежуток разделить на кусочки, до 3*10^4 элементов, что приблизительно равно размеру кэша и работать с ними по-порядку. Тогда мы за минимальное количество загрузок разберемся со всем массивом. Примерно так это будет выглядеть:

1. int CACHE = 30000; // размер кэша  
2. int M = (int)Math.sqrt(N)+1;
3.  
4. int primes = new int[P]; // массив простых чисел до корня из N  
5. boolean segment = new boolean[CACHE]; // вторичное решето  
6. for (int I=M-1; I < N; I+=CACHE) {
7.         Arrays.fill(segment, true);
8.         for (int i= 0; i < P; i++) {
9.             int h = I % primes[i];
10.             int j = h >  0 ? primes[i] - h :  0;
11.             for (; j < CACHE; j+=primes[i])
12.                 segment[j] = false;
13.         }
14.         for (int i= 0; i<CACHE; i++) {
15.             if (segment[i] && (i + I < N)) {
16.                 out.println(i+I); // выводим простое число на экран  
17.             }
18.         }
19. }

Используя этот метод мы очень сильно ускоряем наше решето, практически не усложняя его, для многих задач это ощутимо помогает, например, одна из них: TDPRIMES, на этой задаче можно потренироваться и хорошо увидеть что обычное решето не укладывается в 10с, а сегментированное проходит за 2.8.