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

[Dmitriy40]

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

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

ispseudoprime(x)

Я так понимаю, что это статистическая проверка, а не честный перебор всех простых делителей до корня?

Разумеется не перебор, что Вы! Она по умолчанию использует BPSW тест. Почитайте про тесты простоты по ссылкам оттуда, Вам будет интересно.
Помнится ни одного исключения (ошибки этой функции) не известно:

вики писал(а):

Мощность теста состоит в том, что не существует известных пересечений списков псевдопростых чисел Ферма и псевдопростых чисел Люка.
Ни одно составное число, меньшее $2^{64}$ (что примерно равно $1{,}845 \cdot 10^{19}$), не проходит тест БПСВ[5]. Таким образом, этот тест можно считать детерминированным тестом простоты для чисел, меньших указанной границы. Также пока не известно ни одно составное число, большее границы, которое проходит тест.

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

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

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

[Yadryara]

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

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

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

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

B@R5uk, совсем необязательно о присутствующем говорить в 3-м лице.

Самостоятельно не находил, равно как и другое большое число в вашей же недавней теме. Вы, видимо, не только об этом не знали, но и о многом-многом другом. Как раз для того, чтоб люди больше знали и могли заниматься поиском нового, а не переоткрытием старого велосипеда была создана тема «Сделаем проверку по OEIS обязательной?»

В настоящее время в OEIS примерно $330\;000$ последовательностей. Вдумайтесь в это число. Если вы хотите сделать что-то новое, будьте добры тщательно ознакомиться с уже существующими результатами и алгоритмами.

В данном же случае я просто заглянул в b-file A105975 и взял оттуда 2 последних числа. И проверил их на простоту в Альфе.

[B@R5uk]

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

...заглянул в b-file A105975...

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

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

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

Последовал этому совету, что позволило мне поиграться с размером сегмента решета. Оптимальным оказался размер в диапазоне 20...30 тысяч. Точное значение установить сложно, так как время работы программы от раза к разу скачет. Измерение скорости в зависимости от размера сегмента в этом диапазоне сопоставимо со случайным разбросом. Уменьшение величины сегмента на порядок от оптимального снижает быстродействие в 1,5...2 раза, а увеличение на 2 порядка — на 30%.

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

С проверкой ещё ладно, а вот в инициализации придётся использовать взятие остатка. Какого будет ускорение и будет ли вообще — вопрос.

Я ошибся: вся инициализация заключается в вычитании. А вот сомнение в ускорении оказалось обоснованным. Причём реально ускорение ещё меньше, чем я ожидал — всего порядка 7% (какая жалость...). Решето с постепенным добавлением проверяющих делителей по мере работы получилось такое:

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

import java.util.Arrays;

class Primes {
        
