Быстрый поиск простых чисел

[B@R5uk]

Задался я тут задачей по-быстрому найти все простые числа, которые вмещаются в тип данных int. Оглядываясь назад, глупая, наверное, затея, учитывая что этих чисел на 800 МБ и поиск совсем не быстрый. Но тем не менее.

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

Единственный вопрос в том, как определять окончание проверки кандидата: 1) можно на каждом шаге цикла проверки находить квадрат делителя и сравнивать его с кандидатом: if (primeNominee < primeTester * primeTester)...; 2) а можно заранее посчитать корень primeLimit = (long) Math.sqrt(primeNominee); и на каждом шаге делители сравнивать с найденным порогом: if (primeLimit < primeTester)... На практике второй вариант оказался чуть-чуть быстрее. В результате этих соображений получается такая программа:

Используется синтаксис Java

class Primes {
        
        private static final int primesNumber = 1000000;
        private static final long[] primesList = new long[primesNumber];
        static {
                int primeIndex, testerIndex;
                long primeNominee, primeTester, primeLimit;
                primesList[0] = 2;
                primesList[1] = 3;
                primeIndex = 2;
                primeNominee = 5;
                outerLoop: while (true) {
                        testerIndex = 1;
                        primeLimit = (long) Math.sqrt(primeNominee);
                        while (true) {
                                primeTester = primesList[testerIndex];
                                if (primeLimit < primeTester) {
                                        primesList[primeIndex] = primeNominee;
                                        ++primeIndex;
                                        if (primesNumber == primeIndex) {
                                                break outerLoop;
                                        }
                                        break;
                                }
                                if (0 == primeNominee % primeTester) {
                                        break;
                                }
                                ++testerIndex;
                        }
                        primeNominee += 2;
                }
        }
        
        public static long getPrime(int index) {
                return primesList[index];
        }
        
        public static long getLastPrime() {
                return primesList[primesNumber - 1];
        }
}

public class Primes_Search {
        public static void main(String[] args) {
                long startTime, stopTime;
                System.out.println(String.format("%13d", Integer.MAX_VALUE));
                startTime = System.nanoTime();
                System.out.println(String.format("%13d", Primes.getLastPrime()));
                stopTime = System.nanoTime();
                System.out.println(String.format("\n%8.4f", 1e-9 * (stopTime - startTime)));
        }
}
 

Однако, я был разочарован, так как работает этот алгоритм чрезвычайно медленно (не то, чтобы у меня было с чем сравнивать). Единственные решения для улучшения алгоритма, которые я смог придумать не меняя его сути, заключаются в следующем. Можно перебирать не просто нечётные числа с шагом 2 (то есть не делящиеся нацело на 2), а нечётные с переменным шагом 2, 4, которые не будут делиться не только на 2, но и на 3. Этот подход можно улучшить, если так же исключить числа, делящиеся на 5, взяв переменный шаг 2, 4, 2, 4, 6, 2, 6, 4. Связанный с этим подходом в английском есть термин Wheel factorization.

Кроме того, от медленного вычисления корня можно избавиться, если заметить, что он плавно растёт вместе сплавным ростом кандидата в простое число, а следующий простой делитель в список делителей для проверки добавляется, когда кандидат пересекает очередной квадрат простого числа. Все эти соображения дают такой улучшенный алгоритм:

Используется синтаксис Java

class Primes {
        
        private static final int primesNumber = 1000000;
        private static final long[] primesList = new long[primesNumber];
        static {
                int primeIndex, testerIndex, stepIndex, limitIndex;
                long primeNominee, limitSquare;
                primesList[0] = 2;
                primesList[1] = 3;
                primesList[2] = 5;
                primesList[3] = 7;
                primeIndex = 4;
                primeNominee = 11;
                int[] stepsList = {2, 4, 2, 4, 6, 2, 6, 4};
                stepIndex = 0;
                limitIndex = 3;
                limitSquare = 49;
                outerLoop: while (true) {
                        if (limitSquare == primeNominee) {
                                ++limitIndex;
                                limitSquare = primesList[limitIndex];
                                limitSquare *= limitSquare;
                        } else {
                                testerIndex = 3;
                                while (true) {
                                        if (limitIndex == testerIndex) {
                                                primesList[primeIndex] = primeNominee;
                                                ++primeIndex;
                                                if (primesNumber == primeIndex) {
                                                        break outerLoop;
                                                }
                                                break;
                                        }
                                        if (0 == primeNominee % primesList[testerIndex]) {
                                                break;
                                        }
                                        ++testerIndex;
                                }
                        }
                        primeNominee += stepsList[stepIndex];
                        ++stepIndex;
                        if (8 == stepIndex) {
                                stepIndex = 0;
                        }
                }
        }
        
