База данных простых чисел до ста миллиардов на коленке

image

Самый известный алгоритм для нахождения всех простых чисел, не больших заданного, – решето Эратосфена. Он замечательно работает для чисел до миллиардов, может быть, до десятков миллиардов, если аккуратно написан. Однако каждый, кто любит развлекаться с простыми числами, знает, что их всегда хочется иметь под рукой как можно больше. Как-то раз мне для решения одной задачи на хакерранке понадобилась in-memory база данных простых чисел до ста миллиардов. При максимальной оптимизации по памяти, если в решете Эратосфена представлять нечетные числа битовым массивом, его размер будет около 6 гигабайт, что в память моего ноутбука не влезало. Существует модификация алгоритма, гораздо менее требовательная по памяти (делящая исходный диапазон чисел на несколько кусков и обрабатывающая по одному куску за раз) – сегментированное решето Эратосфена, но она сложнее в реализации, и результат целиком в память все равно не влезет. Ниже предлагаю вашему вниманию алгоритм почти такой же простой, как и решето Эратосфена, но дающий двукратную оптимизацию по памяти (то есть, база данных простых чисел до ста миллиардов будет занимать около 3 гигабайт, что уже должно влезать в память стандартного ноутбука).

Для начала вспомним, как работает решето Эратосфена: в массиве чисел от 1 до N вычеркиваются числа, делящиеся на два, потом числа, делящиеся на 3, и так далее. Оставшиеся после вычеркивания числа и будут простыми:

2 3 . 5 . 7 . 9 . 11 . 13 . 15 . 17 . 19 . 21 . 23 . 25 . 27  . 29  ...
2 3 . 5 . 7 . . . 11 . 13 .  . . 17 . 19 .  . . 23 . 25 .  .  . 29  ...
2 3 . 5 . 7 . . . 11 . 13 .  . . 17 . 19 .  . . 23 .  . .  .  . 29  ...

Код на c#
class EratospheneSieve
{
    private bool[] Data;

    public EratospheneSieve(int length)
    {
        Data = new bool[length];

        for (int i = 2; i < Data.Length; i++) Data[i] = true;

        for (int i = 2; i * i < Length; i++)
        {
            if (!Data[i]) continue;
            for (int d = i * i; d < Length; d += i) Data[d] = false;
        }
    }

    public int Length => Data.Length;

    public bool IsPrime(int n) => Data[n];
}


На моем ноутбуке до миллирарда считает 15 секунд и потребляет гигабайт памяти.

Алгоритм, который мы будем использовать, называется Wheel Factorization. Чтобы понять его суть, напишем натуральные числа табличкой по 30 штук в ряд:

 0  1  2  3  4  5 ... 26 27 28 29
30 31 32 33 34 35 ... 56 57 58 59
60 61 62 63 64 65 ... 86 87 88 89
...

И вычеркнем все числа, делящиеся на 2, 3, и 5:

.  1 . . . . .  7 . . . 11 . 13 . . . 17 . 19 . . . 23 . . . . . 29
. 31 . . . . . 37 . . . 41 . 43 . . . 47 . 49 . . . 53 . . . . . 59
. 61 . . . . . 67 . . . 71 . 73 . . . 77 . 79 . . . 83 . . . . . 89
...

Видно, что в каждом ряду вычеркиваются числа на одних и тех же позициях. Так как мы вычеркивали только составные числа (за исключением, собственно, чисел 2, 3, 5), то простые числа могут быть только на невычеркнутых позициях.

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

image

Техническая реализация


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

class Wheel235
{
    private byte[] Data;

    public Wheel235(long length)
    {
        ...
    }
    public bool IsPrime(long n)
    {
        ...
    }
    ...
}

Для простоты реализации, будем округлять length вверх до ближайшего числа, делящегося на 30. Тогда

class Wheel235
{
    private byte[] Data;

    public Wheel235(long length)
    {
        // ensure length divides by 30
        length = (length + 29) / 30 * 30; 
        Data = new byte[length/30];
        ...
    }