        private static final int sieveSize = 26000;
        private static final int prepSieveSize = 23168;
        private static final int divisorsNumber = 4792;
        private static final int primesNumber = 105097565;
        private static final int[] primesList = new int[primesNumber];
        static {
                int primeIndex, nomineeIndex, sieveIndex, divisorIndex, limitIndex;
                int sieveStart, sieveEnd, primeDivisor;
                long divisorSquare;
                int[] offsetsList = new int[divisorsNumber];
                byte[] sieve;
                primesList[0] = 2;
                primeIndex = 1;
                sieveStart = 3;
                nomineeIndex = 0;
                sieve = new byte[prepSieveSize];
                sieveEnd = sieveStart + (prepSieveSize << 1);
                //      Preparatory sieve loop for testing prime divisors
                while (true) {
                        if (0 == sieve[nomineeIndex]) {
                                primeDivisor = sieveStart + (nomineeIndex << 1);
                                divisorSquare = (long) primeDivisor * primeDivisor;
                                if (sieveEnd <= divisorSquare) {
                                        limitIndex = primeIndex;
                                        break;
                                }
                                sieveIndex = (int) ((divisorSquare - sieveStart) >> 1);
                                for (; prepSieveSize > sieveIndex; sieveIndex += primeDivisor) {
                                        sieve[sieveIndex] = 1;
                                }
                                offsetsList[primeIndex] = sieveIndex - prepSieveSize;
                                primesList[primeIndex] = primeDivisor;
                                ++primeIndex;
                        }
                        ++nomineeIndex;
                }
                //      Collecting remaining testing divisors
                while (true) {
                        if (0 == sieve[nomineeIndex]) {
                                primesList[primeIndex] = sieveStart + (nomineeIndex << 1);
                                ++primeIndex;
                                if (divisorsNumber == primeIndex) {
                                        break;
                                }
                        }
                        ++nomineeIndex;
                }
                sieveStart = sieveEnd;
                sieve = new byte[sieveSize];
                //      Main sieve loop
                outerLoop: while (true) {
                        sieveEnd = sieveStart + (sieveSize << 1);
                        //      Striking out composites with already used testing divisors
                        for (divisorIndex = 1; limitIndex > divisorIndex; ++divisorIndex) {
                                primeDivisor = primesList[divisorIndex];
                                sieveIndex = offsetsList[divisorIndex];
                                for (; sieveSize > sieveIndex; sieveIndex += primeDivisor) {
                                        sieve[sieveIndex] = 1;
                                }
                                offsetsList[divisorIndex] = sieveIndex - sieveSize;
                        }
                        //      Striking out prime squares with new testing divisors
                        while (true) {
                                primeDivisor = primesList[divisorIndex];
                                divisorSquare = (long) primeDivisor * primeDivisor;
                                if (sieveEnd <= divisorSquare) {
                                        limitIndex = divisorIndex;
                                        break;
                                }
                                sieveIndex = (int) ((divisorSquare - sieveStart) >> 1);
                                for (; sieveSize > sieveIndex; sieveIndex += primeDivisor) {
                                        sieve[sieveIndex] = 1;
                                }
                                offsetsList[divisorIndex] = sieveIndex - sieveSize;
                                ++divisorIndex;
                        }
                        //      Collecting primes
                        for (nomineeIndex = 0; sieveSize > nomineeIndex; ++nomineeIndex) {
                                if (0 == sieve[nomineeIndex]) {
                                        primesList[primeIndex] = sieveStart + (nomineeIndex << 1);
                                        ++primeIndex;
                                        if (primesNumber == primeIndex) {
                                                break outerLoop;
                                        }
                                }
                        }
                        sieveStart = sieveEnd;
                        Arrays.fill(sieve, (byte) 0);
                }
        }
        
        public static int getPrime(int index) {
                return primesList[index];
        }
        
        public static int getLastPrime() {
                return primesList[primesNumber - 1];
        }
}

public class Primes_Eratosthenes_3 {
        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();
                }
                //*/
        }
}
 

Осталось ещё протестить идею битовой упаковки флагов простоты. Но в целом, идею по ускорению работы на данном этапе у меня подходят к концу.

[Dmitriy40]

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

В данном же случае я просто заглянул в b-file A105975 и взял оттуда 2 последних числа.

Понятно, спасибо. Абсолютно правильное решение. Поддерживаю.
К большому сожалению в OEIS поиск по b-files не работает, так что далеко не всё там присутствующее можно найти встроенными средствами (без гугла).

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

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

Чтобы скакало меньше оберните всю программу (вычислительную часть) во внешний цикл раз 10-100-1000, чтобы время стало порядка десятков-сотен секунд, тогда скачки будут в доли секунды, т.е. менее процента.

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

Уменьшение величины сегмента на порядок от оптимального снижает быстродействие в 1,5...2 раза, а увеличение на 2 порядка — на 30%.

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

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

[B@R5uk]

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

Конкретно с этим вопросом, конечно, оказалось всё просто. Если взять, например, 4000-й простой делитель (для моего диапазона их всего 4792), то его значение $p_{4000}=37813$. А следующий за ним: $p_{4001}=37831$, расстояние равно 18. Поскольку простой делитель подключается в работу, когда его квадрат попадает в сегмент, то всё дело в расстоянии между квадратами. Его можно оценить как удвоенное значение делителя умножить на расстояние до него, то есть $2\cdot 18\cdot37800=1360800$. Эта величина как раз соответствует 26 эффективным размерам сегмента решета: $\lfloor 1360800 \ (2\cdot 26000)\rfloor=26$. В среднем же статистика расстояний в диапазоне номеров делителей от 3000 до 4000 такая:

     2   134  13,40%
     4   118  11,80%
     6   208  20,80%
     8    76   7,60%
    10    98   9,80%
    12   107  10,70%
    14    40   4,00%
    16    35   3,50%
    18    60   6,00%
    20    27   2,70%
    22    32   3,20%
    24    11   1,10%
    26    11   1,10%
    28     6   0,60%
    30    17   1,70%
Others    20   2,00%
10.374

