Алгоритм нахождения N первых простых чисел

В процессе обучения программированию я столкнулся с задачей поиска 1M первых простых чисел.

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

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

В результате проведенной работы на поиск 1 миллиона простых чисел у меня уходит всего 0,262 секунды, что, конечно же, не предел, но уже впечатляет. Алгоритм реализован на C++.

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

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

10 000 — 2,349 с.
100 000 — 293,552 с.
1 000 000 — не дождался

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

// Listing 1<br/>
 <br/>
#include <iostream><br/>
#include <vector><br/>
#include <ctime><br/>
 <br/>
using namespace std;<br/>
 <br/>
int main()<br/>
{<br/>
    clock_t start = clock();<br/>
 <br/>
    const int AIM = 10000;<br/>
 <br/>
    vector<int> prime_nums;<br/>
    prime_nums.push_back(2);<br/>
 <br/>
    for (int num = 3; prime_nums.size() < AIM; num += 2)<br/>
    {<br/>
        bool isprime = true;<br/>
        for (int i = 3; i <= num / 2; i += 2)<br/>
        {<br/>
            if (num % i == 0)<br/>
            {<br/>
                isprime = false;<br/>
                break;<br/>
            }<br/>
        }<br/>
 <br/>
        if (isprime)<br/>
            prime_nums.push_back(num);<br/>
    }<br/>
 <br/>
    cout << prime_nums[AIM - 1];<br/>
 <br/>
    clock_t end = clock();<br/>
    cout << "\n----------------" << endl;<br/>
    cout << "Time: " << double(end - start) / CLOCKS_PER_SEC << " sec" << endl;<br/>
}


Стало лучше, но миллиона я так и не смог дождаться
10 000 — 0,589 с.
100 000 — 72,662 с.
1 000 000 — после 20 минут ждать перестал

Понятно, что продолжать оптимизировать этот алгоритм можно еще долго, чем я в общем-то и занимался. Можно не проверять числа заканчивающиеся на 5 и 0. Можно составить список самых распространенных делителей и проверять каждое число сначала на них. В какой-то момент мне показалось, что я смогу решить проблему, проверяя новые числа только на уже найденные простые. Значительных результатов это не принесло. Проблема крылась в том, что с каждым новым числом, количество проверок возрастало, да и операция деления — одна из самых тяжелых.

Отправился в поиск. Узнал про решето Эратосфена. Но использовать алгоритм в том виде, как он есть для моей задачи нельзя. Заранее ведь неизвестно*, сколько надо просеять чисел, чтобы найти N простых. Поэтому я решил модифицировать алгоритм под свои нужды следующим образом.

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

// Listing 2<br/>
 <br/>
#include <iostream><br/>
#include <ctime><br/>
#include <vector><br/>
 <br/>
using namespace std;<br/>
 <br/>
int main()<br/>
{<br/>
    clock_t start = clock();<br/>
 <br/>
    const int AIM = 1000000;<br/>
 <br/>
    int startSize = AIM; // стартовый размер массива натуральных чисел<br/>
    int addSize = AIM; // размер дополнительного массива натуральных чисел<br/>
 <br/>
    vector<bool> nums(startSize);<br/>
    vector<int> primeNums(AIM);<br/>
 <br/>
    int foundPrimes = 0;<br/>
 <br/>
    for (int i = 2; i < startSize; i++)<br/>
        nums[i] = true;<br/>
 <br/>
    bool addition = false;<br/>
    int adder = 0;<br/>
    while (true)<br/>
    {<br/>
        if (addition)<br/>
        {<br/>
            nums.resize(nums.size() + addSize, true);<br/>
 <br/>
            // исключим составные числа простыми из предыдущих итераций<br/>
            for (int i = 0; i < foundPrimes; i++)<br/>
            {<br/>
                int cur_num = primeNums[i];<br/>
                if ((addSize + ((nums.size() - addSize) % cur_num)) < cur_num)<br/>
                    continue;<br/>
                for (int j = ((nums.size() - addSize) / cur_num) * cur_num; j < nums.size(); j += cur_num)<br/>
                    nums[j] = false;<br/>
            }<br/>
        }<br/>
        else<br/>
            addition = true;<br/>
 <br/>
 <br/>
        int iter;<br/>
        if (foundPrimes == 0)<br/>
            iter = 2;<br/>
        else<br/>
            iter = primeNums[foundPrimes - 1] + 2;<br/>
 <br/>
        for ( ; iter < nums.size(); iter++)<br/>
        {<br/>
            // выбираем очередное простое число<br/>
            if (nums[iter])<br/>
            {<br/>
                primeNums[foundPrimes] = iter;<br/>
                foundPrimes++;<br/>
                if (foundPrimes == AIM)<br/>
                    break;<br/>
                // просеиваем<br/>
                for (int j = iter + iter; j < nums.size(); j += iter)<br/>
                    nums[j] = false;<br/>
            }<br/>
            else<br/>
                continue;<br/>
 <br/>
        }<br/>
        if (foundPrimes == AIM)<br/>
            break;<br/>
    }<br/>
 <br/>
    cout << "Last prime: " << primeNums[AIM -1];<br/>
 <br/>
    clock_t end = clock();<br/>
    cout << "\n----------------" << endl;<br/>
    cout << "Time: " << double(end - start) / CLOCKS_PER_SEC << " sec" << endl;<br/>
    return 0;<br/>
}


Результаты даже слегка удивили:
10 000 — 0,001 с.
100 000 — 0,021 с.
1 000 000 — 0,262 с.

10 000 000 — 2,999 с.