    public long Length => Data.Length * 30L;

    ...
}

Далее, нам нужны два массива: один переводит индекс 0..29 в номер бита, второй, наоборот, номер бита в позицию в тридцатке. Используя эти массивы и немного битовой арифметики, мы строим однозначное соответствие между натуральными числами и битами в массиве Data. Для простоты мы используем только два метода: IsPrime(n), позволяющий узнать, является ли число простым, и ClearPrime(n), устанавливающий соответствующий бит в 0, маркируя число n как составное:

class Wheel235
{
    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 bool IsPrime(long n)
    {
        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);
    }
}

Заметим, что, если число делится на 2, 3 или 5, то номер бита будет -1, что означает, что число заведомо составное. Ну и, конечно, мы обрабатываем специальный случай n <= 5.

Теперь, собственно, сам алгоритм, который практически дословно повторяет решето Эратосфена:

    public Wheel235(long length)
    {
        length = (length + 29) / 30 * 30;
        Data = new byte[length / 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);
        }
    }

Соответственно, целиком получившийся класс:

выглядит так
class Wheel235
{
    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 Wheel235(long length)
    {
        // ensure length divides by 30
        length = (length + 29) / 30 * 30;
        Data = new byte[length / 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 long Length => Data.Length * 30L;

    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);
    }
}


Есть только одна маленькая проблемка: этот код не работает для чисел, больших примерно 65 миллиардов! При попытке её запустить с такими числами программа падает с ошибкой:

Unhandled Exception: System.OverflowException: Array dimensions exceeded supported range.

Проблема в том, что в c# массивы не могут иметь больше 2^31 элементов (видимо, потому что внутренняя имплементация использует четырехбайтный индекс массива). Один из вариантов – вместо массива byte[] создавать, например, массив long[], но это немного усложнит битовую арифметику. Для простоты мы пойдем другим путем, создав специальный класс, симулирующий нужный массив, держащий внутри два коротких массива. Заодно мы ему дадим возможность сохранять себя на диск, чтобы можно было один раз вычислить простые числа, и потом пользоваться базой повторно:

    class LongArray
    {
        const long MAX_CHUNK_LENGTH = 2_000_000_000L;
        private byte[] firstBuffer;
        private byte[] secondBuffer;

        public LongArray(long length)
        {
            firstBuffer = new byte[Math.Min(MAX_CHUNK_LENGTH, length)];
            if (length > MAX_CHUNK_LENGTH) secondBuffer = new byte[length - MAX_CHUNK_LENGTH];
        }

        public LongArray(string file)
        {
            ...
        }

        public long Length => firstBuffer.LongLength + (secondBuffer == null ? 0 : secondBuffer.LongLength);

        public byte this[long index]
        {
            get => index < MAX_CHUNK_LENGTH ? firstBuffer[index] : secondBuffer[index - MAX_CHUNK_LENGTH];
            set
            {
                if (index < MAX_CHUNK_LENGTH) firstBuffer[index] = value; 
                else secondBuffer[index - MAX_CHUNK_LENGTH] = value;
            }
        }

        public void Save(string file)
        {
            ...
        }

    }


Наконец, объединим все шаги в

окончательный вариант алгоритма
class Wheel235
{
    class LongArray
    {
        const long MAX_CHUNK_LENGTH = 2_000_000_000L;
        private byte[] firstBuffer;
        private byte[] secondBuffer;

        public LongArray(long length)
        {
            firstBuffer = new byte[Math.Min(MAX_CHUNK_LENGTH, length)];
            if (length > MAX_CHUNK_LENGTH) secondBuffer = new byte[length - MAX_CHUNK_LENGTH];
        }

        public LongArray(string file)
        {
            using(System.IO.FileStream stream = System.IO.File.OpenRead(file))
            {
                long length = stream.Length;
                firstBuffer = new byte[Math.Min(MAX_CHUNK_LENGTH, length)];
                if (length > MAX_CHUNK_LENGTH) secondBuffer = new byte[length - MAX_CHUNK_LENGTH];
                stream.Read(firstBuffer, 0, firstBuffer.Length);
                if(secondBuffer != null) stream.Read(secondBuffer, 0, secondBuffer.Length);
            }
        }

