Ищем простые числа до триллиона за тридцать минут10.11.2020 10:47

image

Поиск простых чисел — популярная задача среди программистов, увлекающихся математикой. Самый известный алгоритм, придуманный, по-видимому, больше двух тысяч лет назад, — решето Эратосфена; в настоящее время существует бесчисленное множество его вариантов и оптимизаций.
Сегодня я хотел бы поделиться с вами различными вариантами реализации поиска простых чисел на языке C#, начиная с классических алгоритмов — решета Эратосфена, Сундарама и Аткина, и кончая различными оптимизациями (сегментация, факторизация). Особый упор я делал на простоту: самый быстрый из алгоритмов, который мне удалось получить, содержит 120 строк кода и ищет простые числа до триллиона меньше, чем за 30 минут, а до миллиарда — меньше, чем за секунду (это далеко от производительности лучших из существующих библиотек по поиску простых чисел, но эти библиотеки обычно содержат свыше 4000 строк кода).
В заключение мы применим самую быструю реализацию для поиска максимального расстояния между двумя соседними простыми числами до триллиона. Прежде чем заходить под кат, я предлагаю вам попытаться угадать ответ. Для сравнения, для простых чисел до 100 максимальное растояние равно 8 (между соседними простыми числами 89 и 97), а до тысячи — 20 (между 887 и 907).
Весь исходный код можно найти на гитхабе.

1. Решето Эратосфена


Решето Эратосфена — самый простой и известный алгоритм для поиска простых чисел. Его суть состоит в следующем:
— Выписываем числа 2, 3, …, N;
— Берем первое число в списке $p=2$ и вычеркиваем из списка все числа, кратные ему, начиная с $p^2$: 4, 6, 8, …
— Берем следующее число в списке $p=3$ и вычеркиваем кратные ему 9, 15, …
— Продолжаем процедуру для каждого последующего числа из оставшихся в списке $p=5, 7, 11, ...,$ пока выполняется неравенство $p^2 <= \sqrt{N}$; оставшиеся после этого в списке числа и будут простыми:
вычеркиваем кратные 2: [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, ...]
вычеркиваем кратные 3: [2, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, ...]
вычеркиваем кратные 5: [2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, 35, ...]
и так далее...
Практическая реализация использует битовый массив длины N, с всеми битами изначально установленными в 1, и для каждого составного числа меняет соответствующий бит на 0.

SieveOfEratosthenes.cs 
public class SieveOfEratosthenes : ISieve
{
    private BitArray Data;
    public int Length => Data.Length;

    public SieveOfEratosthenes(int length)
    {
        Data = new BitArray(length);
        Data.SetAll(true);

        for (int p = 2; p * p < length; p++)
        {
            if (Data[p])
            {
                for (int i = p * p; i < Length; i += p)
                {
                    Data[i] = false;
                }
            }
        }
    }

    public void ListPrimes(Action callback)
    {
        for (int i = 2; i < Length; i++)
        {
            if (Data[i]) callback.Invoke(i);
        }
    }
}


Поскольку наша задача — сравнить между собой разные реализации поиска простых чисел, мы наследуем все варианты от интерфейса ISieve, имеющего только один метод: 

public interface ISieve
{
    void ListPrimes(Action callback);
}


Для перечисления простых чисел был выбран способ колбеков, а не более привычные для C# энумераторы, поскольку он на 10–40% быстрее. Этот интерфейс позволяет единообразно тестировать разные реализации и сравнивать их между собой; в частности, мы будем считать количество, сумму и хеш простых чисел до миллиарда:

PerformanceTests 
static void TestSieve(long length, 
                         long expectedCount,
                         long expectedSum,
                         long expectedHash)
{
    var sw = Stopwatch.StartNew();
    ISieve sieve = CreateSieve(length);
    long count = 0, sum = 0, hash = 0;
    Action action = (p) => { count++; sum += p; hash = hash * 31 + p; };
    sieve.ListPrimes(action);
    AssertEquals(expectedCount, count);
    AssertEquals(expectedSum, sum);
    AssertEquals(expectedHash, hash);
    Console.WriteLine($"{typeof(T).Name} up to {length:N0} in {sw.Elapsed}");
}

static void TestSieve_1G() =>
    TestSieve(1_000_000_000, 50847534, 24739512092254535, 1480810120364005255);


Алгоритм: решето Эратосфена
Ищет простые числа до миллиарда: 12.6 секунд.

2. Решето Сундарама


Решето Сундарама — второй по известности алгоритм поиска нечетных простых чисел, придуманный индийским студентом Сундарамом в 1934 году. Его суть состоит в вычеркивании из списка чисел от 1 до $(N-1)/2$ всех чисел вида

$i+j+2ij$

$1<=i<=j$

После этого, для оставшихся в списке чисел $n_i$, последовательность чисел $2n_i+1$ будет представлять все нечетные простые числа до N.
вычеркиваем 1+j+2*1*j: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, ...]
вычеркиваем 2+j+2*2*j: [1, 2, 3, 5, 6, 8, 9, 11, 12, 14, 15, ...]
вычеркиваем 3+j+2*3*j: [1, 2, 3, 5, 6, 8, 9, 10, 11, 14, ..., 24, ...]
и так далее...
последовательность 2n+1: [3, 5, 7, 11, 13, 17, 19, ...]

SieveOfSundaram.cs 
public class SieveOfSundaram : ISieve
{
    private BitArray Data;
    public int Length { get; private set; }

