Простые числа. Решето Эратосфена

Это статья посвящена простым числам и эффективным способам их вычисления. Сразу скажу, что те алгоритмы, которые тут приведены, являются весьма и весьма нетривиальными и самому мне не давались довольно долго, но в итоге я их всё-таки придумал, написал и спешу поделиться со всеми вами. Исходные тексты в статье будут приведены на языках C# и Haskell.

Простое число – это натуральное число больше единицы, которое имеет ровно два делителя: единицу и само это число.

Решето Эратосфена – древний, но при этом весьма эффективный и до сих пор широко используемый алгоритм поиска всех простых чисел не превосходящих некоторого N.

Запишем подряд все числа от 2 до N. Дальше вычеркнем из этого списка все числа кратные 2, исключая саму двойку, потом вычеркнем все числа кратные 3, исключая само число 3, число 4 уже вычеркнуто, вычеркиваем числа кратные 5 и т.д. Продолжаем этот процесс, пока квадрат очередного числа не превысит N.

Самая простая программная реализация этого алгоритма выглядит следующим образом:

bool[] table = new bool[100]; 
int i,j; 
// Отмечаем все числа как простые 
for ( i = 0; i < table.Length; i++ ) 
    table[i] = true; 
// Вычеркиваем лишнее 
for ( i = 2; i*i < table.Length; i++ ) 
    if ( table[i] ) 
        for ( j = 2*i; j < table.Length; j += i ) 
            table[j] = false; 
// Выводим найденное 
for ( i = 2; i < table.Length; i++ ) 
    if ( table[i] ) 
        Console.WriteLine(i); 

Итак, программка очень простая, но что в ней можно улучшить? Приведённая программа выведет список всех простых чисел из первой сотни. Если необходимо вычислить все числа до 1000 достаточно просто изменить размер массива table. Однако в C# объем памяти, занимаемый переменной типа bool равен одному байту. Следовательно, для нахождения всех чисел до одного миллиарда нам потребуется памяти чуть меньше чем гигабайт. А нам ведь хочется найти очень много простых чисел, не только до жалкого миллиарда, правда? :)