        public static long getPrime(int index) {
                return primesList[index];
        }
        
        public static long getLastPrime() {
                return primesList[primesNumber - 1];
        }
}
 

Однако, последние два улучшения являются по сути украшательством, так как прибавка в скорости составляет чуть больше 4%. Очевидно, что основное время выполнения программы набирается на операции деления с остатком ближе к концу работы при проверке больших простых чисел или же чисел с большими простыми множителями. Откидывание чисел, кратных 2, 3 и 5, практически ничего не даёт, так как такие числа в первом примере отсеиваются сразу на первых итерациях внутреннего цикла.

В связи с этим я начинаю задумываться о принципиально других алгоритмах поиска простых чисел. Первый в очереди — решето Эратосфена. Однако, прежде чем переходить к новому алгоритму со всеми его особенностями, хотелось бы оценить сложность наивного алгоритма. С этим у меня проблемы. Внешний цикл имеет количество итераций $\pi \left( N \right)\sim \frac{N}{\ln N}$, где N — это граница сверху на искомые простые числа, а не количество к поиску, как в моих примерах. С внутренним же всё сложнее: количество его итераций скачет туда-сюда и в худшем случае составляет где-то $\pi \left( \sqrt{N} \right)\sim \frac{2\sqrt{N}}{\ln N}$. Хотелось бы теперь просто перемножить эти две величины и получить что-то типа $O\left( \frac{{{N}^{{3}/{2}}}}{{{\ln }^{2}}N} \right)$, но на самом деле вторую величину надо интегрировать (читай: суммировать по всем нечётным n от 3 до N), учитывая при этом её случайным образом пляшущий характер. Как быть?

П.С. По решету Эратосфена рекомендации тоже с радостью принимаются.

[Dmitriy40]

B@R5uk
Решетом Эратосфена все числа до миллиарда находятся за доли секунды. До $2^{32}$ видимо за секунду-две. Зачем ещё убыстрять?
Ощущение что Вы в который раз изобретаете давно известные вещи.
Плюс в сети видел выложены все простые до $10^{12}$, можно просто взять и скачать, это конечно на порядки дольше вычислений, но зато не надо писать программ. Например в A000040 приведена ссылка "Matthew Parker, The first billion primes (7-Zip compressed file) [a large file]" на файл a000040_1B.
Ну и есть Primesieve, которая за те же несколько секунд выдаст вам в stdout все простые в нужном вам диапазоне.

[wrest]

B@R5uk
По решету, гляньте вот сюдой: https://rosettacode.org/wiki/Sieve_of_Eratosthenes
Там на 150+ языках реализация его. Может какие идеи почерпнёте.

[B@R5uk]

Dmitriy40 в сообщении #1435511 писал(а):

Решетом Эратосфена все числа до миллиарда находятся за доли секунды. До $2^{32}$ видимо за секунду-две.

Это если на компьютере памяти много.

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

wrest в сообщении #1435516 писал(а):

По решету, гляньте вот сюдой: https://rosettacode.org/wiki/Sieve_of_Eratosthenes

Для джавы там решения красивые, в том плане, что они компактные или на всю катушку используют языковые возможности типа структур данных. Но это далеко не оптимальные решения для поиска больших массивов простых чисел.

Вообще, главная проблема с решетом Эратосфена (а я ещё не начал его реализовывать), как мне кажется, в расходовании памяти. Сразу весь желаемый диапазон обработать не получится, поэтому придётся разбивать его на блоки. В этом могут возникнуть тонкости и проблемы.