    public SieveOfSundaram(int length)
    {
        Length = length;
        Data = new BitArray((Length + 1) / 2);
        Data.SetAll(true);

        for (int i = 1; i + i + 2 * i * i < Data.Length; i++)
        {
            for (int j = i; i + j + 2 * i * j < Data.Length; j++)
            {
                Data[i + j + 2 * i * j] = false;
            }
        }
    }

    public void ListPrimes(Action callback)
    {
        callback.Invoke(2);

        for (int i = 1; i < Data.Length; i++)
        {
            if (Data[i])
            {
                int p = i * 2 + 1;
                if (p >= Length) break;
                callback.Invoke(p);
            }
        }
    }
}


Заметим, что, в отличие от решета Эратосфена, алгоритм вычеркивает числа для любых комбинаций $(i, j)$, поэтому его теоретическая сложность хуже: $O(N logN)$ против $O(N log log N)$. И действительно, на нашем тесте он работает медленней:

Алгоритм: решето Сундарама
Ищет простые числа до миллиарда: 18 секунд.

3. Решето Аткина


В 2003 году два профессора математики из Иллинойского университета в Чикаго Артур Аткин и Даниэль Бернштейн опубликовали статью, описывающую принципиально новый алгоритм поиска простых чисел, опирающийся на методы алгебраической теории чисел. Этот алгоритм теперь известен как решето Аткина. Его основная идея состоит в том, что существуют квадратичные формы $ax^2 + by^2 = n$ для определенных классов $n$, имеющие нечетное количество решений только когда $n$ — либо простое, либо имеет делитель $p^2$. Это позволяет с помощью квадратичных форм отсеить все составные числа, не делящиеся на квадрат простого числа, а потом способом, аналогичным решету Эратосфена, отсеить кратные квадратам простых чисел.
Первый шаг имеет изящную реализацию на битовых массивах: значение квадратичной формы вычисляется для всех подходящих комбинаций $(x, y)$, и соответствующий бит меняется на противоположный. Таким образом, если изначально все биты были установлены в 0, то в конце будут единичными только биты, переключенные нечетное количество раз, то есть биты для тех чисел $n$, для которых соответствующая квадратичная форма имеет нечетное количество решений!
Практическая реализация рассматривает числа $n$, не делящиеся на 2, 3, 5, и имеющие определенные остатки от деления на 60.

SieveOfAtkin.cs 
public class SieveOfAtkin : ISieve
{
    private BitArray Data;
    public int Length { get; private set; }

    public SieveOfAtkin(int length)
    {
        Length = length;
        Data = new BitArray(Length);
        Data.SetAll(false);

        SieveQuadForm1(Data);
        SieveQuadForm2(Data);
        SieveQuadForm3(Data);

        for (int p = 7; p * p <= length; p++)
        {
            if (Data[p])
            {
                for (int i = p * p; i < Length; i += p * p)
                {
                    Data[i] = false;
                }
            }
        }

    }

    private void SieveQuadForm1(BitArray data)
    {
        for (int x = 1; 4 * x * x < data.Length; x++)
        {
            for (int y = 1; 4 * x * x + y * y < data.Length; y += 2)
            {
                int n = 4 * x * x + y * y;
                int rem = n % 60;
                if (rem == 1 || rem == 13 || rem == 17 || rem == 29 || 
                    rem == 37 || rem == 41 || rem == 49 || rem == 53)
                {
                    data[n] = !data[n];
                }
            }
        }
    }

    private void SieveQuadForm2(BitArray data)
    {
        for (int x = 1; 3 * x * x < data.Length; x += 2)
        {
            for (int y = 2; 3 * x * x + y * y < data.Length; y += 2)
            {
                int n = 3 * x * x + y * y;
                int rem = n % 60;
                if (rem == 7 || rem == 19 || rem == 31 || rem == 43)
                {
                    data[n] = !data[n];
                }
            }
        }
    }

    private void SieveQuadForm3(BitArray data)
    {
        for (int x = 1; 2 * x * x < data.Length; x++)
        {
            for (int y = x - 1; y > 0; y -= 2)
            {
                int n = 3 * x * x - y * y;
                if (n >= data.Length) continue;
                int rem = n % 60;
                if (rem == 11 || rem == 23 || rem == 47 || rem == 59)
                {
                    data[n] = !data[n];
                }
            }
        }
    }

    public void ListPrimes(Action callback)
    {
        callback.Invoke(2);
        callback.Invoke(3);
        callback.Invoke(5);
        for (int i = 7; i < Length; i++)
        {
            if (Data[i])
            {
                callback.Invoke(i);
            }
        }
    }
}


Этот алгоритм асимптотически быстрее, чем решето Эратосфена: $O(N)$ против $O(N log log N)$, но за счет большей сложности в практическом применении часто проигрывает решету Эратосфена при сравнимом количестве оптимизаций.

Алгоритм: решето Аткина
Ищет простые числа до миллиарда: 11.4 секунд.

4. Колёсная факторизация