На этом с задачей поиска 1M первых простых чисел я закончил. Результат меня вполне удовлетворил. Но что можно сделать ещё? Можно реализовать какой-нибудь из известных алгоритмов на основе решета Эратосфена, например, решето Сундарама. Но вряд ли результаты изменяться значительно. Возможно, решето Аткина сможет показать себя лучше — именно этот алгоритм чаще всего используется в современной криптографии.

* — неизвестно мне. halyavin рассказал в первом комментарии, что оказывается давно известно.

UPD
Относительно потребления памяти, описанные алгоритмы ведут себя следующим образом:

1 млн. простых чисел, максимальное ОЗУ
— перебор делителей: ~4Мб
— модифицированное решето Эратосфена: ~6Мб

10 млн. простых чисел, максимальное ОЗУ
— перебор делителей: ~39Мб
— модифицированное решето Эратосфена: ~61Мб
Share post

Similar posts

AdBlock has stolen the banner, but banners are not teeth — they will be back

More
Ads

Comments 29

    +4
    >Заранее ведь неизвестно, сколько надо просеять чисел, чтобы найти N простых.

    Довольно хорошо известно. Можно воспользоваться оценкой в ru.wikipedia.org/wiki/%D0%93%D0%B8%D0%BF%D0%BE%D1%82%D0%B5%D0%B7%D0%B0_%D0%A0%D0%B8%D0%BC%D0%B0%D0%BD%D0%B0, если вы верите в гипотезу Римана. Если не верите, то есть худшие оценки — mathworld.wolfram.com/PrimeNumberTheorem.html.
      0
      Спасибо за ссылки, буду изучать.
        0
        Вот еще одна — primes.utm.edu/howmany.shtml#better. Там есть недавние результаты с конкретными оценками на n-тое простое число.
        0
        На практике удобно использовать неравенство (не асимптотику, верную лишь в пределе!)
        Pn < n (ln n + ln ln n − ½) при n > 20, где Pnn-е простое число.
      • UFO just landed and posted this here
          +6
          Для этого надо все-таки учить алгебру в объеме, превышающем школьную программу.
            +1
            Что? Пусть x = p*q — проверяемое число, где p — простое, а q — возможно составное. Хотя бы одно из этих чисел меньше либо равно \sqrt{x}. Это число мы и найдём в качестве делителя.
              –4
              Простите, но это не доказательство.
                0
                Где именно нестрого?
                  –2
                  Мне кажется тут:
                  Хотя бы одно из этих чисел меньше либо равно \sqrt{x}.
                    +2
                    Если x = p*q, p > \sqrt{x}, q > \sqrt{x}, то умножив почленно получим x < x.
                      0
                      Поддерживаю. Подзабыл доказательство. Каюсь :)
                0
                В любом случае эти факты учатся в ВУЗовской программе.
                  –1
                  Я вот уже получаю второе высшее, но до сих пор ни один преподаватель по математике об этом не говорил.
              0
              Вы абсолютно правы. В окончательном варианте, который я не стал здесь приводить, я сделал именно так. И все же конечный результат по времени — неудовлетворительный.
                0
                Также обычно рассматривают числа виде 6k, 6k+1, 6k+2,6k+3,6k+4,6k+5. Здесь 6=2*3 произведению первых двух простых чисел.
                Из них 6k, 6k+2,6k+3,6k+4 делятся на 2 или 3.
                Поэтому надо проверять на простоту только числа вида 6k+1, 6k+5 (33.3% от всех).

                Это можно обобщить например рассматривая числа вида 30k, 30k+1,…
                Но это становится громоздким при программирование и дает малый выигрыш
                30k+1, 30k+7, 30k+11, 30k+13, 30k+17, 30k+19,30k+23, 30k+29 (26.6% от всех)
                0
                Результаты изменятся.

                Сложность простого решета Эратосфена O(N*log(log(N))).
                Сложность решета Аткина уже O(N/log(log(N)))

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

                Хотя, конечно, нужно еще оценить сложность вашей модификации.
                  +2
                  
                  if (addition)
                  {
                  ..............
                  }
                  else
                       addition = true;
                  

                  Мне кажется, что else не нужен
                    0
                    Когда решил задачу нахождения простых чисел (правда до N) то решил использовать веростностный тест Ферма + алгоритм грубой силы. Вместо алгоритма грубой силы можно использовать другие, но подбор другого алгоритма был не столь важен, сколько реализация вероятстного теста Ферма.

                    Тесты на простоту числа.
                      +2
                      Не вижу в вашем посте именно алгоритма. Вижу его реализацию на одном конкретном языке. А алгоритма не вижу.

                      «Программа хорошо документирована на языке C» — это для башорга, а не статей от которых ожидается хоть слегка научная правильность.
                        –1
                        Поясните, пожалуйста, в чем разница записи алгоритма на языке C++ и на любом другом формальном языке (допустим математическом).
                          0
                          Математический язык имеет строго определённую и простую семантику. Язык C++ имеет очень сложную и местами неопределённую (undefined behavior) семантику.
                            0
                            Спасибо за пояснение, учту на будущее.
                        +1
                        Описание алгоритмов будет неполным без указанием размера требуемой памяти. Ведь компромисс быстродействие/используемая память встречается практически в каждой задаче.
                          0
                          Добавил информацию
                          0
                          Решето Сундарама реализовывать нет смысла — оно медленнее Эратосфена.
                            0
                            удалено
                              0
                              Кажется листинги с кодом сломались.
                                0
                                Листинги смотрятся коряво, пришлось вручную перенести в html-файл.
                                Станно, что мой haswell 4590 просчитал 100тыс за 9.7сек, а 1млн за 123.5 сек, что примерно в 400 раз хуже, чем у автора. Что у него за компьютер?

                                Only users with full accounts can post comments. Log in, please.