Генерация достаточно больших простых чисел

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

• показывать N первых простых чисел;
• выводить простые числа в заданных пределах генерации [A;B];
• показывать количество простых чисел в заданном натуральном интервале значений [A;B];
• искать простое число с порядковым номером N.

Ниже показан шаблонный класс и его тест, консольная программа проверялась в актуальной сборке Visual Studio 2019.

#include <algorithm>
#include <iostream>
#include <cmath>
#include <cstdint>
#include <vector>
#include <limits>
using namespace std;
typedef uint64_t integer; /* тип данных для простых чисел */

template <typename integer> class prime_generator {
private:
 void find_primes (integer); //найти простые числа, начиная с аргумента
 integer count_ = 0; //счётчик чисел
 integer limit_; //предел счёта
 integer index_ = 0;
 integer increment_;
 vector <integer> primes_; //вектор простых чисел
 vector <bool> sieve_; //решето Эратосфена
 integer sieve_limit_ = 0; //лимит для решета
public:
 explicit prime_generator (integer initial_limit = 100, integer increment = 100000);
  //конструктор
 integer next_prime(); //Следующее простое число
 integer count() const { return count_; }
};

template<typename integer> integer next_odd_number(integer n) {
 //Следующее нечётное число
 return n % 2 == 0 ? n + 1 : n;
}

template<typename integer> //шаблон конструктора
 prime_generator<integer>::prime_generator (integer initial_limit, integer increment)
 : limit_(next_odd_number(initial_limit)), increment_(increment) {
 primes_.push_back(2);
 find_primes(3);
}

template<typename integer> //шаблон метода для поиска следующего простого числа
 integer prime_generator<integer>::next_prime() {
 if (index_ == primes_.size()) {
  if (numeric_limits<integer>::max() - increment_ < limit_) return 0;
  int start = limit_ + 2;
  limit_ = next_odd_number(limit_ + increment_);
  primes_.clear();
  find_primes(start);
 }
 ++count_;
 return primes_[index_++];
}

template<typename integer> integer isqrt(integer n) {
 //шаблон для поиска границы проверки делителей
 return next_odd_number(static_cast<integer>(sqrt(n)));
}

template<typename integer>
 void prime_generator<integer>::find_primes(integer start) {
 //поиск простых чисел начиная со start, применяем решето
 index_ = 0;
 integer new_limit = isqrt(limit_);
 sieve_.resize(new_limit / 2);
 for (integer p = 3; p * p <= new_limit; p += 2) {
  if (sieve_[p / 2 - 1]) continue;
  integer q = p * max(p, next_odd_number((sieve_limit_ + p - 1) / p));
  for (; q <= new_limit; q += 2 * p) sieve_[q / 2 - 1] = true;
 }
 sieve_limit_ = new_limit;
 size_t count = (limit_ - start) / 2 + 1;
 vector <bool> composite(count, false);
 for (integer p = 3; p <= new_limit; p += 2) {
  if (sieve_[p / 2 - 1]) continue;
  integer q = p * max(p, next_odd_number((start + p - 1) / p)) - start;
  q /= 2;
  for (; q < count; q += p) composite[q] = true;
 }
 for (integer p = 0; p < count; ++p) {
  if (!composite[p]) primes_.push_back(p * 2 + start);
 }
}

int main() {
 //Начальная инициализация:
 prime_generator <integer> pgen(100, 500000);
 //1. N первых простых чисел
 const integer n = 10;
 cout << "First " << n <<" primes: ";
 for (integer i = 0; i < n; i++) {
  integer p = pgen.next_prime();
  if (i != 0) cout << ", ";
  cout << p;
 }
 //2. простые числа между A и B, используем ранее сгенерированные данные
 integer A = 100, B = 300;
 integer p;
 cout << endl << "Primes between "<< A <<" and " << B << ": ";
 for (integer n = 0;; ) {
  p = pgen.next_prime();
  if (p > B) break;
  if (p >= A) {
   if (n != 0) cout << ", ";
   cout << p;
   n++;
  }
 }
 //3. количество простых чисел в интервале значений [A;B], используем ранее сгенерированные данные
 integer count = 0;
 A = 1000; B = 2000;
 for (;;) {
  p = pgen.next_prime();
  if (p > B) break;
  if (p >= A) count++;
 }
 cout << endl << "Number of primes between " << A << " and " << B << ": " << count << endl;
 //4. простое число с большим порядковым номером N, будет увеличивать решето и медленно
 integer prime;
 int k = 0;
 for (integer n = 1e4; n <= 1e8; n *= 10) { //от 10 тысяч до 100 миллионов по степеням 10
  while (pgen.count() != n) prime = pgen.next_prime();
  cout << n << "th prime (step "<< ++k <<" from 5): " << prime << endl;
 }
 return 0;
}

Конечно, достаточно большое "решето" можно подготовить один раз, а потом пользоваться им, например, сохранив данные в бинарном файле.

Результаты вывода программы таковы (последнее число искалось достаточно долго, так как размер решета динамически изменялся):

First 10 primes: 2, 3, 5, 7, 11, 13, 17, 19, 23, 29
Primes between 100 and 300: 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293
Number of primes between 1000 and 2000: 135
10000th prime (step 1 from 5): 104729
100000th prime (step 2 from 5): 1299709
1000000th prime (step 3 from 5): 15485863
10000000th prime (step 4 from 5): 179424673
100000000th prime (step 5 from 5): 2038074743

Скачать простые числа из диапазона от 2 до 10000000 одним файлом можно в этой моей статье. Номер строки развёрнутого файла, в которой находится простое число, будет соответствовать порядковому номеру этого числа в списке простых.

Проверить на простоту отдельное число можно здесь.

Ищем миллионное простое число коротким кодом без решета, но с помощью динамического массива.

#include <iostream>

int main() {
 auto* primes = new uint32_t[1000000]{ 2, 3, 5 };

 auto isPrime = [&](uint32_t n) {
  for (int i = 0; primes[i] * primes[i] <= n; ++i) {
   if (n % primes[i] == 0) {
    return false;
   }
  }
  return true;
 };

 std::size_t i = 3;
 uint32_t n = 5;
 while (i < 1000000) {
  while (!isPrime(++n));
  primes[i++] = n;
 }
 std::cout << std::endl << primes[999999] << std::endl;

 delete[] primes;
 return 0;
}

Как вариант того же без лямбда-функции:

#include <iostream>


bool isPrime(unsigned n, const unsigned* primes) {
 for (int i = 0; primes[i] * primes[i] <= n; ++i) {
  if (n % primes[i] == 0) {
   return false;
  }
 }
 return true;
};

int main() {
 unsigned* primes = new unsigned[1000000]{ 2, 3, 5 };

 unsigned i = 3;
 unsigned n = 5;
 while (i < 1000000) {
  while (!isPrime(++n, primes));
  primes[i++] = n;
 }
 std::cout << std::endl << primes[999999] << std::endl;

 delete[] primes;
 return 0;
}

16.10.2021, 17:52 [1468 просмотров]