[Dmitriy40]

B@R5uk в сообщении #1435529 писал(а):

Dmitriy40 в сообщении #1435511 писал(а):

Решетом Эратосфена все числа до миллиарда находятся за доли секунды. До $2^{32}$ видимо за секунду-две.

Это если на компьютере памяти много.

Давайте прикинем: до $2^{32}$ будет 203млн простых, это 813МБ (их в памяти хранить не нужно, можно сразу в файл скидывать), для их получения достаточно работать с простыми до $2^{16}$, коих всего $6542$ и которые можно хранить в 16 битных ячейках (раз вам памяти жалко). Для ускорения понадобится второй массив на $6542$ элемента, тоже по 2 байта. Итого, если битовый массив решета удерживать в 32КБ кэша L1 (что даст полумиллионный проверяемый интервал на проход), то $32768 + 6542 \cdot 4 = 58936$, меньше 60КБ. Уж столько памяти есть на любом компьютере.

B@R5uk в сообщении #1435529 писал(а):

В этом могут возникнуть тонкости и проблемы.

И возникнут. Смотрите segmentation алгоритмы, они есть и по великолепной ссылке wrest выше, и вообще.

Вообще, посоветую Вам получше определиться что же хотите: 1) получить максимально быстрый алгоритм получения всех простых до $2^{32}$, 2) получить максимально быстрый алгоритм получения всех простых в миллиардном/миллионном интервале начиная с произвольного места, 3) ещё что-то. Потому что это сильно разные задачи. Например из собственного опыта, скорость генерации простых около $10^{12}$, около $10^{15}$ и около $10^{18}$ различается по сильно разным причинам (а около $10^9$ просто не интересно так как слишком быстро). Т.е. для получения максимальной скорости оптимизировать надо по разному. А для чисел больше $10^{20..21}$ придётся думать что с памятью делать, список необходимых простых уже в неё не помещается. И во всех случаях придётся учитывать особенности кэшей (и не только их объёмы).
Ну и использовать Java для высокоскоростных приложений как минимум странно.

[B@R5uk]

Dmitriy40 в сообщении #1435535 писал(а):

Для ускорения понадобится второй массив на $6542$ элемента, тоже по 2 байта.

Я правильно понимаю, что он нужен, чтобы хранить индексы вычёркиваемых составных чисел для следующего блока решета?

[arseniiv]

(Оффтоп)

Dmitriy40 в сообщении #1435535 писал(а):

Ну и использовать Java для высокоскоростных приложений как минимум странно.

Не сильно странно, она джитится.

[Dmitriy40]

B@R5uk в сообщении #1435544 писал(а):

Dmitriy40 в сообщении #1435535 писал(а):

Для ускорения понадобится второй массив на $6542$ элемента, тоже по 2 байта.

Я правильно понимаю, что он нужен, чтобы хранить индексы вычёркиваемых составных чисел для следующего блока решета?

Ну можно и так сказать. Хотя реально там сидят смещения в массиве решета для каждого простого, фактически $x_i=a \bmod p_i$ ($a$ - величина начала сегмента, $p$ - простое, $x$ - смещение в сегменте следующего кратного этому простому) или как-то так. Чтобы не высчитывать их снова для каждого сегмента (а это до $6541$ совсем не быстрых деления!).

arseniiv
JIT компиляция это великолепно, но я вообще не верю в скорость универсальных инструментов, а типы данных в Java слишком универсальны. Конечно там может быть на них навёрнута оптимизация, но она не может быть статической (на этапе компиляции), т.е. скорость должна быть таки меньше чем на С/С++ (про ассемблер и SSE/AVX ладно уж не буду). Но разумеется лучше бы проверить прямым тестом.

[B@R5uk]

Dmitriy40 в сообщении #1435554 писал(а):

Чтобы не высчитывать их снова для каждого сегмента (а это до $6541$ совсем не быстрых деления!).