В процессе я осознал один интересный факт: для проверяющих простых делителей в диапазоне $\left[ \sqrt[3]{L},\sqrt{L} \right]$ (где L — это верхняя граница искомых простых чисел) составные числа, эффективно вычёркиваемые этими делителями (то есть не вычёркиваемые делителями до $\sqrt[3]{L}$) являются произведениями двух простых чисел, одно из которых — проверяющий делитель $p_n$, а номер m второго простого числа лежит на отрезке $\left[ n,\pi \left( \frac{L}{{{p}_{n}}} \right)\right]$, где $\pi(x)$ — это та самая функция распределения простых чисел. Имея таблицу этих чисел её, кстати, очень легко посчитать:

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

        public static int findPi(int value) {
                int result = Arrays.binarySearch(primesList, value);
                if (0 <= result) {
                        return  1 + result;
                } else {
                        return -1 - result;
                }
        }
 

Это знание предоставляет практическую возможность посчитать количество составных чисел, вычёркиваемых простым делителем $p_n$. Это число равно $\pi \left( \frac{L}{{{p}_{n}}} \right)-n+1$. А просуммировав эти значения в диапазоне от N до 4792 можно посчитать полное число составных чисел, вычёркиваемых простыми делителями начиная с $p_N$.

Однако всё познаётся в сравнении, просто числа ничего не дадут. Этот результат, очевидно, надо сравнивать с количеством обращений ко всем простым делителям от N до 4792, которые производятся в цикле программы. При обработке каждого сегмента решета происходит одно обращение к работающим простым делителям. Делитель $p_n$ начинает работать, когда его квадрат попадает в сегмент, то есть начиная с сегмента $\left\lfloor \frac{p_{n}^{2}-\text{46339}}{2\cdot 26000} \right\rfloor$. Номер же последнего сегмента: $\left\lfloor \frac{{{2}^{31}}-\text{46340}}{2\cdot 26000} \right\rfloor =\text{412968}$. Для расчёта, разумеется, составил программу. Результаты получились такие:

(Оффтоп)

    1300    10657 19588199   84106010
    1325    10891 19172744   83129352
    1350    11149 18766970   82155278
    1375    11393 18370438   81183785
    1400    11657 17982999   80215010
    1425    11897 17605026   79249193
    1450    12109 17234555   78285889
    1475    12347 16871818   77325277
    1500    12553 16516433   76367412
    1525    12791 16167563   75412024
    1550    13007 15825727   74459611
    1575    13249 15490199   73509979
    1600    13499 15161580   72563515
    1625    13729 14839796   71620377
    1650    13967 14523681   70680103
    1675    14251 14214156   69743313
    1700    14519 13911752   68810375
    1725    14737 13614964   67880676
    1750    14947 13323150   66953923
    1775    15199 13036715   66030741
    1800    15401 12755385   65110839
    1825    15647 12478934   64194260
    1850    15881 12207365   63281054
    1875    16103 11940859   62371524
    1900    16381 11679366   61465896
    1925    16649 11422920   60564471
    1950    16903 11171500   59667164
    1975    17137 10924564   58773872
    2000    17389 10682228   57884809
    2025    17609 10443942   56999478
    2050    17891 10210092   56118403
    2075    18119  9980464   55241521
    2100    18313  9754509   54368526
    2125    18553  9532235   53499471
    2150    18899  9314372   52635426
    2175    19181  9101047   51777015
    2200    19423  8891762   50923591
    2225    19609  8685562   50073998
    2250    19891  8483149   49229344
    2275    20117  8284087   48389286
    2300    20357  8088241   47553601
    2325    20627  7895764   46722727
    2350    20897  7706926   45897360
    2375    21139  7521172   45076954
    2400    21383  7338763   44261489
    2425    21599  7159456   43451136
    2450    21841  6983014   42645391
    2475    22079  6809588   41844873
    2500    22307  6638820   41049142
    2525    22613  6471165   40259128
    2550    22817  6306574   39474576
    2575    23059  6144823   38695502
    2600    23321  5985799   37921428
    2625    23599  5829913   37153689
    2650    23827  5676601   36391371
    2675    24049  5525922   35634282
    2700    24281  5377703   34882099
    2725    24571  5232190   34136608
    2750    24877  5089643   33398158
    2775    25147  4949824   32666293
    2800    25391  4812487   31940631
    2825    25643  4677659   31221313
    2850    25913  4545162   30508109
    2875    26161  4415015   29801224
    2900    26399  4287302   29100580
    2925    26683  4162031   28406475
    2950    26891  4039072   27718586
    2975    27127  3918333   27036769
    3000    27449  3800013   26362472
    3025    27739  3684132   25696027
    3050    27947  3570444   25035790
    3075    28211  3458806   24382035
    3100    28499  3349478   23736149
    3125    28687  3242190   23096748
    3150    28933  3136770   22463187
    3175    29207  3033570   21837273
    3200    29443  2932333   21218335
    3225    29753  2833169   20606887
    3250    30059  2736279   20004027
    3275    30293  2641575   19408985
    3300    30559  2548957   18821677
    3325    30841  2458459   18242580
    3350    31091  2369714   17670911
    3375    31319  2282753   17106505
    3400    31601  2197550   16549619
    3425    31883  2114295   16001220
    3450    32159  2033207   15462035
    3475    32377  1953905   14930395
    3500    32609  1876115   14405498
    3525    32887  1800123   13888378
    3550    33113  1725967   13379383
    3575    33377  1653498   12878439
    3600    33613  1582599   12385608
    3625    33851  1513263   11899818
    3650    34157  1445829   11422737
    3675    34381  1380119   10954953
    3700    34649  1315954   10495166
    3725    34897  1253417   10043646
    3750    35159  1192633    9601563
    3775    35447  1133436    9168594
    3800    35759  1075882    8744772
    3825    36007  1020011    8331392
    3850    36277   965702    7926616
    3875    36551   913173    7531620
    3900    36781   862124    7145349
    3925    37013   812533    6767099
    3950    37309   764549    6398110
    3975    37547   718142    6038990
    4000    37813   673209    5688135
    4025    38119   630072    5348258
    4050    38371   588567    5019038
    4075    38677   548729    4700636
    4100    38921   510327    4391478
    4125    39181   473357    4092163
    4150    39419   437736    3801957
    4175    39719   403691    3522098
    4200    39971   371045    3252402
    4225    40231   339837    2992706
    4250    40531   310148    2744584
    4275    40829   281996    2507918
    4300    41081   255214    2281497
    4325    41333   229738    2064962
    4350    41603   205651    1859001
    4375    41851   182971    1662853
    4400    42073   161699    1476770
    4425    42331   141690    1300552
    4450    42557   122849    1133606
    4475    42793   105274     976807
    4500    43051    89005     829984
    4525    43399    74186     695341
    4550    43661    60862     574062
    4575    43969    48926     464536
    4600    44201    38319     365933
    4601    44203    37921     362206
    4602    44207    37524     358483
    4603    44221    37128     354767
    4604    44249    36735     351074
    4605    44257    36346     347429
    4606    44263    35959     343798
    4607    44267    35574     340177
    4608    44269    35190     336563
    4609    44273    34807     332952
    4610    44279    34425     329348
    4611    44281    34044     325754
    4612    44293    33665     322163
    4613    44351    33289     318593
    4614    44357    32920     315122
    4615    44371    32552     311661
    4616    44381    32188     308224
    4617    44383    31826     304804
    4618    44389    31465     301387
    4619    44417    31106     297981
    4620    44449    30750     294622
    4621    44453    30397     291318
    4622    44483    30047     288021
    4623    44491    29700     284775
    4624    44497    29355     281543
    4625    44501    29011     278321
    4626    44507    28669     275106
    4627    44519    28328     271901
    4628    44531    27990     268717
    4629    44533    27653     265553
    4630    44537    27317     262393
    4631    44543    26983     259240
    4632    44549    26650     256097
    4633    44563    26318     252964
    4634    44579    25989     249855
    4635    44587    25663     246774
    4636    44617    25338     243706
    4637    44621    25016     240690
    4638    44623    24696     237681
    4639    44633    24377     234675
    4640    44641    24061     231686
    4641    44647    23747     228711
    4642    44651    23434     225746
    4643    44657    23122     222788
    4644    44683    22812     219841
    4645    44687    22505     216938
    4646    44699    22199     214042
    4647    44701    21895     211167
    4648    44711    21592     208295
    4649    44729    21290     205440
    4650    44741    20992     202616
    4651    44753    20695     199813
    4652    44771    20399     197031
    4653    44773    20107     194280
    4654    44777    19816     191532
    4655    44789    19527     188791
    4656    44797    19241     186071
    4657    44809    18957     183364
    4658    44819    18675     180678
    4659    44839    18395     178009
    4660    44843    18118     175375
    4661    44851    17842     172748
    4662    44867    17568     170134
    4663    44879    17296     167548
    4664    44887    17026     164983
    4665    44893    16758     162432
    4666    44909    16492     159891
    4667    44917    16228     157378
    4668    44927    15965     154878
    4669    44939    15705     152396
    4670    44953    15448     149934
    4671    44959    15194     147497
    4672    44963    14941     145070
    4673    44971    14689     142650
    4674    44983    14438     140244
    4675    44987    14190     137858
    4676    45007    13944     135479
    4677    45013    13700     133135
    4678    45053    13459     130801
    4679    45061    13222     128537
    4680    45077    12987     126287
    4681    45083    12755     124064
    4682    45119    12525     121852
    4683    45121    12300     119702
    4684    45127    12076     117556
    4685    45131    11854     115420
    4686    45137    11633     113291
    4687    45139    11414     111172
    4688    45161    11196     109057
    4689    45179    10981     106980
    4690    45181    10769     104934
    4691    45191    10558     102892
    4692    45197    10350     100867
    4693    45233    10143      98853
    4694    45247     9942      96901
    4695    45259     9742      94974
    4696    45263     9544      93067
    4697    45281     9347      91167
    4698    45289     9153      89299
    4699    45293     8961      87445
    4700    45307     8771      85598
    4701    45317     8583      83775
    4702    45319     8397      81970
    4703    45329     8213      80168
    4704    45337     8031      78383
    4705    45341     7850      76612
    4706    45343     7671      74848
    4707    45361     7493      73088
    4708    45377     7318      71359
    4709    45389     7145      69658
    4710    45403     6974      67978
    4711    45413     6806      66323
    4712    45427     6641      64685
    4713    45433     6479      63071
    4714    45439     6319      61468
    4715    45481     6160      59875
    4716    45491     6005      58356
    4717    45497     5852      56854
    4718    45503     5700      55363
    4719    45523     5549      53882
    4720    45533     5400      52436
    4721    45541     5252      51008
    4722    45553     5106      49594
    4723    45557     4964      48201
    4724    45569     4823      46815
    4725    45587     4685      45450
    4726    45589     4551      44117
    4727    45599     4418      42787
    4728    45613     4286      41475
    4729    45631     4157      40187
    4730    45641     4029      38931
    4731    45659     3904      37692
    4732    45667     3782      36485
    4733    45673     3661      35292
    4734    45677     3541      34109
    4735    45691     3423      32933
    4736    45697     3306      31782
    4737    45707     3191      30642
    4738    45737     3078      29519
    4739    45751     2967      28449
    4740    45757     2857      27404
    4741    45763     2749      26369
    4742    45767     2642      25345
    4743    45779     2536      24328
    4744    45817     2432      23332
    4745    45821     2332      22403
    4746    45823     2234      21481
    4747    45827     2137      20562
    4748    45833     2042      19650
    4749    45841     1948      18749
    4750    45853     1856      17862
    4751    45863     1765      16996
    4752    45869     1677      16148
    4753    45887     1591      15310
    4754    45893     1509      14504
    4755    45943     1428      13709
    4756    45949     1353      13002
    4757    45953     1279      12306
    4758    45959     1206      11617
    4759    45971     1135      10938
    4760    45979     1066      10281
    4761    45989      998       9638
    4762    46021      932       9012
    4763    46027      871       8443
    4764    46049      812       7885
    4765    46051      757       7366
    4766    46061      704       6850
    4767    46073      652       6352
    4768    46091      602       5875
    4769    46093      554       5430
    4770    46099      508       4989
    4771    46103      464       4558
    4772    46133      421       4134
    4773    46141      382       3764
    4774    46147      345       3408
    4775    46153      309       3062
    4776    46171      274       2727
    4777    46181      241       2424
    4778    46183      211       2139
    4779    46187      182       1857
    4780    46199      155       1582
    4781    46219      130       1329
    4782    46229      108       1111
    4783    46237       88        911
    4784    46261       71        725
    4785    46271       57        582
    4786    46273       44        457
    4787    46279       33        335
    4788    46301       23        224
    4789    46307       16        152
    4790    46309       10         91
    4791    46327        5         33
    4792    46337        1          8

