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

[S_WT]

По решету Эратосфена недавно такой код написал, то же чисто из интереса. Здесь никаких особых оптимизаций, классическое решето с просеиванием массива.

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

program eratosfen
implicit none
integer :: a(2:100), j, n
forall (j=2:100) a(j) = j

        n = 2   
        do while (n*n < 100)
           forall (j = n*n : 100 : n) a(j) = 0
           n = n + 1
                do while (a(n)==0)
                  n = n + 1
                end do
        end do

  n = count (a/=0)
  a(2:n+1) = pack (a, a/=0)

print '(10i5)', a(2:n+1)
end program eratosfen
 

Выдает:

    2    3    5    7   11   13   17   19   23   29
   31   37   41   43   47   53   59   61   67   71
   73   79   83   89   97

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

[B@R5uk]

$p_{6057}=59999$ Кто больше?

[B@R5uk]

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

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

import java.util.Arrays;

class Primes {
        
        private static final int sieveSize = 25000;             //      No less than 23167
        private static final int divisorsNumber = 4792;
        private static final int primesNumber = 105097565;
        private static final int[] primesList = new int[primesNumber];
        static {
                int primeIndex, sieveStart, divisorIndex, primeDivisor, sieveIndex;
                int[] offsetsList = new int[divisorsNumber];
                byte[] sieve = new byte[sieveSize];
                primesList[0] = 2;
                primeIndex = 1;
                sieveStart = 3;
                divisorIndex = 0;
                while (true) {
                        if (0 == sieve[divisorIndex]) {
                                primeDivisor = sieveStart + (divisorIndex << 1);
                                sieveIndex = (primeDivisor * primeDivisor - sieveStart) >> 1;
                                for (; sieveSize > sieveIndex; sieveIndex += primeDivisor) {
                                        sieve[sieveIndex] = 1;
                                }
                                offsetsList[primeIndex] = sieveIndex - sieveSize;
                                primesList[primeIndex] = primeDivisor;
                                ++primeIndex;
                                if (divisorsNumber == primeIndex) {
                                        break;
                                }
                        }
                        ++divisorIndex;
                }
                for (sieveIndex = divisorIndex + 1; sieveSize > sieveIndex; ++sieveIndex) {
                        if (0 == sieve[sieveIndex]) {
                                primesList[primeIndex] = sieveStart + (sieveIndex << 1);
                                ++primeIndex;
                        }
                }
                outerLoop: while (true) {
                        sieveStart += sieveSize << 1;
                        Arrays.fill(sieve, (byte) 0);
                        for (divisorIndex = 1; divisorsNumber > divisorIndex; ++divisorIndex) {
                                primeDivisor = primesList[divisorIndex];
                                sieveIndex = offsetsList[divisorIndex];
                                for (; sieveSize > sieveIndex; sieveIndex += primeDivisor) {
                                        sieve[sieveIndex] = 1;
                                }
                                offsetsList[divisorIndex] = sieveIndex - sieveSize;
                        }
                        for (sieveIndex = 0; sieveSize > sieveIndex; ++sieveIndex) {
                                if (0 == sieve[sieveIndex]) {
                                        primesList[primeIndex] = sieveStart + (sieveIndex << 1);
                                        ++primeIndex;
                                        if (primesNumber == primeIndex) {
                                                break outerLoop;
                                        }
                                }
                        }
                }
        }
        
        public static long getPrime(int index) {
                return primesList[index];
        }
        
        public static long getLastPrime() {
                return primesList[primesNumber - 1];
        }
}

public class Primes_Eratosthenes_2 {
        public static void main(String[] args) {
                int k, l, m;
                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\n", 1e-9 * (stopTime - startTime)));
                m = 99900;
                for (k = 0; 20 > k; ++k) {
                        for (l = 0; 5 > l; ++l) {
                                System.out.print(String.format("%8d", Primes.getPrime(m)));
                                ++m;
                        }
                        System.out.println();
                }
        }
}
 

   2147483647
   2147483647

  9,5549

При этом я решил пока ограничиться простыми числами в диапазоне $2...2^{31}-1$, что как раз соответствует максимальному значению int в Java (разработчики посчитали, что беззнаковые целые не нужны). Теперь я знаю, что int в Java заканчивается простым числом Мерсенна.

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