        public long Length => firstBuffer.LongLength + (secondBuffer == null ? 0 : secondBuffer.LongLength);

        public byte this[long index]
        {
            get => index < MAX_CHUNK_LENGTH ? firstBuffer[index] : secondBuffer[index - MAX_CHUNK_LENGTH];
            set
            {
                if (index < MAX_CHUNK_LENGTH) firstBuffer[index] = value; 
                else secondBuffer[index - MAX_CHUNK_LENGTH] = value;
            }
        }

        public void Save(string file)
        {
            using(System.IO.FileStream stream = System.IO.File.OpenWrite(file))
            {
                stream.Write(firstBuffer, 0, firstBuffer.Length);
                if (secondBuffer != null) stream.Write(secondBuffer, 0, secondBuffer.Length);
            }
        }

    }

    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 LongArray Data;

    public Wheel235(long length)
    {
        // ensure length divides by 30
        length = (length + 29) / 30 * 30;
        Data = new LongArray(length / 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 Wheel235(string file)
    {
        Data = new LongArray(file);
    }

    public void Save(string file) => Data.Save(file);

    public long Length => Data.Length * 30L;

    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);
    }

}


В результате мы получили класс, который может вычислить простые числа, не превосходящие ста миллиардов, и сохранить результат на диск (у меня этот код на ноутбуке работал 50 минут; размер созданного файла получился примерно 3 гигабайта, как и ожидалось):

        long N = 100_000_000_000L;
        Wheel235 wheel = new Wheel235(N);
        wheel.Save("./primes.dat");

А потом зачитать с диска и использовать:

        DateTime start = DateTime.UtcNow;
        Wheel235 loaded = new Wheel235("./primes.dat");
        DateTime end = DateTime.UtcNow;
        Console.WriteLine($"Database of prime numbers up to {loaded.Length:N0} loaded from file to memory  in {end - start}");
        long number = 98_000_000_093L;
        Console.WriteLine($"{number:N0} is {(loaded.IsPrime(number) ? "prime" : "not prime")}!");
AdBlock похитил этот баннер, но баннеры не зубы — отрастут

Подробнее
Реклама