В первом столбце находится номер простого числа, с которого начинается подсчёт, во втором — само простое число, в третьем — полное количество вычёркиваемых составных чисел и в четвёртом — количество обращений к простому числу, которое происходит в цикле программы. Если найти отношение данных в 4-ом столбце к данным в 3-ем, то узнаем среднее число обращений к простому делителю на одно вычёркивание составного числа. Получается такой интересный график:


Как видно, пока не так уж всё и плохо. Не больше десятка обращений (надо учитывать, что это интегральная картинка, фактически надо смотреть на самую первую точку). Каждое обращение к простому числу состоит из чтения из массива, сравнение с переменной, вычитания и записи в массив, если не было вычёркивания. Если вычёркивание было, то добавляются ещё одно сравнение, сложение и запись в массив решета. Я не знаю, как точно считать сложность этих действий, но получается 4-7 операций на обращение.

Всё это я вычислял применительно к идее реализации вычёркивания через очередь с приоритетом. Все простые делители, начиная с некоторого $p_N$, достаточно редко участвуют в вычёркивании. Поскольку вычёркиваемые составные числа являются произведением простого делителя $p_n$ на другой делитель $p_m\ge p_n$, то каждое такие составное число, соответствующее делителю $p_n$ для $n\in \left[ N,\pi \left( \sqrt{L} \right) \right]$ начиная с $p_n^2$ можно добавить в очередь с приоритетом и извлекать по одному при необходимости. Извлечённое составное число $p_n p_m$, соответствующее делителю $p_n$, вычёркивается из решета, а вместо него в очередь добавляется число $p_n p_{m+1}$ и работа продолжается. При этом $N>\pi \left( \sqrt[3]{L} \right)$, поскольку для меньших простых делителей второй множитель вычёркиваемого числа уже может быть составным.

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

