squarefree_count.hpp¶
Count square-free positive integers at most n. Mobius inversion handles divisors up to D = sqrt(n / n^(1/5)) directly. The remaining terms share only O(n^(1/5)) distinct square-root quotients; their Mertens values are recovered in descending order with the quotient-block recurrence, reducing both time and memory from O(sqrt(n)) to O(n^(2/5)).
Verified by counting_squarefrees.
\[
\displaystyle #\{n\le N:\mu(n)^2=1\}
\]
Implementation¶
#ifndef NOYA_SQUAREFREE_COUNT_HPP
#define NOYA_SQUAREFREE_COUNT_HPP 1
/// @complexity Time: O(n^(2/5) log log n).
/// Space: O(n^(2/5)).
#include <cassert>
#include <cmath>
#include <cstdint>
#include <algorithm>
#include <limits>
#include <vector>
namespace noya {
namespace squarefree_count_detail {
inline std::uint64_t floor_square_root(std::uint64_t value) {
std::uint64_t root = std::sqrt(static_cast<long double>(value));
while ((root + 1) <= value / (root + 1)) {
root++;
}
while (root > 0 && root > value / root) {
root--;
}
return root;
}
inline std::uint64_t floor_fifth_root(std::uint64_t value) {
auto at_most = [&](std::uint64_t root) {
unsigned __int128 product = 1;
for (int exponent = 0; exponent < 5; exponent++) {
product *= root;
if (product > value) {
return false;
}
}
return true;
};
std::uint64_t low = 0;
std::uint64_t high = 65536;
while (high - low > 1) {
std::uint64_t middle = (low + high) / 2;
(at_most(middle) ? low : high) = middle;
}
return low;
}
inline std::vector<std::int8_t> mobius_table(int limit) {
std::vector<int> least_prime(limit + 1);
std::vector<int> primes;
std::vector<std::int8_t> mobius(limit + 1);
if (limit >= 1) {
mobius[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;
}
for (int prime : primes) {
if (prime > least_prime[value] || value > limit / prime) {
break;
}
least_prime[value * prime] = prime;
mobius[value * prime] =
value % prime == 0 ? 0 : std::int8_t(-mobius[value]);
}
}
return mobius;
}
} // namespace squarefree_count_detail
/// @brief Count square-free positive integers at most n. Mobius inversion
/// handles divisors up to `D = sqrt(n / n^(1/5))` directly. The remaining
/// terms share only `O(n^(1/5))` distinct square-root quotients; their Mertens
/// values are recovered in descending order with the quotient-block
/// recurrence, reducing both time and memory from O(sqrt(n)) to O(n^(2/5)).
inline std::uint64_t squarefree_count(std::uint64_t n) {
if (n == 0) {
return 0;
}
std::uint64_t split =
std::max<std::uint64_t>(1, squarefree_count_detail::floor_fifth_root(n));
std::uint64_t direct_limit =
squarefree_count_detail::floor_square_root(n / split);
assert(direct_limit <=
std::uint64_t(std::numeric_limits<int>::max()));
auto mobius =
squarefree_count_detail::mobius_table(int(direct_limit));
std::vector<std::int64_t> mertens(direct_limit + 1);
for (std::uint64_t value = 1; value <= direct_limit; value++) {
mertens[value] = mertens[value - 1] + mobius[value];
}
__int128 direct_sum = 0;
for (std::uint64_t value = 1; value <= direct_limit; value++) {
direct_sum += __int128(mobius[value]) * (n / value / value);
}
std::vector<std::int64_t> large_mertens;
large_mertens.reserve(split > 0 ? split - 1 : 0);
__int128 large_sum = 0;
for (std::uint64_t index = split; index-- > 1;) {
std::uint64_t argument =
squarefree_count_detail::floor_square_root(n / index);
std::uint64_t small_root =
squarefree_count_detail::floor_square_root(argument);
std::int64_t value = 1;
std::uint64_t quotient_boundary = argument / (small_root + 1);
for (std::uint64_t divisor = 1; divisor <= quotient_boundary; divisor++) {
value -= std::int64_t(argument / divisor - argument / (divisor + 1)) *
mertens[divisor];
}
for (std::uint64_t divisor = 2; divisor <= small_root; divisor++) {
std::uint64_t quotient = argument / divisor;
if (quotient <= direct_limit) {
value -= mertens[quotient];
} else {
std::uint64_t represented_index = index * divisor * divisor;
assert(represented_index < split);
value -= large_mertens[split - represented_index - 1];
}
}
large_mertens.push_back(value);
large_sum += value;
}
__int128 result = direct_sum + large_sum -
__int128(split - 1) * mertens[direct_limit];
assert(result >= 0);
return std::uint64_t(result);
}
} // namespace noya
#endif // NOYA_SQUAREFREE_COUNT_HPP
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <limits>
#include <vector>
/// @complexity Time: O(n^(2/5) log log n).
/// Space: O(n^(2/5)).
namespace noya {
namespace squarefree_count_detail {
inline std::uint64_t floor_square_root(std::uint64_t value) {
std::uint64_t root = std::sqrt(static_cast<long double>(value));
while ((root + 1) <= value / (root + 1)) {
root++;
}
while (root > 0 && root > value / root) {
root--;
}
return root;
}
inline std::uint64_t floor_fifth_root(std::uint64_t value) {
auto at_most = [&](std::uint64_t root) {
unsigned __int128 product = 1;
for (int exponent = 0; exponent < 5; exponent++) {
product *= root;
if (product > value) {
return false;
}
}
return true;
};
std::uint64_t low = 0;
std::uint64_t high = 65536;
while (high - low > 1) {
std::uint64_t middle = (low + high) / 2;
(at_most(middle) ? low : high) = middle;
}
return low;
}
inline std::vector<std::int8_t> mobius_table(int limit) {
std::vector<int> least_prime(limit + 1);
std::vector<int> primes;
std::vector<std::int8_t> mobius(limit + 1);
if (limit >= 1) {
mobius[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;
}
for (int prime : primes) {
if (prime > least_prime[value] || value > limit / prime) {
break;
}
least_prime[value * prime] = prime;
mobius[value * prime] =
value % prime == 0 ? 0 : std::int8_t(-mobius[value]);
}
}
return mobius;
}
} // namespace squarefree_count_detail
/// @brief Count square-free positive integers at most n. Mobius inversion
/// handles divisors up to `D = sqrt(n / n^(1/5))` directly. The remaining
/// terms share only `O(n^(1/5))` distinct square-root quotients; their Mertens
/// values are recovered in descending order with the quotient-block
/// recurrence, reducing both time and memory from O(sqrt(n)) to O(n^(2/5)).
inline std::uint64_t squarefree_count(std::uint64_t n) {
if (n == 0) {
return 0;
}
std::uint64_t split =
std::max<std::uint64_t>(1, squarefree_count_detail::floor_fifth_root(n));
std::uint64_t direct_limit =
squarefree_count_detail::floor_square_root(n / split);
assert(direct_limit <=
std::uint64_t(std::numeric_limits<int>::max()));
auto mobius =
squarefree_count_detail::mobius_table(int(direct_limit));
std::vector<std::int64_t> mertens(direct_limit + 1);
for (std::uint64_t value = 1; value <= direct_limit; value++) {
mertens[value] = mertens[value - 1] + mobius[value];
}
__int128 direct_sum = 0;
for (std::uint64_t value = 1; value <= direct_limit; value++) {
direct_sum += __int128(mobius[value]) * (n / value / value);
}
std::vector<std::int64_t> large_mertens;
large_mertens.reserve(split > 0 ? split - 1 : 0);
__int128 large_sum = 0;
for (std::uint64_t index = split; index-- > 1;) {
std::uint64_t argument =
squarefree_count_detail::floor_square_root(n / index);
std::uint64_t small_root =
squarefree_count_detail::floor_square_root(argument);
std::int64_t value = 1;
std::uint64_t quotient_boundary = argument / (small_root + 1);
for (std::uint64_t divisor = 1; divisor <= quotient_boundary; divisor++) {
value -= std::int64_t(argument / divisor - argument / (divisor + 1)) *
mertens[divisor];
}
for (std::uint64_t divisor = 2; divisor <= small_root; divisor++) {
std::uint64_t quotient = argument / divisor;
if (quotient <= direct_limit) {
value -= mertens[quotient];
} else {
std::uint64_t represented_index = index * divisor * divisor;
assert(represented_index < split);
value -= large_mertens[split - represented_index - 1];
}
}
large_mertens.push_back(value);
large_sum += value;
}
__int128 result = direct_sum + large_sum -
__int128(split - 1) * mertens[direct_limit];
assert(result >= 0);
return std::uint64_t(result);
}
} // namespace noya