Я сначала сделал с вычислением остатков, а потом понял, что эту часть программы можно оптимизировать, если ввести дополнительный массив. С одной стороны я был потрясён этим, так что задумался о том, переделать ли или оставить как проходной вариант. А с другой — был огорчён, что вы, судя по всему, это уже знали, но не упомянули, посоветовав копаться в чужом коде. Как видите, в изобретении велосипедов куча фана. Проходной код решил пока оставить неоптимизированным, чтобы было. (В том числе для сравнения быстродействия).

Dmitriy40 в сообщении #1435554 писал(а):

фактически $x_i=a \bmod p_i$ ($a$ - величина начала сегмента, $p$ - простое, $x$ - смещение в сегменте следующего кратного этому простому) или как-то так.

На самом деле там более сложная формула. Я с ней долго и упорно боролся, прежде чем мне удалось её одолеть. Как-то так должно выглядеть:

$${{x}_{i}}=\left\{ \begin{array}{*{35}{l}} 0, & 0=a\bmod {{p}_{i}} \\ {{p}_{i}}-a\bmod {{p}_{i}}, & 0\ne a\bmod {{p}_{i}} \\ \end{array} \right.$$

Величина $x_i$ должна быть такой, чтобы при прибавлении к a получалось число, кратное $p_i$. В моём же случае, программа работает только с нечётными числами, поэтому эта формула ещё сложнее:

$${{x}_{i}}=\left\{ \begin{array}{*{35}{l}} 0, & 0=\left( a-{{p}_{i}} \right)\bmod \left( 2{{p}_{i}} \right) \\ 2{{p}_{i}}-\left( a-{{p}_{i}} \right)\bmod \left( 2{{p}_{i}} \right), & 0\ne \left( a-{{p}_{i}} \right)\bmod \left( 2{{p}_{i}} \right) \\ \end{array} \right.$$

При добавлении $x_i$ к a должно получаться не просто число, кратное $p_i$, оно должно быть ещё и нечётным.

-- 17.01.2020, 06:25 --

Плюс тут я нашёл ещё одну тонкость. Оригинальное решето Эратосфена подразумевает, что простые числа находятся и сразу же используются в процессе вычёркивания. В случае же сегментированного решета так будет происходить только с самым первым сегментом. Со всеми остальными сегментами простые числа будут браться из ранее подготовленного массива, а найденные после вычёркивания, если и будут использоваться, то только в следующих сегментах. Чтобы не заморачиваться с процедурой, которая будет одинаково хорошо работать в обоих случаях, я исключил первый блок, заменив его наивным поиском простых чисел вплоть до нужной величины (зависит от размера сегмента).

[Dmitriy40]

B@R5uk в сообщении #1435580 писал(а):

На самом деле там более сложная формула.

Я знаю. Но Вам что, сразу готовый код дать что ли? Я вообще-то этот код на асме переписывал (чтобы получить выигрыш от коротких делений когда это возможно). Впрочем вот тут его написал и на С.

B@R5uk в сообщении #1435580 писал(а):

А с другой — был огорчён, что вы это, судя по всему, знали, но не упомянули, посоветовав копаться в чужом коде.

Кажется я советовал копаться в алгоритмах, не коде, ну да ладно. А так да, знал и не упомянул, как наверняка и кучу других вещей, видимо посчитал замену повторных вычислений на табличку (массив) достаточно очевидной, ведь это же известный способ оптимизации по скорости.

Что ещё не упомянул: сегменты желательно делать размером с кэш (L1 или лучше L2) или хотя бы не больше. Новые $x_i$ в массив не сохранять до прохода всего сегмента (записи в память относительно дорогие), пользоваться временной переменной. Лет 6 назад я даже где-то здесь целый список улучшений решета приводил, пунктов 13 кажется было, из них сам реализовал лишь 7-9, не могу найти то сообщение (где-то рядом с этим должно быть). Нечётные и сегменты там были в первой половине списка.

[B@R5uk]

Мда! То, что раньше работало почти 5 минут, теперь делается чуть больше, чем за секунду. На всякий случай: $p_{10'000'000}=179'424'673$ ?

-- 17.01.2020, 06:54 --

Dmitriy40 в сообщении #1435584 писал(а):

Но Вам что, сразу готовый код дать что ли?

Ни в коем случае! Лучше идеи, а реализацию я сам сделаю.

[Dmitriy40]

О, нашёл тот список! Только там везде времена уже устарели, я недавно ускорил свою программу почти на порядок, с 12-16 секунд до 2 секунд.
Простые по номеру и наоборот можно проверить тут: "The 10,000,000th prime is 179,424,673."
Насчёт идей почитайте описание алгоритмов у primesieve, она одна из самых быстрых в мире.

[B@R5uk]

Dmitriy40, спасибо за ссылки. Теперь я знаю, откуда вы знаете сколько простых чисел до заданного числа! В приведённом вами списке оптимизаций есть вещи, о которых я задумывался ещё даже не начав реализовывать решето. Например, пункт 12. Пока, я правда, не придумал, как поступать с большими делителями, которые попадают в решето один раз или даже могут даже не попасть в решето вообще.

На данный момент программа получилась такая:

Используется синтаксис Java

class Primes {
        
        private static final int primesNumber = 85000000;
        private static final long[] primesList = new long[primesNumber];
        private static final int sieveSize = 30000;
        static {
                int primeIndex, testerIndex;
                long primeNominee, sieveStart, primeTester, doubleTester, primeLimit, testerSquare, tmpValue, sieveIndex;
                byte[] sieve = new byte[sieveSize];
                primesList[0] = 2;
                primesList[1] = 3;
                primeIndex = 2;
                primeNominee = 5;
                outerLoop1: while (true) {
                        testerIndex = 1;
                        primeLimit = (long) Math.sqrt(primeNominee);
                        while (true) {
                                primeTester = primesList[testerIndex];
                                if (primeLimit < primeTester) {
                                        primesList[primeIndex] = primeNominee;
                                        ++primeIndex;
                                        if (sieveSize + 2 < primeNominee * (primeNominee - 1)) {
                                                break outerLoop1;
                                        }
                                        break;
                                }
                                if (0 == primeNominee % primeTester) {
                                        break;
                                }
                                ++testerIndex;
                        }
                        primeNominee += 2;
                }
                sieveStart = primeNominee + 2;
                outerLoop2: while (true) {
                        for (sieveIndex = 0; sieveSize > sieveIndex; sieveIndex += 2) {
                                sieve[(int) sieveIndex] = 0;
                        }
                        testerIndex = 1;
                        while (true) {
                                primeTester = primesList[testerIndex];
                                testerSquare = primeTester * primeTester;
                                if (sieveStart + sieveSize <= testerSquare) {
                                        break;
                                }
                                doubleTester = 2 * primeTester;
                                if (sieveStart <= testerSquare) {
                                        sieveIndex = testerSquare - sieveStart;
                                } else {
                                        tmpValue = (sieveStart - primeTester) % doubleTester;
                                        if (0 == tmpValue) {
                                                sieveIndex = 0;
                                        } else {
                                                sieveIndex = doubleTester - tmpValue;
                                        }
                                }
                                while (sieveSize > sieveIndex) {
                                        sieve[(int) sieveIndex] = 1;
                                        sieveIndex += doubleTester;
                                }
                                ++testerIndex;
                        }
                        for (sieveIndex = 0; sieveSize > sieveIndex; sieveIndex += 2) {
                                if (0 == sieve[(int) sieveIndex]) {
                                        primesList[primeIndex] = sieveStart + sieveIndex;
                                        ++primeIndex;
                                        if (primesNumber == primeIndex) {
                                                break outerLoop2;
                                        }
                                }
                        }
                        sieveStart += sieveSize;
                }
        }
        
        public static long getPrime(int index) {
                return primesList[index];
        }
        
        public static long getLastPrime() {
                return primesList[primesNumber - 1];
        }
}

public class Primes_Eratosthenes {
        public static void main(String[] args) {
                int k, l, m;
                long startTime, stopTime;
                System.out.println(String.format("%13d", 2L * Integer.MAX_VALUE + 1));
                startTime = System.nanoTime();
                System.out.println(String.format("%13d", Primes.getLastPrime()));
                stopTime = System.nanoTime();
                System.out.println(String.format("\n%8.4f\n", 1e-9 * (stopTime - startTime)));
                m = 0;
                for (k = 0; 25 > k; ++k) {
                        for (l = 0; 7 > l; ++l) {
                                System.out.print(String.format("%5d", Primes.getPrime(m)));
                                ++m;
                        }
                        System.out.println();
                }
        }
}
 

Результат работы:

   4294967295
   1717783147

 14,6777

Можно сказать, что исходная задача практически решена.Не много не хватило памяти, надо увеличивать выделяемую под виртуальную машину объём. Ну и плюс я ещё хочу проделать различные оптимизации.

[Dmitriy40]

Какой-то у Вас первый цикл странный, перебор делителей, почему не решето? У меня он примерно такой (перевёл на С):

Используется синтаксис C

#define be1 65536 //До какого числа искать простые
#define be2 256 //Корень из be1
int p[be1];//Список простых
int np = 0;//Количество простых в p[]
for (int i = 2; i < be1; i++) p[i] = 1;
for (int i = 2; i < be1; i++) {
        if (p[i] == 1) {
                p[np++] = i;
                if (i < be2) for (int k =  i * i; k < be1; k += i) p[k] = 0;
        }
}

Вычисляет простые до 100 млн за 2 секунды, до миллиона за единицы мс. Ваш тоже конечно быстрый для десятка тысяч простых, но он как бы не проще этого и точно медленнее.

[B@R5uk]

Dmitriy40 в сообщении #1435594 писал(а):

Какой-то у Вас первый цикл странный, перебор делителей, почему не решето?

Потому что проходной вариант. Там ещё большой недостаток в том, что алгоритм фактически использует половину ячеек решета из-за того, что работает только с нечётными числами. Остальная половина присутствует, но не используется вообще. Сейчас вот думаю, что это можно и нужно устранить. Ну и само собой разумеется нужно ввести дополнительный массив для индексов вычёркивания.

-- 17.01.2020, 12:41 --

А ещё я думаю, как бы по-компактнее хранить получающиеся простые числа, сохранив при этом их функциональность. Потому что поиск простых чисел — это только вспомогательная задача, основная задача будет полученный массив в своей работе использовать. Идей пока две (вернее в сути в своей это она и та же идея). В обоих случаях простые числа будут храниться блоками некоторого размера: начало блока 8 байт, размер блока и разности по 1 байту на каждое простое. По одной идее можно хранить в байтах дифференциальные разности между соседними простыми числами. В этом случае блоки можно сделать по-длиннее. По второй идее можно хранить разности между простыми и началом блока. В этом случае функциональность лучше (особенно для поиска), но размеры блоков будут меньше и памяти потребуется больше.

Я даже в связи с этим изучил как растут промежутки между соседними простыми:

(Оффтоп)

         2            2     1
         3            3     2
         4            5     2
         5            7     4
         7           13     4
         9           19     4
        10           23     6
        12           31     6
        16           47     6
        17           53     6
        19           61     6
        22           73     6
        24           83     6
        25           89     8
        31          113    14
        63          293    14
        67          317    14
       100          523    18
       155          887    20
       190         1129    22
       218         1327    34
      1060         8467    34
      1184         9551    36
      1533        12853    36
      1664        14107    36
      1832        15683    44
      2226        19609    52
      2811        25471    52
      3386        31397    72
     14358       155921    86
     29041       338033    86
     30803       360653    96
     31546       370261   112
     40934       492113   114
    103521      1349533   118
    104072      1357201   132
    118506      1561919   132
    149690      2010733   148
    325853      4652353   154
    733589     11113933   154
    983016     15203977   154
   1094422     17051707   180
   1319946     20831323   210
   2850175     47326693   220
   6957877    122164747   222
  10539433    189695659   234
  10655463    191912783   248
  20684333    387096133   250
  22749705    428045491   250
  23163299    436273009   282
  64955635   1294268491   288
  72507381   1453168141   292

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