pekempeyのブログ

競技プログラミングに関する話題を書いていきます。

prime counting function

以下の問題を解く際に「n 以下の素数の個数」を高速に求める必要が出てくる。

http://codeforces.com/contest/665/problem/F

素数列を\(p_0,p_1,\ldots\)とする。

\(\phi(x,p_i)\)を\(x\)以下の正整数のうち\(p_0,p_1,\ldots,p_i\)を約数に持たないものの個数と定義する。このとき以下の漸化式が成り立つ。

\[ \phi(x,p_i)=\phi(x,p_{i-1})-\phi\left(\frac{x}{p_i},p_{i-1}\right) \]

さらに\(\pi(x,p_i)=\phi(x,p_i)+i+1\)と定義する。\(x \le p_i^2\)のとき\(\pi(x,p_i)=\pi(x)\)になることは試し割りの素数判定の原理を考えれば容易にわかる。2,3,5,7,...で割り切れない正整数の個数を数えているのであまり良く考えないと 1 を素数に含めてしまうので注意

これを最初の漸化式に当てはめると以下の漸化式が導ける。

\[ \pi(x,p_i)=\pi(x,p_{i-1})-\pi\left(\frac{x}{p_i},p_{i-1}\right)+i+1 \]

\( \pi(x,p_i) \) のテーブルを更新するとき、\( \pi(1..p_i^2-1,p_i) \) の値は変わらないことを使うと、DPテーブルの更新は大幅に高速化できる。

この手法を使うとpi(N/i)側が求められるという利点がある。1000000以下の素数で計算しても1sec程度である。

#include <iostream>
#include <ctime>
#include <vector>
#include <algorithm>

using namespace std;

const long long N = 5e5;
bool tab[N + 1];
long long pi[N + 1];
long long hi[N + 1];
vector<long long> primes;

int main() {
	long long n;
	cin >> n;
	// 1 is prime !!
	for (int i = 1; i <= N; i++) {
		pi[i] = i;
		hi[i] = n / i;
		tab[i] = true;
	}
	for (int i = 2; i <= N; i++) {
		if (tab[i]) {
			primes.push_back(i);
			for (int j = i * 2; j <= N; j += i) {
				tab[j] = false;
			}
		}
	}
	for (int i = 0; i < primes.size(); i++) {
		long long p = primes[i];
		// n/j>=p*p
		// j<=n/(p*p)
		for (int j = 1; j <= min(N, n / (p * p)); j++) {
			long long k = n / j / p;
			hi[j] -= (k <= N ? pi[k] : hi[j * p]) - (i + 1);
		}
		for (int j = N; j >= p * p; j--) {
			pi[j] -= pi[j / p] - (i + 1);
		}
	}
	long long ans = 0;
	for (long long p : primes) {
		if (p * p * p <= n) {
			ans++;
		}
		ans += max(0LL, hi[p] - pi[p]);
	}
	cout << ans << endl;
}