dujiao_sieve.hpp¶
Summatory Mobius and Euler-phi functions with a linear-sieve prefix and quotient-block memoization for large arguments.
Verified by sum_of_totient_function.
\[
\displaystyle M(x)=\sum_{i=1}^x \mu(i),\quad \Phi(x)=\sum_{i=1}^x \varphi(i)
\]
Implementation¶
#ifndef NOYA_DUJIAO_SIEVE_HPP
#define NOYA_DUJIAO_SIEVE_HPP 1
/// @complexity Time: O(L) preprocessing; quotient-block memoized queries use the standard O(n^(2/3)) regime.
/// Space: O(L + sqrt(n)) cached states.
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <unordered_map>
#include <vector>
namespace noya {
/// @brief Summatory Mobius and Euler-phi functions with a linear-sieve prefix
/// and quotient-block memoization for large arguments.
struct dujiao_sieve {
std::vector<int> primes;
std::vector<int> least_prime;
std::vector<std::int64_t> mobius_prefix;
std::vector<__int128> phi_prefix;
std::unordered_map<std::uint64_t, std::int64_t> mobius_cache;
std::unordered_map<std::uint64_t, __int128> phi_cache;
explicit dujiao_sieve(int limit) { prepare(limit); }
/// @brief Rebuild the direct prefix tables through limit.
void prepare(int limit) {
assert(limit >= 0);
primes.clear();
least_prime.assign(limit + 1, 0);
std::vector<std::int64_t> mobius(limit + 1);
std::vector<std::uint64_t> phi(limit + 1);
if (limit >= 1) {
mobius[1] = 1;
phi[1] = 1;
}
for (int value = 2; value <= limit; value++) {
if (least_prime[value] == 0) {
least_prime[value] = value;
primes.push_back(value);
mobius[value] = -1;
phi[value] = value - 1;
}
for (int prime : primes) {
if (prime > least_prime[value] ||
std::uint64_t(value) * prime > std::uint64_t(limit)) {
break;
}
int product = value * prime;
least_prime[product] = prime;
if (value % prime == 0) {
mobius[product] = 0;
phi[product] = phi[value] * prime;
} else {
mobius[product] = -mobius[value];
phi[product] = phi[value] * (prime - 1);
}
}
}
mobius_prefix.assign(limit + 1, 0);
phi_prefix.assign(limit + 1, 0);
for (int value = 1; value <= limit; value++) {
mobius_prefix[value] = mobius_prefix[value - 1] + mobius[value];
phi_prefix[value] = phi_prefix[value - 1] + phi[value];
}
mobius_cache.clear();
phi_cache.clear();
}
/// @brief Return sum_{i=1}^n mu(i).
std::int64_t summatory_mobius(std::uint64_t n) {
if (n < mobius_prefix.size()) {
return mobius_prefix[n];
}
if (auto iterator = mobius_cache.find(n); iterator != mobius_cache.end()) {
return iterator->second;
}
std::int64_t result = 1;
for (std::uint64_t left = 2, right; left <= n; left = right + 1) {
std::uint64_t quotient = n / left;
right = n / quotient;
result -= std::int64_t(right - left + 1) * summatory_mobius(quotient);
}
return mobius_cache.emplace(n, result).first->second;
}
/// @brief Return sum_{i=1}^n phi(i) as a signed 128-bit integer.
__int128 summatory_phi(std::uint64_t n) {
if (n < phi_prefix.size()) {
return phi_prefix[n];
}
if (auto iterator = phi_cache.find(n); iterator != phi_cache.end()) {
return iterator->second;
}
__int128 result = __int128(n) * (n + 1) / 2;
for (std::uint64_t left = 2, right; left <= n; left = right + 1) {
std::uint64_t quotient = n / left;
right = n / quotient;
result -= __int128(right - left + 1) * summatory_phi(quotient);
}
return phi_cache.emplace(n, result).first->second;
}
};
} // namespace noya
#endif // NOYA_DUJIAO_SIEVE_HPP
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <unordered_map>
#include <vector>
/// @complexity Time: O(L) preprocessing; quotient-block memoized queries use the standard O(n^(2/3)) regime.
/// Space: O(L + sqrt(n)) cached states.
namespace noya {
/// @brief Summatory Mobius and Euler-phi functions with a linear-sieve prefix
/// and quotient-block memoization for large arguments.
struct dujiao_sieve {
std::vector<int> primes;
std::vector<int> least_prime;
std::vector<std::int64_t> mobius_prefix;
std::vector<__int128> phi_prefix;
std::unordered_map<std::uint64_t, std::int64_t> mobius_cache;
std::unordered_map<std::uint64_t, __int128> phi_cache;
explicit dujiao_sieve(int limit) { prepare(limit); }
/// @brief Rebuild the direct prefix tables through limit.
void prepare(int limit) {
assert(limit >= 0);
primes.clear();
least_prime.assign(limit + 1, 0);
std::vector<std::int64_t> mobius(limit + 1);
std::vector<std::uint64_t> phi(limit + 1);
if (limit >= 1) {
mobius[1] = 1;
phi[1] = 1;
}
for (int value = 2; value <= limit; value++) {
if (least_prime[value] == 0) {
least_prime[value] = value;
primes.push_back(value);
mobius[value] = -1;
phi[value] = value - 1;
}
for (int prime : primes) {
if (prime > least_prime[value] ||
std::uint64_t(value) * prime > std::uint64_t(limit)) {
break;
}
int product = value * prime;
least_prime[product] = prime;
if (value % prime == 0) {
mobius[product] = 0;
phi[product] = phi[value] * prime;
} else {
mobius[product] = -mobius[value];
phi[product] = phi[value] * (prime - 1);
}
}
}
mobius_prefix.assign(limit + 1, 0);
phi_prefix.assign(limit + 1, 0);
for (int value = 1; value <= limit; value++) {
mobius_prefix[value] = mobius_prefix[value - 1] + mobius[value];
phi_prefix[value] = phi_prefix[value - 1] + phi[value];
}
mobius_cache.clear();
phi_cache.clear();
}
/// @brief Return sum_{i=1}^n mu(i).
std::int64_t summatory_mobius(std::uint64_t n) {
if (n < mobius_prefix.size()) {
return mobius_prefix[n];
}
if (auto iterator = mobius_cache.find(n); iterator != mobius_cache.end()) {
return iterator->second;
}
std::int64_t result = 1;
for (std::uint64_t left = 2, right; left <= n; left = right + 1) {
std::uint64_t quotient = n / left;
right = n / quotient;
result -= std::int64_t(right - left + 1) * summatory_mobius(quotient);
}
return mobius_cache.emplace(n, result).first->second;
}
/// @brief Return sum_{i=1}^n phi(i) as a signed 128-bit integer.
__int128 summatory_phi(std::uint64_t n) {
if (n < phi_prefix.size()) {
return phi_prefix[n];
}
if (auto iterator = phi_cache.find(n); iterator != phi_cache.end()) {
return iterator->second;
}
__int128 result = __int128(n) * (n + 1) / 2;
for (std::uint64_t left = 2, right; left <= n; left = right + 1) {
std::uint64_t quotient = n / left;
right = n / quotient;
result -= __int128(right - left + 1) * summatory_phi(quotient);
}
return phi_cache.emplace(n, result).first->second;
}
};
} // namespace noya