Instantly share code, notes, and snippets.

# math314/optimized_prime_sum.cpp Created Jun 4, 2016

What would you like to do?
 #include using namespace std; using ll = long long; const int MOD = int(1e9) + 7; //素数列挙 bool prime[1000001]; //10^6 vector 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; } } ll prime_sum(ll n) { int n2 = (int)sqrt(n); while (ll(n2 + 1) * (n2 + 1) <= n) n2++; const int k = lower_bound(prs.begin(), prs.end(), n2 + 1) - prs.begin() - 1; vector n_div_p; for (int i = 1; i <= n2; i++) { n_div_p.push_back(i); n_div_p.push_back(n / i); } sort(n_div_p.begin(), n_div_p.end()); n_div_p.erase(unique(n_div_p.begin(), n_div_p.end()), n_div_p.end()); vector ndp_table(n2 + 1); for (auto x : n_div_p) if (x <= n2) ndp_table[(int)x]++; for(int i = 0; i < n2; i++) ndp_table[i + 1] += ndp_table[i]; vector dp; for (auto x : n_div_p) { if (x % 2 == 0) dp.push_back((x / 2 % MOD)*((x + 1) % MOD) % MOD - 1); else dp.push_back((x % MOD)*((x + 1) / 2 % MOD) % MOD - 1); } for(int i = 0; i <= k; i++) { const ll p = prs[i]; for (int j = int(n_div_p.size()) - 1; j >= 0; j--) { const ll x = n_div_p[j]; if (x < p * p) break; const ll np = x / p; int n_div_p_idx; if (np <= n2) n_div_p_idx = ndp_table[(int)np] - 1; else n_div_p_idx = int(n_div_p.size()) - ndp_table[(int)(n / np)]; assert(n_div_p[n_div_p_idx] == np); (dp[j] -= p * (dp[n_div_p_idx] - prime_sum_small_part_memo[p - 1])) %= MOD; } } return (dp.back() % MOD + MOD) % MOD; }