Так же я поигрался с размером сегмента решета. Эффект не значительный, где-то 5%-10%. Причём самое парадоксальное, уменьшение размера приводит к увеличению скорости работы. К сожалению, размер сегмента у меня ограничен снизу в связи с тем, что при обработке первого сегмента решета происходит поиск проверяющих простых делителей и инициализация массива смещений. Не знаю, актуален ли для Java вопрос использования кэшей процессора и связано ли это изменение во времени работы с ним. У меня, если что, Intel Core 2 Duo 2GHz P7350 с L1 кэшем 32KБ на каждое из двух ядер и L2 кэшем 3МБ.

[Dmitriy40]

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

и разности по 1 байту на каждое простое.

Учтите что разности между соседними простыми (очень полезные данные для оценки, сохраните и сверяйтесь) растут и для числа $436273009$ становятся больше $256$. Если хранить половину разности (всё равно она чётная), то переполнение байта будет на числе $304599508537$, которое к моему и вашему счастью в 140 раз больше $2^{31}$. Однако даже 300млрд не так уж много, это всего минут 10 работы хорошего генератора ... А все простые до $2^{31}$, которых $105097565$, займут 100МБ. И даже если их хранить в прямом виде, как int32, то это 400МБ, вполне терпимо для компьютера.

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

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

Тут вопрос как именно их использовать. Нужно произвольное обращение к списку или можно его пройти последовательно (пусть даже несколько раз, но не сотни). Если произвольное, то хранить массив (разностей или прямо, не суть), если можно пройти последовательно, то массив не хранить, обрабатывать по мере получения.
Если хранить разности, то можно завести дополнительный массив, в котором хранить точную величину каждого сотого (тысячного) простого, для ускорения прохода по разностям.
Если часты повторные обращения к близким простым, то при запросе одного соседние к нему просчитать и пересохранить в отдельный массив-кэш, из которого их вынуть значительно быстрее в следующий раз. Можно задействовать и стандартный тип хэштаблиц.
Но повторю, лучше заранее оценить время обработки и время генерации, возможно генерация займёт менее 1% общего времени и можно её перезапускать и не хранить все простые.
А ещё можно один раз их просчитать, сохранить в файл (уж 400МБ найти на диске не проблема вообще), а потом пользоваться отображением файла в память. Оно работает быстро, удобно и кэшируется самой ОС и не требует наличия свободной физической памяти такого же объёма, в отличии от чтения файла в память (или расчёта).

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

Причём самое парадоксальное, уменьшение размера приводит к увеличению скорости работы.

Странно, у меня тенденция наоборот, скорость растёт при увеличении размера сегмента. Но возможно у меня программа оптимальнее и тормозит в другом месте чем у Вас, потому и тенденция другая. Или для малых чисел тенденция меняется, что тоже возможно, я до $10^{15}$ вообще не оптимизировал счёт.

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

Не знаю, актуален ли для Java вопрос использования кэшей процессора

У Вас Java выполняется не на процессоре? А если на нём, то актуален. Но у вас сегмент и так меньше L1, ещё меньше делать имеет смысл только если скорость реально растёт.

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

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

Заведите для работы первого цикла отдельный массив, уж лишняя сотня КБ, тем более временно, совершенно не жалко.

[B@R5uk]

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

для числа $436273009$ становятся больше $256$.

Так я потому и писал, что числа надо хранить блоками: оригинальное значение простого числа в виде long + массив приращений в виде byte. Очевидно, что после каждого такого большого скачка будет начинаться новый блок. Тут скорее вопрос в другом: как много изолированных простых чисел? То есть таких простых чисел, которые отделены и от большего, и от меньшего соседа большим промежутком. Сейчас, наверное, прямо и займусь этим вопросом.

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

Но у вас сегмент и так меньше L1, ещё меньше делать имеет смысл только если скорость реально растёт.