Колёсная факторизация — вид оптимизации решета Эратосфена, позволяющий заранее исключить из рассмотрения числа, заведомо не являющиеся простыми. Например, любое простое число, большее 6, имеет остаток от деления на 6 либо 1, либо 5: $p = 6k+1$ или $p=6k+5$: при любых других остатках число будет делиться либо на 2, либо на 3. Аналогично, простое число, большее $2*3*5=30$, может иметь только один из 8 остатков от деления на 30: $(1, 7, 11, 13, 17, 19, 23, 29)$. Еще одно преимущество оптимизации — уменьшение потребения памяти: вместо одного бита на число для решета Эратосфена или 0.5 бит на число для решета Сундарама, факторизация с базой $(2, 3, 5)$ требует примерно 0.27 бит на число, поскольку позволяет хранить 30 чисел в 8 битах. При аккуратном применении оптимизация позволяет ускорить просеивание в несколько раз, однако неоптимизированные варианты обычно работают медленно за счет большей сложности. В прошлой публикации я подробно описывал алгоритм, однако он не был оптимизирован по производительности, поэтому здесь я приведу более простой и эффективный вариант.

Wheel235.cs 
public class Wheel235 : ISieve
{
    private static int[] BIT_TO_INDEX = new int[] { 1, 7, 11, 13, 17, 19, 23, 29 };

    private static int[] INDEX_TO_BIT = new int[] {
            -1, 0,
            -1, -1, -1, -1, -1, 1,
            -1, -1, -1, 2,
            -1, 3,
            -1, -1, -1, 4,
            -1, 5,
            -1, -1, -1, 6,
            -1, -1, -1, -1, -1, 7,
        };

    private byte[] Data;
    public long Length { get; private set; }

    public Wheel235(long length)
    {
        Length = length;
        Data = new byte[(length + 29) / 30];
        for (long i = 0; i < Data.Length; i++) Data[i] = byte.MaxValue;

        for (long i = 7; i * i < Length; i++)
        {
            if (!IsPrime(i)) continue;
            for (long d = i * i; d < Length; d += i) ClearPrime(d);
        }
    }

    public bool IsPrime(long n)
    {
        if (n >= Length) throw new ArgumentException("Number too big");
        if (n <= 5) return n == 2 || n == 3 || n == 5;
        int bit = INDEX_TO_BIT[n % 30];
        if (bit < 0) return false;
        return (Data[n / 30] & (1 << bit)) != 0;
    }

    private void ClearPrime(long n)
    {
        int bit = INDEX_TO_BIT[n % 30];
        if (bit < 0) return;
        Data[n / 30] &= (byte)~(1 << bit);
    }

    public void ListPrimes(Action callback)
    {
        callback.Invoke(2);
        callback.Invoke(3);
        callback.Invoke(5);

        for (long position = 0; position < Data.Length; position++)
        {
            for (int bit = 0; bit < 8; bit++)
            {
                if ((Data[position] & (1 << bit)) != 0)
                {
                    long p = position * 30 + BIT_TO_INDEX[bit];
                    if (p <= 5) continue;
                    if (p >= Length) return;
                    callback.Invoke(p);
                }
            }
        }
    }
}


Алгоритм: (2,3,5)-факторизация
Ищет простые числа до миллиарда: 9.6 секунд.

5. Оптимизированное решето Эратосфена


В оригинальном решете Эратосфена мы создавали битовый массив для всех чисел, хотя очевидно, что простые числа, кроме двойки, не могут быть четными. Поэтому мы оптимизируем хранение, создавая массив, в котором каждый бит будет соответствовать нечетным числам. Соответственно, нечетное число $n$ будет храниться в бите с номером $(n -1)/2$. Кроме того, при просеивании чисел, кратных нечетным простым $p$, мы можем пропускать четные числа $p^2 + p, p^2 + 3p, p^2+5p, ...$ и обрабатывать только нечетные числа $p^2 + 2p, p^2+4p, p^2+6p, ...$

OptimizedSieveOfEratosthenes.cs 
public class OptimizedSieveOfEratosthenes : ISieve
{
    private BitArray Data;
    public int Length { get; private set; }

    public OptimizedSieveOfEratosthenes(int length)
    {
        Length = length;
        Data = new BitArray(Length / 2 + 1);
        Data.SetAll(true);

        int maxFactor = (int)Math.Sqrt(Length);

        for (int p = 3; p <= maxFactor; p += 2)
        {
            if (Data[p / 2])
            {
                for (int i = p * p; i < Length; i += 2 * p)
                {
                    Data[i / 2] = false;
                }
            }
        }
    }

    public void ListPrimes(Action callback)
    {
        callback.Invoke(2);

        for (int i = 1; i < Length / 2; i++)
        {
            if (Data[i])
            {
                long p = i * 2 + 1;
                if (p >= Length) break;
                callback.Invoke(p);
            }
        }
    }
}


Эта оптимизация дает нам практически двукратный выигрыш в производительности, и, кроме того, требует 0.5 бита на натуральное число, как и решето Сундарама.

Алгоритм: Оптимизированное решето Эратосфена
Ищет простые числа до миллиарда: 6.6 секунд.

6. Оптимизированное решето Сундарама


Предыдущая реализация подозрительно похожа на решето Сундарама. По крайней мере, простые числа $p$ и там, и там хранятся в бите с номером $(p - 1)/2$. При этом решето Сундарама вычеркивает числа вида $n=i+j+2ij$. Можно заметить, что из этого следует

$2n+1=(2i+1)(2j+1),$

то есть, оно делает в точности то же самое, что и оптимизированное решето Эратосфена, только вычисления производятся над индексами битов, а не над самими числами. Поэтому сразу становится видно различие: решето Эратосфена вычеркивает только числа, кратные простым, а решето Сундарама — кратные всем числам! Очевидная оптимизация решета Сундарама — обрабатывать только те числа $i$, для которых $2i+1$ — простое число: 

    if (!Data[i]) continue;


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