Комментарии 21

    0

    А почему бы не использовать диск как swap?

      +1
      При генерации простых чисел будет много random writes, боюсь, будет слишком медленно. Лучше тогда использовать сегментированное решето Эратосфена, вычислять в памяти по кускам и скидывать найденные куски на диск. Тогда можно будет генерить простые числа, в принципе, до бесконечности.
      А при использовании готовой базы — да, конечно, вместо того, чтобы загружать базу в память, можно просто читать с диска по нужному смещению.
      +2
      красиво
        +6
        Можно сжать, если проверять на простоту с помощью МТФ — хранить только те, которые не являются простыми но проходят проверку. Или Тест Агравала — Каяла — Саксены
          0
          Интересная идея! Насколько я понимаю, все же, скорость проверки не удастся сохранить O(1).
          Впрочем, далеко не во всех задачах требуется константная скорость проверки, тогда, действительно, можно пожертвовать скоростью ради уменьшения размера базы. В предельном случае, если я правильно понимаю, можно вообще ничего не хранить, и применять AKS-тест каждый раз заново.
            0

            Если Вы ограничиваете себе верхнюю планку чисел, то сложность в любом алгоритме будет константной.

              0
              Верхняя планка чисел всегда ограничена количеством памяти на компьютере и разрядностью. Поправьте меня, если я ошибаюсь, но стандартная de-facto вычислительная модель, если не оговорено обратное — машина Тьюринга, что эквивалентно компьютеру с бесконечным количеством памяти.
                +1

                Для оценки сложности логичнее использовать RAM-машину с бесконечной памятью, а не машину Тюринга, потому что даже простое обращение к ячейке N в RAM-машине имеет сложность O(1), а МТ — O(N).


                В любом случае ссылаться на конечность данных реальных компьютеров при анализе алгоритмов — верный признак недопонимания. С другой стороны, и константу перед главным членом в реальности тоже надо учитывать.

                  –2

                  Ой, а в чем именно у меня верный признак недопонимания? Где именно ошибка в моем комментарии, что


                  Если Вы ограничиваете себе верхнюю планку чисел, то сложность в любом алгоритме будет константной.

                  Укажите, будьте столь любезны, понимающий Вы наш.

                    0
                    логичнее использовать RAM-машину

                    Понятно, спасибо. Я почему-то думал, что MT — это обобщенное название класса машин, эквивалентных классической MT с одномерной лентой.
              0
              Можно сжать, если проверять на простоту с помощью МТФ — хранить только те, которые не являются простыми но проходят проверку.

              Нашел в Викпедии, что Baillie–PSW primality test, который является комбинацией двух вероятностных тестов на простоту (Ферма и Лукаса), вообще не имеет исключений для 64-битных чисел. Больше того, даже для произвольно больших чисел за сорок лет ни одного исключения так и не нашли. Так что Ваша идея тут кстати: размер базы будет нулевой практически всегда. А если вдруг кто-нибудь найдет составное число, проходящее этот тест, он даже получит премию в 30 долларов!
              А из более поздних исследований нашел эту статью 2015 года — они утверждают, что их алгоритм самый быстрый (скорость проверки — чуть больше миллиона 64-битных чисел в секунду) из класса алгоритмов, не требующих предварительных вычислений.
              Забавно, но я у них нашел такую цитату:
              If we expect that we’ll need to test a significantly large number of small integers for primality, the best solution might be precomputing all possible answers
              and then answering each query in constant time. The canonical implementation
              uses one bit for each odd number, i.e., 2b/16 bytes of memory if the valid inputs are b-bit integers

              То есть, они тоже считают, что для некоторых случаев битовая таблица — лучшее решение, но даже они говорят о бите на каждое нечетное число! Хотя Wheel factorization широко известен, почему-то даже ученые его игнорируют при практическом применении
                0

                Я пытался написать максимально быстрое решето, и упаковка 30 чисел в байт проигрывает по скорости хранению только нечетных из-за необходимости вычисления остатка.
                Сработал вариант хранения восьми отдельных массивов для отдельных остатков.

                  +1
                  Можете подробнее про Ваш вариант рассказать?
                  Оптимизация таких числодробилок иногда бывает очень нетривиальна. У меня, например, скорость увеличивается в полтора раза, только если закомментарить строку if (n >= Length) throw new ArgumentException(«Number too big»);
                  Возможно, она как-то препятствует jit-у инлайнить метод, не изучал вопрос.

                  Кстати, получить итератор по всем простым числам в возрастающем порядке (самый частый use-case) можно эффективнее, чем проверять на простоту числа последовательно от 1 до ста миллионов (используя, как раз, массив BIT_TO_INDEX), просто я старался сделать код как можно проще, даже ценой ухудшения производительности.
                    +1
                    Немного кода
                    /// <summary>
                    /// Test only numbers like Size*k + Offsets[i]
                    /// </summary>
                    private struct Wheel {
                        public int       Size;
                        public int[]     Offsets;
                        public int[]     Deltas;
                        public DivMod[,] DivMods;
                    }
                    private static T EratospheneWheel<T>(T consumer, int bound, Wheel wheel)
                        where T : struct, IConsumer {
                        int low    = (int) Math.Sqrt(bound);
                        int length = bound / wheel.Size + 1;
                        var sieves = new BitArray[wheel.Offsets.Length];
                        for (int i = 0; i < sieves.Length; i++) {
                            sieves[i] = new BitArray(length, true);
                        }
                        int k = 0, num = wheel.Offsets[1], sieve = 1;
                        while (num <= low) {
                            if (sieves[sieve][k]) {
                                consumer.Consume(num);
                                for (int other = 0; other < sieves.Length; other++) {
                                    var divMod = wheel.DivMods[sieve, other];
                                    var bits   = sieves[divMod.Idx];
                                    var start = k * (num + wheel.Offsets[other]) + divMod.Div;
                                    for (int idx = start; idx < length; idx += num) {
                                        bits[idx] = false;
                                    }
                                }
                            }
                            num += wheel.Deltas[sieve];
                            sieve++;
                            if (sieve == sieves.Length) {
                                ++k;
                                sieve = 0;
                            }
                        }
                        while (num <= bound) {
                            if (sieves[sieve][k]) {
                                consumer.Consume(num);
                            }
                            num += wheel.Deltas[sieve];
                            sieve++;
                            if (sieve == sieves.Length) {
                                ++k;
                                sieve = 0;
                            }
                        }
                        return consumer;
                    }

                    Это самый быстрый вариант, который у меня получился. Наилучший размер колеса — 210.
                    Этот код суммирует простые числа до миллиарда за 5 секунд.
                    И конечно в числодробилке не надо использовать ничего сверх минимума. Только массивы, только хардкор. Try и Exception сильно бьют по призводительности.
                    Деления и остатки вычисляются заранее, но возможно компилятор сможет оптимизировать деление на константу.

                    0
                    Когда я писал на Java, 64-битовый массив немного проигрывал по скорости массиву boolean, хотя и экономил 8х память, из-за необходимости считать арифметику.
                      0
                      А исполнялась программа под 64-битную ява машину или 32-битную? В 32-битных процессах длинная арифметика медленнее, так как 64-битные регистры процессора недоступны.
                      Хотя я согласен, битовая арифметика будет замедляет алгоритм в любом случае.
                      Теоретически, если искать простые числа не в одном большом массиве на гигабайты, а разделить его на куски по несколько мегабайт, чтобы весь битовый массив поместился в кеш третьего уровня, может, тогда ускорение доступа к памяти может исправить замедление от битовых операций.
                0

                А в чём трудности с битовой арифметикой при использовании ulong вместо byte? Насколько я понимаю, она остаётся той же самой, просто массивы INDEX_TO_BIT и BIT_TO_INDEX увеличиваются в 8 раз (с соответствующим заполнением значений, шаблон такого заполнения легко вычисляется), а деление на 30 превращается в деление на 240.

                  +1
                  Вы абсолютно правы, у меня в изначальной версии был именно ulong[]. Просто я не упомянул еще некоторые неудобства:
                  1) массив байт можно сохранить в файл одним вызовом File.Write
                  2) в .NET Framerork массив ulong-ов все равно просто так не создать
                  3) если мы захотим считать простые числа, например, до триллиона (вдруг у нас есть 30 гб памяти), массив ulong сам станет больше двух миллиардов элементов.
                  0
                  А не проверяли сжимаемость файла архиваторами?
                    0
                    Сейчас проверил (deflate, normal compression) сжался с 3.2 гб до 2.1, я, честно говоря, не ожидал
                    0
                    Вернусь немного, в некотором смысле, в прошлое:
                    • Начинаем с 0.
                      На каждое натуральное число нужен 1 бит
                    • Далее 2: можно не хранить чётные, экономим 50% («правильный» остаток — 1, «неправильный» — 0), сама двойка в уме.
                    • Далее к двойке присоединяется тройка, остатки от деления на 6, соответственно, 1 и 5; на каждые 6 чисел (натуральных) требуется 2 бита, ещё полуторакратное улучшение, 2 и 3 в уме.
                    • Тут к двойке с тройкой присоединяется пятёрка: делитель 30, остатков 8, соответственно, 30 чисел в восьми битах (по сравнению с предыдущим вариантом выигрыш всего 20%: до этого те же 30 чисел занимали 10 бит), 2, 3 и 5 в уме.
                    • Тут на сцену могла бы выйти семёрка, но 210 чисел в 48-ми битах — это дополнительный выигрыш всего около 14% (уже 35 чисел на те же 8 бит), а код ещё более сложный и он должен работать, а значит тратить процессорное и наше время.

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

                    Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                    Самое читаемое