Сегмент то да, но кроме сегмента ещё же есть массив смещений, и его размер — почти 19 КБ. И обращение к каждому его элементу происходит на каждой итерации цикла два раза. $19+23=42>32$ — не влезает. У меня, кстати, появились две идеи, как с этими смещениями и делителями можно поступить.

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

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

-- 17.01.2020, 19:40 --

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

Как много изолированных простых чисел?.. Сейчас, наверное, прямо и займусь этим вопросом.

Что-то такое получается:

(Оффтоп)

         2            3     1     2
         3            5     2     2
         4            7     2     4
         5           11     4     2
         6           13     2     4
         7           17     4     2
         8           19     2     4
         9           23     4     6
        12           37     6     4
        15           47     4     6
        16           53     6     6
        24           89     6     8
        37          157     6     6
        40          173     6     6
        47          211    12    12
       240         1511    12    12
       283         1847    16    14
       327         2179    18    24
       368         2503    26    18
       549         3967    20    22
      1337        11027    24    20
      1400        11657    24    20
      1663        14107    20    36
      1865        16033    26    24
      1939        16787    24    24
      2065        18013    24    28
      2700        24281    30    36
      3995        37747    30    34
      4059        38501    40    42
      5949        58831    42    58
     18281       203713    44    48
     18526       206699    48    50
     34871       413353    54    58
     98104      1272749    70    62
    162604      2198981    72    80
    314604      4479581    72    80
    355155      5102953    90    76
    722697     10938023   102    96
    826235     12623189   100    98
   1546100     24662467    98   102
   4258995     72546283   140   108
   8040878    142414669   116   190
   9123683    162821917   164   124
   9170830    163710121   150   136
  17567977    325737821   172   156
  57161802   1131241763   160   170
  88474734   1791752797   168   162

Формат вывода: номер простого числа - само число - расстояние до предыдущего соседа - расстояние до следующего соседа.

[Dmitriy40]

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

Так я потому и писал, что числа надо хранить блоками: оригинальное значение простого числа в виде long + массив приращений в виде byte.

Не люблю нерегулярных массивов. Лучше регулярный массив приращений и столь же регулярный массив точных значений. До 300млрд вам байтов хватит.
Кроме того, полезно знать формулу $\pi (N) \approx N/\ln N$, для чисел около $2^{31}$ она даёт средний интервал между простыми всего лишь 20-21, можно прикинуть среднюю частоту больших интервалов (вспомнив что средняя частота простых близнецов порядка $1/\ln^2 N$).
Вам лучше не рекорды интервалов посчитать, а гистограмму величины интервала, какие чаще, какие реже, и насколько. Потом думать.
Но мне всё равно кажется странным желание сэкономить 400МБ памяти или диска.

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

Сегмент то да, но кроме сегмента ещё же есть массив смещений, и его размер — почти 19 КБ. И обращение к каждому его элементу происходит на каждой итерации цикла два раза. $19+23=42>32$ — не влезает.

Забыли про сами простые, к ним тоже обращение постоянно. Но это не совсем корректный расчёт, к обоим массивам (и простых и смещений) обращение строго последовательное, а такое обращение быстрое (процессоры имеют блоки предвыборки последовательных данных и буфера отложенной записи) и практически не тратит объём кэша, ведь старые значения уже не нужны и потому не тормозят вычисления даже когда вытесняются из кэша по первому требованию. Даже объём самого сегмента решета тоже не обязан быть строго меньше объёма кэша, я когда игрался размерами, видел что даже тройное превышение размера L2 (но меньше L3) практически не уменьшало скорость.
Ещё наблюдение: упаковка флагов простоты в битовый массив давало примерно 3-х кратное (если не ошибаюсь) ускорение при той же длине сегмента в байтах, несмотря на усложнение операций во внутреннем цикле, увеличение длины проверяемого интервала в 8 раз пересиливало замедление операций в цикле.
Впрочем это всё лучше перепроверять на конкретном компьютере.

[Yadryara]

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

$p_{6057}=59999$ Кто больше?

В каком смысле, кто больше? По большому количеству девяток подряд соскучились? Так это, годовалый ребёнок знает, что $999999999...999999999997999...9999999$ — простое.

51 девятка, затем 7-ка или 4-ка, и ещё 44 девятки.

[Andrey A]

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