OptimizedSieveOfSundaram.cs 
public class OptimizedSieveOfSundaram : ISieve
{
    private BitArray Data;
    public int Length { get; private set; }

    public OptimizedSieveOfSundaram(int length)
    {
        Length = length;
        Data = new BitArray((Length + 1) / 2);
        Data.SetAll(true);

        for (int i = 1; i + i + 2 * i * i < Data.Length; i++)
        {
            if (!Data[i]) continue;

            for (int j = i; i + j + 2 * i * j < Data.Length; j++)
            {
                Data[i + j + 2 * i * j] = false;
            }
        }
    }

    public void ListPrimes(Action callback)
    {
        callback.Invoke(2);

        for (int i = 1; i < Data.Length; i++)
        {
            if (Data[i])
            {
                int p = i * 2 + 1;
                if (p >= Length) break;
                callback.Invoke(p);
            }
        }
    }
}


Производительность, что не удивительно, совпадает с предыдущим алгоритмом:

Алгоритм: Оптимизированное решето Сундарама
Ищет простые числа до миллиарда: 6.5 секунд.

7. Сегментированное решето Эратосфена


Все предыдущие алгоритмы хранили весь результат в памяти, что ограничивает их применение небольшими диапазонами. Например, при поиске простых чисел до триллиона, решето Эратосфена будет требовать 120 гигабайт оперативной памяти; кроме того, его скорость будет деградировать, так как просеивание (вычеркивание составных чисел) бегает по всему массиву, не эффективно используя кеши процессора.
Чтобы обойти это ограничение, вспомним, что для нахождения простых чисел до $N$ нам достаточно вычеркнуть числа, кратные простым, не превосходящим $\sqrt{N}$. Это позволяет нам сначала найти простые числа до $\sqrt{N}$, а потом разбить отрезок $[2..N]$ на меньшие отрезки (сегменты), и в каждом вычеркнуть кратные простым по отдельности. Нахождение первого числа, кратного $p$, в сегменте требует немного арифметики: 

long first = p * p;
if (first < segmentStart)
{
    first += (segmentStart - first + p - 1) / p * p;
}


вычеркиваем кратные 2: [ 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90]
вычеркиваем кратные 3: [ 81, 83, 85, 87, 89]
вычеркиваем кратные 5: [83, 85, 89]
вычеркиваем кратные 7 (нет в интервале): [83, 89]
Этот подход называется сегментированное решето, и с его применением требования к памяти снижаются с $O(N)$ до $O(\sqrt{N})$. На практике размер интервалов выбирают либо равным $\sqrt{N}$, либо зависящим от размера кеша процессора.

SegmentedSieveOfEratosthenes.cs 
public class SegmentedSieveOfEratosthenes : ISieve
{
    private long Length;
    private int FirstChunkLength;
    private long[] FirstPrimes;

    public SegmentedSieveOfEratosthenes(long length)
    {
        Length = length;
        FirstChunkLength = (int)Math.Sqrt(length) + 1;
        SieveOfEratosthenes sieve = new SieveOfEratosthenes(FirstChunkLength);
        List firstPrimes = new List();
        sieve.ListPrimes(firstPrimes.Add);
        FirstPrimes = firstPrimes.ToArray();
    }

    private void SieveSegment(BitArray segmentData, long segmentStart, long segmentEnd)
    {
        segmentData.SetAll(true);

        foreach (long p in FirstPrimes)
        {
            long first = p * p;
            if (first < segmentStart)
            {
                first += (segmentStart - first + p - 1) / p * p;
            }

            for (long i = first; i < segmentEnd; i += p)
            {
                segmentData[(int)(i - segmentStart)] = false;
            }
        }
    }

    public void ListPrimes(Action callback)
    {
        foreach (long p in FirstPrimes) callback.Invoke(p);

        int segmentLength = FirstChunkLength;
        BitArray segmentData = new BitArray(segmentLength);
        long segmentStart = FirstChunkLength;
        long segmentEnd = Math.Min(segmentStart + segmentLength, Length);
        while (segmentStart < Length)
        {
            SieveSegment(segmentData, segmentStart, segmentEnd);
            for (long i = segmentStart; i < segmentEnd; i++)
            {
                if (segmentData[(int)(i - segmentStart)]) callback.Invoke(i);
            }
            segmentStart = segmentEnd;
            segmentEnd = Math.Min(segmentStart + segmentLength, Length);
        }
    }
}


За счет лучшего использования кеша алгоритм работает быстрее обычного решета Эратосфена на тех же данных:

Алгоритм: Сегментированное решето Эратосфена
Ищет простые числа до миллиарда: 8.6 секунд.

8. Оптимизированное сегментированное решето


Данная оптимизация относится к «народной» — она широко известна, но, похоже, не имеет своего названия. Цель оптимизации — избавиться от арифметического выражения, вычисляющего первое число в сегменте, кратное простому числу $p$. Так как сегментированное решето обрабатывает сегменты подряд в порядке возрастания, нетрудно убедиться, что нужное число будет равно первому числу, кратному $p$, выхожящему за граниы предыдущего сегмента. Таким образом, нам достаточно для каждого простого числа $p<=\sqrt{N}$ хранить текущее кратное ему число. В упрощенном виде это можно представить так: 

long[] PrimeMultiples = FirstPrimes.Select(p => p * p).ToArray();
...
while(PrimeMultiples[i] < segmentEnd)
{
    segmentData[PrimeMultiples[i] - segmentStart] = false;
    PrimeMultiples[i] += FirstPrimes[i];
}