[B@R5uk]

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

Есть другая методика, пользоваться не смещением в сегменте, а просто истинной величиной числа

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

То есть будут ли две операции ниже выполнятся одинаково быстро?

array[index] = data;

array[index + subindex] = data;

[Dmitriy40]

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

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

На архитектуре x86/x64 (а другие я знаю значительно хуже) совсем уж точно одинаково не будут, хоть и есть адресация с тремя операндами, но такой адрес вычисляется чуточку медленнее двух- и одно-операндного. На самом деле даже вдвое медленнее (1 такт против 0.5 такта), плюс ещё и параллелится с другими операциями хуже.
Но дело в том, что операция вычисления адреса лишь одна из нескольких, и потому даже замедление её аж вдвое не приведёт к заметному замедлению всего цикла. Гораздо больше влияния оказывает произвольный доступ к памяти сегмента для записи флага составного числа, латентность кэша L1 по чтению порядка 3 тактов (по записи нулевая из-за наличия буферов отложенной записи, но они не бесконечные, да и очередь запросов на чтение кэшлинии из L2 в L1 тоже не слишком большая (даже не знаю какая)). Или ошибки предсказания ветвлений (вычёркивать текущее число или нет). На всём этом фоне лишние полтакта теряются.

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

То есть будут ли две операции ниже выполнятся одинаково быстро?

На ассемблере — практически да (см. выше), что скомпилит ваш компилятор ЯВУ — вопрос, надо смотреть конкретно.

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

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

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

Другое дело что адресации [reg-reg] в процессоре нет и значит не получится прямо выполнить код sieve[a-base]=0, но операция вычитания занимает максимум 1 такт (а в потоке вплоть до 0.25 такта), и потому для малых чисел (менее длины сегмента) лучше сначала вычесть и получить смещение в сегменте, потом всё вычеркнуть, а потом добавить (тоже за 0.25-1 такт) обратно для восстановления правильной величины. Для больших чисел (более длины сегмента) вычитание возможно совместится со сравнением и дополнительного времени не потребует, а добавление для восстановления вообще не нужно.

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

[B@R5uk]

Dmitriy40, спасибо за разъяснение.

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

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

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

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

Можно и не к каждой, а лишь к первой и суммирование после последней

[Dmitriy40]

B@R5uk
Можно вообще по разному хранить первые смещения (для простых меньше длины сегмента) и остальные, первые лишь как смещение внутри сегмента, вторые прямо, ведь это лишь вопрос их инициализации и не более, примерно так (цикл вычёркивания):

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

for (i = 1; i < prime_small_end; i++) {//Цикл лишь по простым меньше длины сегмента
        for (k = d[i]; k < seg_length; k += prime[i]) seg[k] = 0;
        d[i] = k - seg_length;
}
for (; i < prime_end; i++)//Цикл по оставшимся простым
        if (k = d[i] - seg_start < seg_length) {
                seg[k] = 0;
                d[i] += prime[i];
        }

Если же хранить все одинаково, в прямом виде, то первый цикл будет таким:

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

