#include "nth_prime.h" #include #include #include #include #include namespace prime { using namespace std; number_t nth(size_t n, size_t sieve_size) { if (n < 1) throw domain_error("invalid prime number index"); if (n == 1) return 2; vector primes = {3, 5, 7, 11, 13, 17, 19, 23, 29, 31}; if (n - 2 < primes.size()) return primes[n - 2]; primes.reserve(n - 1); if (auto_sieve_size == sieve_size) { sieve_size = min(size_t(32 * 1024), size_t(0.5*(n * (log(n) + log(log(n)))) + 0.5)); } vector sieve(sieve_size, false); vector next(n - 1, 0); sieve[0] = true; number_t sieve_first = 1; auto process_prime = [&](size_t pn) { number_t i = next[pn]; for (; i < sieve_size; i += primes[pn]) { sieve[i] = true; } next[pn] = i - sieve_size; }; auto create_next = [&sieve_first](number_t p) -> number_t { return (p*p - sieve_first)/2; }; transform(primes.cbegin(), primes.cend(), next.begin(), create_next); for (auto p : primes) sieve[(p - sieve_first) / 2] = true; while (true) { for (size_t pn = 0; pn < primes.size(); ++pn) { process_prime(pn); } for (size_t i = 0; i < sieve_size; ++i) { if (sieve[i]) continue; number_t p = sieve_first + 2*i; size_t pn = primes.size(); primes.push_back(p); if (pn == n - 2) return p; next[pn] = create_next(p); sieve[i] = true; process_prime(pn); } number_t next_sieve_first = sieve_first + 2*sieve_size; if (next_sieve_first < sieve_first) throw runtime_error("failed to reach prime number"); sieve_first = next_sieve_first; memset(&sieve[0], false, sieve.size()*sizeof(sieve[0])); } } }