Полная реализация с небольшими дополнительными оптимизациями позволяет ускорить предыдущий алгоритм примерно на 20%

OptimizedSegmentedSieve.cs 
public class OptimizedSegmentedSieve : ISieve
{
    const int BUFFER_LENGTH = 64 * 1024;
    private long Length;
    private long[] FirstPrimes;
    private long[] PrimeMultiples;

    public OptimizedSegmentedSieve(long length)
    {
        Length = length;
        int firstChunkLength = (int)Math.Sqrt(length) + 1;
        SieveOfEratosthenes sieve = new SieveOfEratosthenes(firstChunkLength);
        List firstPrimes = new List();
        sieve.ListPrimes(firstPrimes.Add);
        FirstPrimes = firstPrimes.ToArray();
        PrimeMultiples = FirstPrimes.Select(p => p * p).ToArray();
    }

    private void SieveSegment(BitArray segmentData, long segmentStart, long segmentEnd)
    {
        segmentData.SetAll(true);
        long segmentLength = segmentEnd - segmentStart;
        for (int i = 0; i < PrimeMultiples.Length; i++)
        {
            long current = PrimeMultiples[i] - segmentStart;
            long prime = FirstPrimes[i];
            if (current >= segmentLength) continue;

            while (current < segmentLength)
            {
                segmentData[(int)current] = false;
                current += prime;
            }
            PrimeMultiples[i] = segmentStart + current;
        }
    }

    public void ListPrimes(Action callback)
    {
        BitArray segmentData = new BitArray(BUFFER_LENGTH);
        long segmentStart = 2;
        long segmentEnd = Math.Min(segmentStart + BUFFER_LENGTH, Length);
        while (segmentStart < Length)
        {
            SieveSegment(segmentData, segmentStart, segmentEnd);
            for (int i = 0; i < segmentEnd - segmentStart; i++)
            {
                if (segmentData[i]) callback.Invoke(segmentStart + i);
            }
            segmentStart = segmentEnd;
            segmentEnd = Math.Min(segmentStart + BUFFER_LENGTH, Length);
        }
    }
}


Алгоритм: Оптимизированное сегментированное решето
Ищет простые числа до миллиарда: 7.1 секунд.

9. Сегментированная факторизация с базой 2


Раньше мы рассматривали оптимзацию, позволяющую хранить только нечетные числа, и ничего нам не мешает применить ее к предыдущему алгоритму. Только теперь массив PrimeMultiples будет хранить не $p^2 + kp$, а $((p^2 + 2kp) - 1) / 2$: 

long[] PrimeMultiples = FirstPrimes.Select(p => (p * p - 1) / 2).ToArray();
...
while(PrimeMultiples[i] < segmentEnd)
{
    segmentData[PrimeMultiples[i] - segmentStart] = false;
    PrimeMultiples[i] += FirstPrimes[i];
}


Интересно, что основной цикл исключения составных чисел не изменился, разница только в инициализации массива PrimeMultiples и в преобразовании индекса битов в число по знакомой нам формуле $p=2i+1$.

SegmentedWheel2.cs 
public class SegmentedWheel2 : ISieve
{
    const int BUFFER_LENGTH = 64 * 1024;
    private long Length;
    private long[] FirstPrimes;
    private long[] PrimeMultiples;

    public SegmentedWheel2(long length)
    {
        Length = length;
        int firstChunkLength = (int)Math.Sqrt(length) + 1;
        SieveOfEratosthenes sieve = new SieveOfEratosthenes(firstChunkLength);
        List firstPrimes = new List();
        sieve.ListPrimes(firstPrimes.Add);
        FirstPrimes = firstPrimes.Skip(1).ToArray();
        PrimeMultiples = FirstPrimes.Select(p => (p * p - 1) / 2).ToArray();
    }

    private void SieveSegment(BitArray segmentData, long segmentStart, long segmentEnd)
    {
        segmentData.SetAll(true);
        long segmentLength = segmentEnd - segmentStart;
        for (int i = 0; i < PrimeMultiples.Length; i++)
        {
            long current = PrimeMultiples[i] - segmentStart;
            long prime = FirstPrimes[i];
            if (current >= segmentLength) continue;

            while (current < segmentLength)
            {
                segmentData[(int)current] = false;
                current += prime;
            }
            PrimeMultiples[i] = segmentStart + current;
        }
    }

    public void ListPrimes(Action callback)
    {
        callback.Invoke(2);
        BitArray segmentData = new BitArray(BUFFER_LENGTH);
        long max = Length / 2;
        long segmentStart = 1;
        long segmentEnd = Math.Min(segmentStart + BUFFER_LENGTH, max);
        while (segmentStart < max)
        {
            SieveSegment(segmentData, segmentStart, segmentEnd);
            for (int i = 0; i < segmentEnd - segmentStart; i++)
            {
                if (segmentData[i])
                {
                    callback.Invoke((segmentStart + i) * 2 + 1);
                }
            }
            segmentStart = segmentEnd;
            segmentEnd = Math.Min(segmentStart + BUFFER_LENGTH, max);
        }
    }
}


Эта оптимизация ускоряет поиск простых чисел больще, чем вдвое:

Алгоритм: Сегментированная факторизация с базой 2
Ищет простые числа до миллиарда: 3.4 секунды.

10. Сегментированная (2,3)-факторизация