for (i = 1; i < prime_small_end; i++) {//Цикл лишь по простым меньше длины сегмента
        for (k = d[i] - seg_start; k < seg_length; k += prime[i]) seg[k] = 0;
        d[i] = k + seg_start;
}

Думаю и то и другое будет работать одинаково по скорости.
И не надо никаких сложных методов адресации.

PS. Все константы рассчитываются заранее, до цикла вычёркивания.
PS. Написал из головы, без проверки, лишь для иллюстрации идеи. Но работать должно, наверное ...

-- 19.01.2020, 21:17 --

Операцию чтения prime[i] в k += prime[i] в первом цикле надо тоже выносить за цикл и в цикле добавлять уже просто переменную, но с этим должен справиться и компилятор сам.

[B@R5uk]

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

Можно вообще по разному хранить первые смещения (для простых меньше длины сегмента) и остальные

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

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

Сейчас вот это попробую:

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

Можно и не к каждой, а лишь к первой и суммирование после последней

Попробовал. Стало значительно хуже!

То есть из вот этих двух вариантов кода:

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

        divisorOffset = offsetsList[divisorIndex];
        if (sieveHalfEnd > divisorOffset) {
                primeDivisor = primesList[divisorIndex];
                do {
                        sieve[divisorOffset - sieveHalfStart] = 1;
                        divisorOffset += primeDivisor;
                } while (sieveHalfEnd > divisorOffset);
                offsetsList[divisorIndex] = divisorOffset;
        }
 

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

        divisorOffset = offsetsList[divisorIndex];
        if (sieveHalfEnd > divisorOffset) {
                primeDivisor = primesList[divisorIndex];
                divisorOffset -= sieveHalfStart;
                do {
                        sieve[divisorOffset] = 1;
                        divisorOffset += primeDivisor;
                } while (sieveSize > divisorOffset);
                offsetsList[divisorIndex] = divisorOffset + sieveHalfStart;
        }
 

именно первый выполняется быстрее. Причём я заметил, что со вторым вариантом кода время выполнения программы скачет значительно сильнее (до 10%). Без понятия почему так происходит. Прям вообще чудеса.

На самом деле я покривил душой, когда сказал, что ускорения не заметно. Если собрать статистику, то видно, что прибавка к скорости составляет 2,5%:

9.3049      9.5430
9.2755      9.5157
9.2793      9.5881
9.3070      9.5450
9.1764      9.4593
9.2855      9.5012
9.2812      9.5025
9.3159      9.6040
9.2288      9.3900
9.2664      9.5121
========================
9.2721      9.5161
0.0416      0.0615

Здесь в первом столбце время работы новой программы, во втором — старой.

Код новой программы такой:

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

import java.util.Arrays;

class Primes {
        