Первая и очевидная оптимизация – это использовать для хранения информации об одном числе одного бита. В этом случае в один гигабайт войдут все числа не превосходящие 8`589`934`592. Для работы с такими числами уже не подходят 32-х битные переменные. Очень хорошо! Переходим на 64-х битные. Сама программа:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66

using System; 
using System.Collections.Generic; 
using System.Text; 
 
namespace sieve 
{ 
    class Program 
    { 
        static uint[] table; 
 
        static void reset_bit(long bit_number) 
        { 
            long i = bit_number / 32; 
            long p = bit_number % 32; 
            table[i] = table[i] & (uint)(~(1 << (int)p)); 
        } 
 
        static bool get_bit(long bit_number) 
        { 
            long i = bit_number / 32; 
            long p = bit_number % 32; 
            if ( (table[i] & (uint)(1 << (int)p)) != 0 ) 
                return true; 
            return false; 
        } 
 
        static void Main(string[] args) 
        { 
            DateTime dt = DateTime.Now; 
            // One Gb! 
            // table = new uint[1024*1024*1024 / 4]; 
 
            table = new uint[1024*1024*10 / 4]; 
            long bit_count = (long)table.Length * (long)32; 
 
            long i,j; 
 
            for ( i = 0; i < table.Length; i++ ) 
                table[i] = 0xffffffff; 
 
            for ( i = 2; i*i < bit_count; i++ ) 
                if ( get_bit(i) ) 
                    for ( j = i*2; j < bit_count; j += i ) 
                        reset_bit(j); 
 
            TimeSpan calc_time = DateTime.Now - dt; 
            dt = DateTime.Now; 
 
            long prime_count = 0; 
            long max = 0; 
 
            for ( i = 2; i < bit_count; i++ ) 
                if ( get_bit(i) ) 
                { 
                    //Console.WriteLine(i); 
                    prime_count++; 
                    max = i; 
                } 
 
            Console.WriteLine("Всего найдено: {0}", prime_count); 
            Console.WriteLine("Максимальное: {0}", max); 
            Console.WriteLine("Время на вычисление: {0}", calc_time); 
            Console.WriteLine("Время на вывод: {0}", DateTime.Now - dt); 
        } 
    } 
} 

Кому хватит для жизни, достигнутых восьми с половиной миллиардов, могут дальше не читать. А со всеми кому мало, переходим к самому интересному. Если внимательно посмотреть, можно увидеть, что сразу после первой итерации оказываются вычеркнутыми все четные числа. Таким образом, буквально после первых нескольких итераций будет вычеркнуто больше половины массива. Существует ли способ не хранить эти числа в массиве, тем самым существенно сэкономив память? Я не буду приводить программу, хранящую только нечетные числа, я сразу обобщу задачу следующим образом: реализовать алгоритм решета, таким образом, чтобы массив не включал чисел, кратных некоторому количеству первых простых (например: 2,3,5,7).

Для исследования разных закономерностей распределения чисел, не делящихся на первые несколько простых, я использовал язык Haskell. Здесь я приведу только ту цепочку рассуждений, и покажу только тот код, которые привели к верному решению. На самом деле я предпринял много попыток.

Итак, запускаем hugs или ghci и начинаем колдовать :) .

Для начала нам нужно определить список чисел, из которого исключены кратные некоторому множеству чисел:

nums a b = [x | x <- [a..], not (isDiv x b)]  
    where  
    isDiv num lst = 
        case lst of 
        []     -> False 
        (h:t)   -> if num `mod` h == 0 then True else isDiv num t 
 
> take 10 (nums 3 [2]) 
[3,5,7,9,11,13,15,17,19,21] 
> take 10 (nums 5 [2,3]) 
[5,7,11,13,17,19,23,25,29,31] 

Для удобства я определю несколько вариантов, с которыми и буду работать:

sieve_2_3 = nums 5   [2,3] 
sieve_2_5 = nums 7   [2,3,5] 
sieve_2_7 = nums 11  [2,3,5,7] 

Нам понадобиться работать с разностями между ближайшими числами в последовательностях. Функция преобразования некоторой последовательности в список разностей между её соседними элементами:

rzn lst = [b - a | (a,b) <- zip lst (tail lst)] 
 
> take 10 (rzn $ nums 3 [2]) 
[2,2,2,2,2,2,2,2,2,2] 
> take 10 (rzn sieve_2_3) 
[2,4,2,4,2,4,2,4,2,4] 
> take 10 (rzn sieve_2_5) 
[4,2,4,2,4,6,2,6,4,2] 

Ага! Уже что-то есть, можно заметить периодически повторяющиеся последовательности. Пишем функцию поиска повторяющихся последовательностей:

get_period lst = 
    let 
        is_all_eq lst = 
            let  
                fs = head lst 
            in 
                foldl1 (&&) (map (== fs) (tail lst)) 
        slice_period p lst = [take p lst] ++ (slice_period p (drop p lst)) 
        find_period p = 
            if is_all_eq $ take 5 (slice_period p lst) then p 
                else find_period (p+1) 
        found_period = find_period 1 
    in 
        take found_period lst 

Я не буду на ней подробно останавливаться. Просто предположим, что она есть и корректно работает :) . Кому интересно, я думаю, смогут в ней разобраться. Проверяем с помощью этой функции наши последовательности:

> get_period $ rzn $ nums 3 [2] 
[2] 
> get_period $ rzn $ sieve_2_3 
[2,4] 
> get_period $ rzn $ sieve_2_5 
[4,2,4,2,4,6,2,6] 
> get_period $ rzn $ sieve_2_7 
[2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4, 
8,6,4,6,2,4,6,2,6,6,4,2,4,6,2,6,4,2,4,2,10,2,10] 

Теперь разберемся, в какую сторону нам предстоит двигаться дальше. Собственно для полного счастья нам необходимо научиться делать две взаимообратные вещи: во-первых, по номеру числа в массиве определять само это число (ведь число в массиве это всего один бит), а во-вторых, по числу определять его позицию в массиве.

Самое время перейти к конкретным примерам. Будем разбираться с массивом без чисел кратных двум и трем.

Первые 10 чисел этого списка выглядят так:

5,7,11,13,17,19,23,25,29,31 

Его период:

2,4 

Следовательно, этот массив состоит из подряд идущих пар, числа внутри пары отличаются на 2, разница между первыми числами подряд идущих пар равна 6. При этом не забываем, что элемент с индексом 0 равен 5.

Имеем функцию нахождения индекса по числу:

find_idx_by_num_2_3 num =  
    let  
        a = num-5 
        d = a `div` 6 
        m = a `mod` 6 
    in 
        d*2 + (if m < 2 then 0 else 1) 

Обратите внимание, мы проверяем остаток от деления на 6, если он больше 2 (разницы между числами в паре), то, значит, мы имеем дело со вторым числом в паре и, следовательно, должны к индексу прибавить единицу.

> find_idx_by_num_2_3 5 
0 
> find_idx_by_num_2_3 7 
1 
> find_idx_by_num_2_3 25 
7 

Обратная функция. Нахождение числа по индексу:

find_num_by_idx_2_3 i = 
    let 
        a = i `div` 2 
        b = i `mod` 2 
    in 
        a * 6 + (if b == 1 then 2 else 0) + 5 

Тут мы уже делим и находим остаток от деления аргумента на 2 – длину периода для нашего массива.

> find_num_by_idx_2_3 0 
5 
> find_num_by_idx_2_3 1 
7 
> find_num_by_idx_2_3 7 
25 

Ну что ж, в частном случае задача решена. Осталась самая малость – обобщить полученный результат до списка с произвольным числом вычеркнутых первых простых.

Индекс по числу:

find_idx_by_num num period first_number = 
    let  
        period_sum = foldl1 (+) period 
        period_len = integer_len period 
        a = num - first_number 
        d = a `div` period_sum 
        m = a `mod` period_sum 
        find_inc lst sum len n = 
            case lst of 
            []     -> len 
            (h:t)     -> 
                if n < sum then len 
                else find_inc t (sum+h) (len+1) n 
    in 
        d*period_len + (find_inc period 0 0 m) - 1 

Здесь наибольший интерес представляет функция find_inc, она находит положение числа внутри периода: суммируем поочередно числа из периода, пока полученное число не станет больше чем искомое число (остаток от деления аргумента num, функции find_idx_by_num на сумму всех чисел периода).

Число по индексу:

find_num_by_idx i period first_number = 
    let 
        period_sum = foldl1 (+) period 
        period_len = integer_len period 
        a = i `div` period_len 
        b = i `mod` period_len 
        sm n list 
            | n == 0    = 0 
            | otherwise    = 
                case list of 
                []    -> 0 
                (h:t)    -> h + sm (n-1) t 
    in 
        a * period_sum + (sm b period) + first_number 

Здесь функция sm суммирует первые n чисел списка list.

Самые глазастые заметили использование некой неописанной функции integer_len, исправляюсь:

integer_len l = 
    case l of 
    []    -> 0 
    (h:t)    -> 1 + integer_len t 

Это обычный length, отличающийся от стандартного типом возвращаемого значения, в данном случае это Integer вместо Int. Haskell язык со строгой типизацией, и вольностей между Int и Integer он не позволяет.

Проверяем, проверяем, проверяем:

> take 10 sieve_2_7 
[11,13,17,19,23,29,31,37,41,43] 
> find_idx_by_num 19 (get_period $ rzn $ sieve_2_7) 11 
3 
> find_num_by_idx 3 (get_period $ rzn $ sieve_2_7) 11 
19 

Ну что ж, теоретическую часть на этом можно считать законченной, способ работы с нужным форматом массива найден, с чем я себя и вас поздравляю! Ура! :)

Переходим к практической реализации всего вышенаколдованного.

Сразу приходит в голову примерно такая реализация: берем очередной ненулевой бит нашего массива, по его номеру находим число, которому он соответствует, дальше начинаем перебирать все кратные ему числа и вычеркивать их из массива с помощью функции поиска номера числа по самому этому числу. Звучит вроде не плохо и так уже можно сделать работающую реализацию, однако нужно учесть, что каждый поиск числа по индексу и поиск индекса по числу содержит внутри себя перебор периода для нашего формата, хотелось бы от этой операции избавиться.

Снова уходим в глубину Haskell’я (последний раз, клянусь :) ). Посмотрим, как распределены числа кратные некоторому N в наших массивах. Функция находит индексы всех чисел кратных num в списке list:

div_idx list num = [b | (a,b) <- zip list [0..], a `mod` num == 0] 

К примеру, индексы чисел кратных 5 в массиве без двоек и троек:

> take 10 (div_idx sieve_2_3 5) 
[0,7,10,17,20,27,30,37,40,47] 

Ну что, знакомая картина, явно наша любимая периодически повторяющаяся последовательность разностей между соседними числами:

> get_period $ rzn $ div_idx sieve_2_3 5 
[7,3] 

Конечный алгоритм будет выглядеть следующим образом: берем очередной единичный бит, находим соответствующее ему число. Находим N+1 ближайших кратных ему чисел (при этом исключаем числа кратные простым из списка вычеркнутых) и находим разности между их индексами. Дальше начинаем последовательно обнулять индексы всех кратных найденному чисел. К примеру, мы обнаружили единичный бит с номером 0 соответствующий числу 5 в массиве без 2 и 3, мы знаем, что индексы чисел кратных 5 изменяются с периодом [7,3] (см. выше), соответственно, обнуляем биты с номерами 7,10,17,20,27 и т.д.

В этом месте следовало бы сказать, что и это ещё не всё, однако этого (пока?) не произойдет. Дальнейший путь ускорения я вижу в выявлении математического закона, формирующего период разностей индексов чисел кратных некоторому N. Господа математики, если у вас есть какие-то идеи на этот счет, милости прошу писать комментарии тут или мне на e-mail.

Я не буду тут приводить полный вариант финальной программы на C#, остановлюсь разве только на C# реализации функций получения индекса по числу и числа по индексу:

public Int64 get_index_by_number(Int64 num) 
{ 
    Int64 a = num - first_number; 
    Int64 d = a / period_sum; 
    Int64 m = a % period_sum; 
 
    int sum = 0; 
    int i; 
 
    for ( i = 0; i < period.Length; i++ ) 
    { 
        if ( m < sum ) break; 
        sum += period[i]; 
    } 
     
    i--; 
 
    return d*period.Length + i; 
} 
 
public Int64 get_number_by_index(Int64 idx) 
{ 
    Int64 a = idx / period.Length; 
    Int64 b = idx % period.Length; 
 
    Int64 sum = 0; 
    for ( int i = 0; i < b; i++ ) 
        sum += period[i]; 
 
    return a * period_sum + sum + first_number; 
} 

Здесь period, period_sum и first_number определены в файле gen.cs, этот файл генерируеться программой на Haskell, один из вариантов этого файла:

namespace PrimeNumbers 
{ 
    partial class PNTable 
    { 
        protected static readonly int[] period = {4,2,4,2,4,6,2,6}; 
        protected static readonly int[] base_primes = {2,3,5}; 
        protected const long period_sum = 30; 
        protected const long first_number = 7; 
    } 
} 

Генератор для таких файлов находится в папке haskell внутри приложенного к статье файла. Он получает из командной строки число первых простых, которые должны быть исключены и выводит в стандартный поток вывода сгенерированный C# код.

Теперь посмотрим в цифрах, чего мы всё-таки добились:

Классический вариант, для хранения числа используется один бит. Размер массива один гигабайт:

Всего найдено простых чисел: 393`615`806 
Максимальное простое число: 8`589`934`583 
Время на вычисление: 00:41:44 