В предыдущей реализации мы исключили из рассмотрения числа, делящиеся на два; теперь исключим числа, делящиеся на 2 и 3. Тогда оставшиеся числа делятся на два класса: дающие 1 или 5 в остатке от деления на 6. Поэтому вместо одного битового массива используем два: PrimeMultiples_6kPlus1 и PrimeMultiples_6kPlus5; они будут хранить числа вида $(n - 1) / 6$ и $(n - 5) / 6$ для чисел, кратных $p$, имеющих 1 или 5 в остатке от деления на 6. Дополнительные усилия нужны, чтобы найти первое число вида $p^2+kp$, имеющее нужный остаток от деления на 6, но это решается простым перебором $k=0, 1, 2, 3, 4, 5$.

SegmentedWheel23.cs 
public class SegmentedWheel23 : ISieve
{
    const int BUFFER_LENGTH = 128 * 1024;

    private long Length;
    private long[] FirstPrimes;
    private long[] PrimeMultiples_6kPlus1;
    private long[] PrimeMultiples_6kPlus5;

    public SegmentedWheel23(long length)
    {
        Length = length;
        int firstChunkLength = (int)Math.Sqrt(length) + 1;
        SieveOfEratosthenes sieve = new SieveOfEratosthenes(firstChunkLength);
        List firstPrimes = new List();
        sieve.ListPrimes(firstPrimes.Add);
        FirstPrimes = firstPrimes.Skip(2).ToArray();
        PrimeMultiples_6kPlus1 = new long[FirstPrimes.Length];
        PrimeMultiples_6kPlus5 = new long[FirstPrimes.Length];
        for (int j = 0; j < FirstPrimes.Length; j++)
        {
            long prime = FirstPrimes[j];
            long val1 = prime * prime;
            while (val1 % 6 != 1) val1 += 2 * prime;
            PrimeMultiples_6kPlus1[j] = (val1 - 1) / 6;
            long val2 = prime * prime;
            while (val2 % 6 != 5) val2 += 2 * prime;
            PrimeMultiples_6kPlus5[j] = (val2 - 5) / 6;
        }
    }

    private void SieveSegment(BitArray segmentData_6kPlus1, BitArray segmentData_6kPlus5, long segmentStart, long segmentEnd)
    {
        segmentData_6kPlus1.SetAll(true);
        segmentData_6kPlus5.SetAll(true);
        long segmentLength = segmentEnd - segmentStart;
        for (int j = 0; j < FirstPrimes.Length; j++)
        {
            long prime = FirstPrimes[j];

            long val1 = PrimeMultiples_6kPlus1[j] - segmentStart;
            while (val1 < segmentLength)
            {
                segmentData_6kPlus1[(int)val1] = false;
                val1 += prime;
            }
            PrimeMultiples_6kPlus1[j] = val1 + segmentStart;

            long val2 = PrimeMultiples_6kPlus5[j] - segmentStart;
            while (val2 < segmentLength)
            {
                segmentData_6kPlus5[(int)val2] = false;
                val2 += prime;
            }
            PrimeMultiples_6kPlus5[j] = val2 + segmentStart;
        }
    }

    public void ListPrimes(Action callback)
    {
        callback.Invoke(2);
        callback.Invoke(3);
        callback.Invoke(5);
        BitArray segmentData_6kPlus1 = new BitArray(BUFFER_LENGTH);
        BitArray segmentData_6kPlus5 = new BitArray(BUFFER_LENGTH);
        long max = (Length + 5) / 6;
        long segmentStart = 1;
        long segmentEnd = Math.Min(segmentStart + BUFFER_LENGTH, max);
        while (segmentStart < max)
        {
            SieveSegment(segmentData_6kPlus1, segmentData_6kPlus5, segmentStart, segmentEnd);
            for (int i = 0; i < segmentEnd - segmentStart; i++)
            {
                if (segmentData_6kPlus1[i])
                {
                    long p = (segmentStart + i) * 6 + 1;
                    if (p >= Length) break;
                    callback.Invoke(p);
                }
                if (segmentData_6kPlus5[i])
                {
                    long p = (segmentStart + i) * 6 + 5;
                    if (p >= Length) break;
                    callback.Invoke(p);

                }
            }
            segmentStart = segmentEnd;
            segmentEnd = Math.Min(segmentStart + BUFFER_LENGTH, max);
        }
    }
}


Алгоритм: Сегментированная (2,3)-факторизация
Ищет простые числа до миллиарда: 2.6 секунды.

11. Сегментированная (2,3,5)-факторизация


Применим ту же идею для чисел, не делящихся на 2, 3 и 5. Тогда у нас будет не два битовых массива, а восемь, по одному на каждый остаток от деления на 30: $(1, 7, 11, 13, 17, 19, 23, 29)$. В остальном код не особо изменится, и даст нам еще небольшое ускорение.