        private static final int sieveSize = 26000;
        private static final int prepSieveSize = 23168;
        private static final int divisorsNumber = 4792;
        private static final int primesNumber = 105097565;
        private static final int[] primesList = new int[primesNumber];
        static {
                int primeIndex, nomineeIndex, divisorOffset, divisorIndex, limitIndex;
                int sieveHalfStart, sieveHalfEnd, primeDivisor;
                int[] offsetsList = new int[divisorsNumber + 1];
                byte[] sieve;
                primesList[0] = 2;
                primeIndex = 1;
                sieveHalfStart = 1;
                nomineeIndex = 0;
                sieve = new byte[prepSieveSize];
                sieveHalfEnd = sieveHalfStart + prepSieveSize;
                //      Preparatory sieve loop for testing prime divisors
                while (true) {
                        if (0 == sieve[nomineeIndex]) {
                                primeDivisor = (sieveHalfStart + nomineeIndex << 1) + 1;
                                divisorOffset = primeDivisor * primeDivisor >> 1;
                                if (sieveHalfEnd <= divisorOffset) {
                                        limitIndex = primeIndex;
                                        break;
                                }
                                for (; sieveHalfEnd > divisorOffset; divisorOffset += primeDivisor) {
                                        sieve[divisorOffset - sieveHalfStart] = 1;
                                }
                                offsetsList[primeIndex] = divisorOffset;
                                primesList[primeIndex] = primeDivisor;
                                ++primeIndex;
                        }
                        ++nomineeIndex;
                }
                //      Collecting remaining testing divisors
                while (true) {
                        if (0 == sieve[nomineeIndex]) {
                                primeDivisor = (sieveHalfStart + nomineeIndex << 1) + 1;
                                primesList[primeIndex] = primeDivisor;
                                offsetsList[primeIndex] = primeDivisor * primeDivisor >> 1;
                                ++primeIndex;
                                if (divisorsNumber == primeIndex) {
                                        break;
                                }
                        }
                        ++nomineeIndex;
                }
                offsetsList[divisorsNumber] = Integer.MAX_VALUE;
                sieveHalfStart = sieveHalfEnd;
                sieve = new byte[sieveSize];
                //      Main sieve loop
                outerLoop: while (true) {
                        sieveHalfEnd = sieveHalfStart + sieveSize;
                        //      Striking out composites with already used testing divisors
                        for (divisorIndex = 1; limitIndex > divisorIndex; ++divisorIndex) {
                                divisorOffset = offsetsList[divisorIndex];
                                if (sieveHalfEnd > divisorOffset) {
                                        primeDivisor = primesList[divisorIndex];
                                        do {
                                                sieve[divisorOffset - sieveHalfStart] = 1;
                                                divisorOffset += primeDivisor;
                                        } while (sieveHalfEnd > divisorOffset);
                                        offsetsList[divisorIndex] = divisorOffset;
                                }
                        }
                        //      Striking out prime squares with new testing divisors
                        while (true) {
                                divisorOffset = offsetsList[divisorIndex];
                                if (sieveHalfEnd > divisorOffset) {
                                        primeDivisor = primesList[divisorIndex];
                                        do {
                                                sieve[divisorOffset - sieveHalfStart] = 1;
                                                divisorOffset += primeDivisor;
                                        } while (sieveHalfEnd > divisorOffset);
                                        offsetsList[divisorIndex] = divisorOffset;
                                } else {
                                        limitIndex = divisorIndex;
                                        break;
                                }
                                ++divisorIndex;
                        }
                        //      Collecting primes
                        for (nomineeIndex = 0; sieveSize > nomineeIndex; ++nomineeIndex) {
                                if (0 == sieve[nomineeIndex]) {
                                        primesList[primeIndex] = (sieveHalfStart + nomineeIndex << 1) + 1;
                                        ++primeIndex;
                                        if (primesNumber == primeIndex) {
                                                break outerLoop;
                                        }
                                }
                        }
                        sieveHalfStart = sieveHalfEnd;
                        Arrays.fill(sieve, (byte) 0);
                }
        }
        
        public static int getPrimesNumber() {
                return primesNumber;
        }
        
        public static int getLastPrime() {
                return primesList[primesNumber - 1];
        }
        
        public static int getPrime(int index) {
                return primesList[index];
        }
        
        public static int findPi(int value) {
                int result = Arrays.binarySearch(primesList, value);
                if (0 <= result) {
                        return  1 + result;
                } else {
                        return -1 - result;
                }
        }
}

public class Primes_Eratosthenes_4 {
        public static void main(String[] args) {
                displayTime();
                displayA007053();
                //System.out.println(46349 * 46349 >> 1 & 0x7FFFFFFF);
        }
        
        private static void displayTime() {
                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)));
        }
        
        private static void displayA007053() {
                int k, lastPrime, limit;
                for (k = 0; 3 > k; ++k) {
                        System.out.println(String.format("%6d %d", k, Primes.findPi(k + 1)));
                }
                limit = 3;
                lastPrime = Primes.getLastPrime();
                while (true) {
                        limit = (limit << 1) + 1;
                        if (lastPrime < limit) {
                                break;
                        }
                        System.out.println(String.format("%6d %d", k, Primes.findPi(limit)));
                        ++k;
                        if (lastPrime == limit) {
                                break;
                        }
                }
        }
}
 

[Dmitriy40]

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

именно первый выполняется быстрее.

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

[B@R5uk]

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

Вы уверены что разница в скорости значительно превышает погрешность измерения?

Да! Во втором случае разброс большой, но нижняя граница времени больше среднего для первого случая.

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

но в Java то вроде как JIT-компиляция...

Я вот тут думаю: кроме всего прочего время работы моей программы так же включает время работы JIT-компилятора? Интересно бы знать, в каких случая он запускается.

[Dmitriy40]

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

Интересно бы знать, в каких случая он запускается.

Всегда. Без него оценивать скорость просто бессмысленно, она будет маленькой. У Вас же мало кода, выполняемого только один раз, да он и почти не влияет на скорость. И уж точно цикл вычёркивания чуть ли не самый нагруженный, так что он уже после первого прохода весь будет в скомпилированном виде. Ну или поставьте запись времени непосредственно вокруг цикла вычёркивания и проверьте.