prime_sum.cpp

#include <bits/stdc++.h>

using namespace std;
using ll = long long;

//素数列挙
bool prime[1000001]; //10^6
vector<int> prs;
void init_prime() {
	memset(prime, 1, sizeof(prime));
	prime[0] = prime[1] = false;
	for (int i = 2; i < sizeof(prime); i++) if (prime[i])
		for (int j = i * 2; j < sizeof(prime); j += i) prime[j] = false;
	for (int i = 2; i < sizeof(prime); i++) if (prime[i]) prs.push_back(i);
}

ll prime_sum_small_part_memo[sizeof(prime)];
void init_prime_sum_small_part() {
	for (int i = 0; i < sizeof(prime) - 1; i++) {
		prime_sum_small_part_memo[i + 1] = prime_sum_small_part_memo[i];
		if (prime[i + 1]) prime_sum_small_part_memo[i + 1] += i + 1;
	}
}

unordered_map<ll, ll> prime_sum_memo[80000];
//2 <= x <= prs[k]までの値で篩をした後に、残っている数のsum
ll prime_sum(ll n, int k) {
	if (k == -1) {
		return n * (n + 1) / 2 - 1; // 篩をしていないので、2~nまでの総和
	}
	ll p = prs[k];
	//枝刈り
	if (n < sizeof(prime) && n < p * p) {
		//愚直に計算した結果を返せば良い
		return prime_sum_small_part_memo[(int)n];
	}
	//メモ化
	auto it = prime_sum_memo[k].find(n);
	if (it != prime_sum_memo[k].end()) return it->second;
	auto& ret = prime_sum_memo[k][n];
	// n < p^2の時、pが最小の素因数になるような 2 <= x <= n なる x は存在しない
	//    -> prime_sum(n, k - 1) と一致する
	if (n < p * p) {
		return ret = prime_sum(n, k - 1);
	}
	// p^2 <= n の場合を考える
	// (2 <= x <= pまでの値で篩をした後に、残っているn以下の数のsum)
	//   = (2 <= x < pまでの値で篩をした後に、残っているn以下の数のsum) - { (pが最小の素因数となるn以下の数) - p }
	// { (pが最小の素因数となるn以下の数) - p } = p * (pが最小の素因数となるn/p以下の数)
	//   = p * { (2 <= x <= pまでの値で篩をした後に、残っているn/p以下の数のsum) - (p未満の素数のsum) }
	return ret = prime_sum(n, k - 1) - p * (prime_sum(n / p, k - 1) - prime_sum_small_part_memo[p - 1]);
}


ll prime_sum(ll n) {
	int n2 = (int)sqrt(n);
	while (ll(n2 + 1) * (n2 + 1) <= n) n2++;
	return prime_sum(n, n2);
}

int main() {
	init_prime();
	init_prime_sum_small_part();

	cout << prime_sum(20) << endl;
	cout << prime_sum(1e2) << endl;
	cout << prime_sum(1e5) << endl;
	cout << prime_sum(1e6) << endl;
	// cout << prime_sum(1e10) << endl;

	return 0;
}