SegmentedWheel235.cs 
public class SegmentedWheel235 : ISieve
{
    const int BUFFER_LENGTH = 128 * 1024;
    const int WHEEL = 30;
    const int WHEEL_PRIMES_COUNT = 3;
    private long[] WheelRemainders = { 1, 7, 11, 13, 17, 19, 23, 29 };
    private long[] SkipPrimes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29 };

    private long Length;
    private long[] FirstPrimes;
    private long[][] PrimeMultiples;

    public SegmentedWheel235(long length)
    {
        Length = length;
        int firstChunkLength = (int)Math.Sqrt(length) + 1;
        SieveOfEratosthenes sieve = new SieveOfEratosthenes(firstChunkLength);
        List firstPrimes = new List();
        sieve.ListPrimes(firstPrimes.Add);
        FirstPrimes = firstPrimes.Skip(WHEEL_PRIMES_COUNT).ToArray();
        PrimeMultiples = new long[WheelRemainders.Length][];
        for (int i = 0; i < WheelRemainders.Length; i++)
        {
            PrimeMultiples[i] = new long[FirstPrimes.Length];
            for (int j = 0; j < FirstPrimes.Length; j++)
            {
                long prime = FirstPrimes[j];
                long val = prime * prime;
                while (val % WHEEL != WheelRemainders[i]) val += 2 * prime;
                PrimeMultiples[i][j] = (val - WheelRemainders[i]) / WHEEL;
            }
        }
    }

    private void SieveSegment(BitArray[] segmentDatas, long segmentStart, long segmentEnd)
    {
        for (int i = 0; i < segmentDatas.Length; i++)
        {
            BitArray segmentData = segmentDatas[i];
            segmentData.SetAll(true);
            long segmentLength = segmentEnd - segmentStart;
            for (int j = 0; j < PrimeMultiples[i].Length; j++)
            {
                long current = PrimeMultiples[i][j] - segmentStart;
                long prime = FirstPrimes[j];
                if (current >= segmentLength) continue;

                while (current < segmentLength)
                {
                    segmentData[(int)current] = false;
                    current += prime;
                }
                PrimeMultiples[i][j] = segmentStart + current;
            }
        }
    }

    public void ListPrimes(Action callback)
    {
        foreach (long prime in SkipPrimes) if (prime < Length) callback.Invoke(prime);

        BitArray[] segmentDatas = new BitArray[WheelRemainders.Length];
        for (int i = 0; i < segmentDatas.Length; i++)
        {
            segmentDatas[i] = new BitArray(BUFFER_LENGTH);
        }
        long max = (Length + WHEEL - 1) / WHEEL;
        long segmentStart = 1;
        long segmentEnd = Math.Min(segmentStart + BUFFER_LENGTH, max);
        while (segmentStart < max)
        {
            SieveSegment(segmentDatas, segmentStart, segmentEnd);
            for (int i = 0; i < segmentEnd - segmentStart; i++)
            {
                long offset = (segmentStart + i) * WHEEL;
                for (int j = 0; j < segmentDatas.Length; j++)
                {
                    if (segmentDatas[j][i])
                    {
                        long p = offset + WheelRemainders[j];
                        if (p >= Length) break;
                        callback.Invoke(p);
                    }
                }
            }
            segmentStart = segmentEnd;
            segmentEnd = Math.Min(segmentStart + BUFFER_LENGTH, max);
        }
    }
}


Алгоритм: Сегментированная (2,3,5)-факторизация
Ищет простые числа до миллиарда: 2.3 секунды.

12. Оптимизированная сегментированная (2,3,5)-факторизация


Факторизация позволяет уменьшить количество чисел, которые мы вычеркиваем, но вычисление индекса байта и бита в битовом массиве (как при просеивании, так и при повторном проходе по массиву для перечисления простых чисел) все равно занимает большую часть времени. Поэтому воспользуемся уже знакомым нам трюком, пользуясь тем, что в факторизации с базой (2,3,5) восемь остатков: вместо восьми битовых массивов используем один массив байт, в котором каждый из восьми бит отвечает за свой остаток от деления на 30. Это позволит нам в массиве PrimeMultiples держать индекс байта, а индекс бита будет постоянным. Тогда внутренний цикл просеивания становится исключительно простым: 

while (current < segmentLength)
{
    segmentData[current] &= mask;
    current += prime;
}


При таком подходе повторный проход по массиву уже занимает больше времени, поскольку требует больше битовых операций. Однако любой байт может принимать всего 256 вариантов, поэтому мы можем заранее вычислить все 256 комбинаций остатков, и вообще избавиться от битовых операций. Например:
байт = 255 = 0b11111111, остатки: [1, 7, 11, 13, 17, 19, 23, 29]
байт = 9 = 0b00001001, остатки: [1, 13]
байт = 135 = 0b10000111, остатки: [1, 7, 11, 29]
Тогда получение простых чисел сводится к простому циклу: 

for (int i = 0; i < segmentData.Length; i++)
{
    byte data = segmentData[i];
    int[] offsets = OffsetsPerByte[data];
    for (int j = 0; j < offsets.Length; j++)
    {
        long p = (segmentStart + i) * 30 + offsets[j];
        if (p >= Length) break;
        callback.Invoke(p);
    }
}


Эта оптимизация ускоряет предыдущую реализацию больше, чем в два раза, позволяя найти простые числа до миллиарда меньше, чем за секунду!

OptimizedSegmentedWheel235.cs 
public class OptimizedSegmentedWheel235 : ISieve
{
    const int BUFFER_LENGTH = 200 * 1024;
    const int WHEEL = 30;
    const int WHEEL_PRIMES_COUNT = 3;