Вариант решета с исключением чисел кратных 2,3,5,7,11,13. Размер массива один гигабайт:

Всего найдено простых чисел: 1`907`461`679 
Максимальное простое число: 44`783`981`911 
Время на вычисление: 00:58:58 

С помощью улучшенного алгоритма удалось найти почти в 5 раз больше простых чисел, при этом затраты по времени увеличились только в полтора раза. Очень хороший результат!

Я должен заметить, что дальше исключения последовательности 2,3,5,7,11,13 продвигаться нет смысла. Во-первых затраты времени на вычисление массива с разностями индексов для каждого очередного числа становятся очень большими. Во-вторых эффект от каждого очередного числа сильно снижается. Рассмотрим, как меняется диапазон чисел доступных для работы. Будем брать сотый член разных последовательностей:

> (nums 1 []) !! 100 
101 
> (nums 3 [2]) !! 100 
203 
> (nums 5 [2,3]) !! 100 
305 
> (nums 7 [2,3,5]) !! 100 
379 
> (nums 11 [2,3,5,7]) !! 100 
443 
> (nums 13 [2,3,5,7,11]) !! 100 
491 
> (nums 17 [2,3,5,7,11,13]) !! 100 
529 
> (nums 19 [2,3,5,7,11,13,17]) !! 100 
569 
> (nums 23 [2,3,5,7,11,13,17,19]) !! 100 
593 
> (nums 29 [2,3,5,7,11,13,17,19,23]) !! 100 
601 

Видно, что наибольший эффект имеет вычеркивание двойки, чуть худший – тройки и чем дальше, тем меньший эффект от вычеркивания очередного простого числа.

Также надо сказать, что использовать массив с вычеркнутыми 2,3,5,7,11,13 можно только на больших объемах памяти, по крайней мере в моей реализации, иначе не будет хватать чисел для вычисления массива с разностями для очередного числа и этот массив будет заполнен некорректно.

Мне было очень интересно решать эту задачу, надеюсь, вам также было интересно это все читать. Я постарался показать на примере очень интересное свойство функциональных языков: код живет буквально на кончиках пальцев, программу можно крутить и разворачивать как угодно, можно очень быстро поверять свои идеи и смотреть на возникающие закономерности. И, кроме всего, приятно чувствовать себя осилившим долго не дававшуюся сложную задачу :) .

Полные исходные тексты всех приведенных в статье программ находятся в прилагаемом файле.