... как много изолированных простых чисел? То есть таких простых чисел, которые отделены и от большего, и от меньшего соседа большим промежутком.

Если $p$ простое, и $P=(p-1)!-1$ также простое, то вокруг $P$ большие промежутки. Сколь угодно большие, зависит только от величины $p$. Насчет годовалого ребенка не знаю, но в Кванте об этом писали.

[Dmitriy40]

Andrey A в сообщении #1435745 писал(а):

Если $p$ простое, и $P=(p-1)!-1$ также простое, то вокруг $P$ большие промежутки. Сколь угодно большие, зависит только от величины $p$.

Это слишком слабое утверждение, промежутки слишком медленно растут, это ещё если $P$ будет простым для больших факториалов. Вот для $p \le 200$:

? forprime(p=1,200, x=(p-1)!-1; if(!ispseudoprime(x),next); a=precprime(x-1); b=nextprime(x+1); printf("%u: %u-%u-%u\n",p,x-a,x,b-x))
5: 4-23-6
7: 10-719-8
13: 12-479001599-30
31: 78-265252859812191058636308479999999-108
167: 1038-9.0037e297-444

Последнюю строку сжал округлением. Слева и справа интервалы до ближайшего простого.

B@R5uk
Рекорды интервала для уединённых простых есть в A023186, сами рекордные интервалы в A023187 (странно что не синхронные, забыл видимо добавить).

[B@R5uk]

Dmitriy40, а вы эту статью видели? Сам читал запоем. Она правда древняя: 2017 год...

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

Не люблю нерегулярных массивов. Лучше регулярный массив приращений и столь же регулярный массив точных значений.

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

[Dmitriy40]

Yadryara в сообщении #1435741 писал(а):

Так это, годовалый ребёнок знает, что $999999999...999999999997999...9999999$ — простое.
51 девятка, затем 7-ка или 4-ка, и ещё 44 девятки.

А почему именно это число? С тем же количеством цифр есть ещё 10 простых. Наибольшее из них $10^{96}-1-10^{12}$.
И полно с ещё большим количеством девяток. Например 9 простых длиной 100 цифр из 99 девяток и одной другой цифры, наибольшее $10^{100}-1-10^4$.

[B@R5uk]

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

А почему именно это число?

Может потому что он его самостоятельно нашёл и проверил? Впрочем, всё всё равно интересное замечание. Я, например, не знал, что люди таким заморочились.

[Andrey A]

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

Это слишком слабое утверждение...

Мне казалось достаточное, чтобы расхотелось искать рекорды ))

[Dmitriy40]

B@R5uk
Нет, статью не видел. Вам там полезен как минимум Рис.13 и Рис.21, как раз статистика интервалов между простыми.
Меня мало интересует статистика самих простых, мне интереснее алгоритмы и программы их генерации.

-- 17.01.2020, 23:03 --

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

Я, например, не знал, что люди таким заморочились.

Я тоже не знал. Запустил PARI/GP и сварганил за пару минут программку проверки, она за секунду выдала десятки страниц чисел ... Ограничил одним значением количества цифр и отсортировал по возрастанию. Всего минут 5 работы. Чем PARI/GP и привлекателен, просто (но медленно!) работать с очень большими числами.
Программка:

? for(k=100,100,forstep(i=k-1,0,-1,forstep(d=8,1,-1, x=10^k-1-d*10^i; if(ispseudoprime(x),printf("10^%u-1-%u*10^%u\n",k,d,i)))))
10^100-1-4*10^83
10^100-1-8*10^75
10^100-1-7*10^62
10^100-1-1*10^47
10^100-1-1*10^40
10^100-1-8*10^10
10^100-1-4*10^8
10^100-1-5*10^5
10^100-1-1*10^4
time = 31 ms.

[B@R5uk]

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

ispseudoprime(x)

Я так понимаю, что это статистическая проверка, а не честный перебор всех простых делителей до корня? Просто я натыкался на статью Prime and Prejudice. Primality Testing Under Adversarial, в которой люди обсуждают как они конструировали составные числа, которые гарантированно или почти гарантированно проходили каждый из известных популярных тестов на простоту.