    private static long[] WheelRemainders = { 1, 7, 11, 13, 17, 19, 23, 29 };
    private static long[] SkipPrimes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29 };
    private static byte[] Masks = { 1, 2, 4, 8, 16, 32, 64, 128 };
    private static int[][] OffsetsPerByte;

    private long Length;
    private long[] FirstPrimes;
    private long[][] PrimeMultiples;

    static OptimizedSegmentedWheel235()
    {
        OffsetsPerByte = new int[256][];
        List offsets = new List();
        for (int b = 0; b < 256; b++)
        {
            offsets.Clear();
            for (int i = 0; i < WheelRemainders.Length; i++)
            {
                if ((b & Masks[i]) != 0) offsets.Add((int)WheelRemainders[i]);
            }
            OffsetsPerByte[b] = offsets.ToArray();
        }
    }

    public OptimizedSegmentedWheel235(long length)
    {
        Length = length;
        int firstChunkLength = (int)Math.Sqrt(length) + 1;
        SieveOfEratosthenes sieve = new SieveOfEratosthenes(firstChunkLength);
        List firstPrimes = new List();
        sieve.ListPrimes(firstPrimes.Add);
        FirstPrimes = firstPrimes.Skip(WHEEL_PRIMES_COUNT).ToArray();
        PrimeMultiples = new long[WheelRemainders.Length][];
        for (int i = 0; i < WheelRemainders.Length; i++)
        {
            PrimeMultiples[i] = new long[FirstPrimes.Length];
            for (int j = 0; j < FirstPrimes.Length; j++)
            {
                long prime = FirstPrimes[j];
                long val = prime * prime;
                while (val % WHEEL != WheelRemainders[i]) val += 2 * prime;
                PrimeMultiples[i][j] = (val - WheelRemainders[i]) / WHEEL;
            }
        }
    }

    private void SieveSegment(byte[] segmentData, long segmentStart, long segmentEnd)
    {
        for (int i = 0; i < segmentData.Length; i++) segmentData[i] = 255;
        long segmentLength = segmentEnd - segmentStart;

        for (int i = 0; i < WheelRemainders.Length; i++)
        {
            byte mask = (byte)~Masks[i];
            for (int j = 0; j < PrimeMultiples[i].Length; j++)
            {
                long current = PrimeMultiples[i][j] - segmentStart;
                if (current >= segmentLength) continue;
                long prime = FirstPrimes[j];

                while (current < segmentLength)
                {
                    segmentData[current] &= mask;
                    current += prime;
                }

                PrimeMultiples[i][j] = segmentStart + current;
            }
        }
    }

    public void ListPrimes(Action callback)
    {
        foreach (long prime in SkipPrimes) if (prime < Length) callback.Invoke(prime);

        long max = (Length + WHEEL - 1) / WHEEL;
        byte[] segmentData = new byte[BUFFER_LENGTH];
        long segmentStart = 1;
        long segmentEnd = Math.Min(segmentStart + BUFFER_LENGTH, max);
        while (segmentStart < max)
        {
            SieveSegment(segmentData, segmentStart, segmentEnd);
            for (int i = 0; i < segmentData.Length; i++)
            {
                long offset = (segmentStart + i) * WHEEL;
                byte data = segmentData[i];
                int[] offsets = OffsetsPerByte[data];
                for (int j = 0; j < offsets.Length; j++)
                {
                    long p = offset + offsets[j];
                    if (p >= Length) break;
                    callback.Invoke(p);
                }
            }
            segmentStart = segmentEnd;
            segmentEnd = Math.Min(segmentStart + BUFFER_LENGTH, max);
        }
    }
}


Алгоритм: Оптимизированная сегментированная (2,3,5)-факторизация
Ищет простые числа до миллиарда: 0.93 секунды.

Максимальное расстояние между простыми числами до триллиона


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

static void Main(string[] args)
{
    var sieve = new OptimizedSegmentedWheel235(1_000_000_000_000);
    long maxDistance = 0;
    long lastPrime = 0;
    Action action = (p) =>
    {
        long dist = p - lastPrime;
        if (dist > maxDistance)
        {
            Console.WriteLine($"{p} - {lastPrime} = {dist}");
            maxDistance = dist;
        }
        lastPrime = p;
    };
    sieve.ListPrimes(action);
}


Результат работы программы:
Простые числа до 1,000,000,000,000
Время работы: 26 минут 30 секунд
Максимальное расстояние между простыми числами: 540
между 738,832,927,927 и 738,832,928,467

Лучшие библиотеки


Рассказ был бы не полон без сравнения с лучшими из существующих библиотек по поиску простых чисел.
Самая известная библиотека — primegen, основанная на сильно оптимизированном решете Аткина. Она написана на C самим Даниэлем Бернштейном (соавтором Аткина) 20 лет назад и оптимизирована под Пентиум, поэтому для эффективного ее использования требуется поменять размер буфера в файле conf-words:
2048
This is B/32: the number of 32-bit words of space used in the primegen
inner loop. This should fit into the CPU's level-1 cache.

2048 works well on a Pentium-100.
3600 works well on a Pentium II-350.
4004 works well on an UltraSPARC-I/167.


На моем ноутбуке наилучшая производительность достигается при установке буфера 65536 байт. В библиотеку входит утилита primegaps, которая как раз и ищет максимальное расстояние между простыми числами. Время ее работы для простых чисел до триллиона — 19 минут.

Но на данный момент, пожалуй, самая быстрая библиотека — primesieve, написанная на c++ и включающая огромное количество оптимизаций, включая различные стратегии просеивания для маленьких, средних и больших простых чисел, автоматическое определение размера кешей процессора и многопоточность. На моем ноутбуке программа ищет простые числа до триллиона в 8 потоков за 70 секунд. Автор библиотеки Kim Walish подробно описывает оптимизации в коде, поэтому эта библиотека отлично подходит для изучения самых лучших на сегодняшний день оптимизаций и для практического применения.

© Habrahabr.ru