Squarefree Count
(math/squarefree_count.hpp)
- View this file on GitHub
- Last update: 2026-07-16 21:30:39+09:00
- Include:
#include "math/squarefree_count.hpp"
Overview
count_squarefree(n) returns the number of square-free positive integers not
greater than n. An integer is square-free when it is not divisible by the
square of any prime.
The starting identity is
\[Q(n)=\sum_{d\leq\sqrt n}\mu(d) \left\lfloor\frac{n}{d^2}\right\rfloor.\]Evaluating every term would require $O(\sqrt n)$ time. The implementation chooses a fifth-root split, sieves the Möbius function through $O(n^{2/5})$, and obtains the remaining grouped terms from selected values of the Mertens prefix sum.
Interface
uint64_t count_squarefree(uint64_t n);
| Function | Description | Complexity |
|---|---|---|
count_squarefree(n) |
Counts square-free integers in [1, n]. Returns 0 for n == 0. |
$O(n^{2/5})$ time and memory |
All roots are computed with exact integer arithmetic. The implementation
supports the full uint64_t input range for which its internal sieve limit fits
in int; in particular, it directly supports the Library Checker bound
$n\leq10^{18}$. The returned count never exceeds n.
Example
#include "math/squarefree_count.hpp"
#include <cassert>
int main() {
assert(m1une::math::count_squarefree(10) == 7);
assert(m1une::math::count_squarefree(100) == 61);
}
Depends on
Required by
Verified with
Code
#ifndef M1UNE_MATH_SQUAREFREE_COUNT_HPP
#define M1UNE_MATH_SQUAREFREE_COUNT_HPP 1
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <limits>
#include <vector>
#include "integer_arithmetic.hpp"
namespace m1une {
namespace math {
namespace squarefree_count_detail {
inline std::vector<int> mobius_prefix(int limit) {
std::vector<std::int8_t> mobius(limit + 1, 0);
if (limit >= 1) mobius[1] = 1;
{
std::vector<int> primes;
std::vector<bool> composite(limit + 1, false);
for (int value = 2; value <= limit; value++) {
if (!composite[value]) {
primes.push_back(value);
mobius[value] = -1;
}
for (int prime : primes) {
if (prime > limit / value) break;
const int product = prime * value;
composite[product] = true;
if (value % prime == 0) {
mobius[product] = 0;
break;
}
mobius[product] = static_cast<std::int8_t>(-mobius[value]);
}
}
}
std::vector<int> prefix(limit + 1, 0);
for (int value = 1; value <= limit; value++) {
prefix[value] = prefix[value - 1] + mobius[value];
}
return prefix;
}
} // namespace squarefree_count_detail
// Returns the number of square-free positive integers not greater than n.
inline uint64_t count_squarefree(uint64_t n) {
if (n == 0) return 0;
const uint64_t split = std::max<uint64_t>(1, floor_kth_root(n, 5U));
const uint64_t sieve_limit_u64 = isqrt(n / split);
assert(
sieve_limit_u64 <=
static_cast<uint64_t>(std::numeric_limits<int>::max())
);
const int sieve_limit = static_cast<int>(sieve_limit_u64);
const std::vector<int> mertens =
squarefree_count_detail::mobius_prefix(sieve_limit);
__int128_t direct_sum = 0;
for (int divisor = 1; divisor <= sieve_limit; divisor++) {
direct_sum += static_cast<__int128_t>(
n / static_cast<uint64_t>(divisor) /
static_cast<uint64_t>(divisor)
) * (mertens[divisor] - mertens[divisor - 1]);
}
std::vector<std::int64_t> large_mertens;
large_mertens.reserve(split - 1);
std::int64_t large_mertens_sum = 0;
for (uint64_t index = split; index-- > 1;) {
const uint64_t argument = isqrt(n / index);
const uint64_t square_root = isqrt(argument);
std::int64_t value = 1;
const uint64_t small_quotient_limit =
argument / (square_root + 1);
for (uint64_t quotient = 1;
quotient <= small_quotient_limit;
quotient++) {
const uint64_t multiplicity =
argument / quotient - argument / (quotient + 1);
value -= static_cast<std::int64_t>(multiplicity) *
mertens[static_cast<int>(quotient)];
}
for (uint64_t divisor = 2; divisor <= square_root; divisor++) {
const uint64_t quotient = argument / divisor;
if (quotient <= sieve_limit_u64) {
value -= mertens[static_cast<int>(quotient)];
} else {
const uint64_t previous_argument =
index * divisor * divisor;
assert(previous_argument < split);
const uint64_t position = split - previous_argument - 1;
assert(position < large_mertens.size());
value -= large_mertens[position];
}
}
large_mertens.push_back(value);
large_mertens_sum += value;
}
const std::int64_t grouped_sum =
large_mertens_sum -
static_cast<std::int64_t>(split - 1) * mertens[sieve_limit];
const __int128_t answer = direct_sum + grouped_sum;
assert(answer >= 0);
assert(answer <= std::numeric_limits<uint64_t>::max());
return static_cast<uint64_t>(answer);
}
} // namespace math
} // namespace m1une
#endif // M1UNE_MATH_SQUAREFREE_COUNT_HPP#line 1 "math/squarefree_count.hpp"
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <limits>
#include <vector>
#line 1 "math/integer_arithmetic.hpp"
#line 5 "math/integer_arithmetic.hpp"
#include <concepts>
#line 7 "math/integer_arithmetic.hpp"
#include <optional>
#include <type_traits>
namespace m1une {
namespace math {
namespace integer_arithmetic_detail {
template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
constexpr std::optional<T> checked_multiply(T first, T second) {
constexpr T minimum = std::numeric_limits<T>::min();
constexpr T maximum = std::numeric_limits<T>::max();
if constexpr (std::unsigned_integral<T>) {
if (second != 0 && maximum / second < first) return std::nullopt;
} else {
if (0 < first) {
if (0 < second) {
if (maximum / second < first) return std::nullopt;
} else if (second < minimum / first) {
return std::nullopt;
}
} else if (first < 0) {
if (0 < second) {
if (first < minimum / second) return std::nullopt;
} else if (second < maximum / first) {
return std::nullopt;
}
}
}
return T(first * second);
}
template <std::unsigned_integral T>
constexpr bool kth_power_leq(T base, unsigned exponent, T limit) {
assert(exponent > 0);
if (base <= 1) return base <= limit;
const T multiplication_limit = limit / base;
T product = 1;
for (unsigned i = 0; i < exponent; i++) {
if (product > multiplication_limit) return false;
product *= base;
}
return true;
}
} // namespace integer_arithmetic_detail
// Returns floor(sqrt(value)) exactly, without floating-point arithmetic.
template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
constexpr T isqrt(T value) {
if constexpr (std::signed_integral<T>) assert(0 <= value);
if (value <= 1) return value;
T low = 1;
T high = value / 2 + 1;
while (low < high) {
T middle = low + (high - low + 1) / 2;
if (middle <= value / middle) {
low = middle;
} else {
high = middle - 1;
}
}
return low;
}
template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
constexpr T floor_sqrt(T value) {
return isqrt(value);
}
// Returns ceil(sqrt(value)) exactly, without floating-point arithmetic.
template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
constexpr T ceil_sqrt(T value) {
T result = isqrt(value);
if (result == 0) return 0;
if (result != 0 && value / result == result && value % result == 0) {
return result;
}
return result + 1;
}
// Returns floor(value^(1 / degree)) exactly, without floating-point arithmetic.
template <std::integral T, std::integral Degree>
requires(
!std::same_as<std::remove_cv_t<T>, bool>
&& !std::same_as<std::remove_cv_t<Degree>, bool>
)
constexpr T floor_kth_root(T value, Degree degree) {
if constexpr (std::signed_integral<T>) {
assert(0 <= value);
if (value < 0) return T();
}
assert(0 < degree);
if (degree <= 0) return T();
if (value <= 1 || degree == 1) return value;
if (degree == 2) return isqrt(value);
using U = std::make_unsigned_t<T>;
using UDegree = std::make_unsigned_t<Degree>;
constexpr int digits = std::numeric_limits<U>::digits;
const UDegree unsigned_degree = static_cast<UDegree>(degree);
if (unsigned_degree >= static_cast<UDegree>(digits)) return T(1);
const unsigned exponent = static_cast<unsigned>(unsigned_degree);
const U unsigned_value = static_cast<U>(value);
int bit_width = 0;
for (U remaining = unsigned_value; remaining != 0; remaining >>= 1) {
bit_width++;
}
const int root_bits =
(bit_width + static_cast<int>(exponent) - 1) /
static_cast<int>(exponent);
U low = 1;
U high = U(1) << root_bits;
while (high - low > 1) {
const U middle = low + (high - low) / 2;
if (
integer_arithmetic_detail::kth_power_leq(
middle, exponent, unsigned_value
)
) {
low = middle;
} else {
high = middle;
}
}
return static_cast<T>(low);
}
// Returns base^exponent, or nullopt when the result does not fit in T.
template <std::integral T, std::unsigned_integral Exponent>
requires(
!std::same_as<std::remove_cv_t<T>, bool>
&& !std::same_as<std::remove_cv_t<Exponent>, bool>
)
constexpr std::optional<T> checked_ipow(T base, Exponent exponent) {
T result = 1;
while (exponent != 0) {
if (exponent & 1) {
auto product =
integer_arithmetic_detail::checked_multiply(result, base);
if (!product.has_value()) return std::nullopt;
result = *product;
}
exponent >>= 1;
if (exponent != 0) {
auto square =
integer_arithmetic_detail::checked_multiply(base, base);
if (!square.has_value()) return std::nullopt;
base = *square;
}
}
return result;
}
template <std::integral T, std::unsigned_integral Exponent>
requires(
!std::same_as<std::remove_cv_t<T>, bool>
&& !std::same_as<std::remove_cv_t<Exponent>, bool>
)
constexpr std::optional<T> checked_integer_pow(T base, Exponent exponent) {
return checked_ipow(base, exponent);
}
// Returns base^exponent. The result must be representable by T.
template <std::integral T, std::unsigned_integral Exponent>
requires(
!std::same_as<std::remove_cv_t<T>, bool>
&& !std::same_as<std::remove_cv_t<Exponent>, bool>
)
constexpr T ipow(T base, Exponent exponent) {
std::optional<T> result = checked_ipow(base, exponent);
assert(result.has_value());
return result.value_or(T());
}
template <std::integral T, std::unsigned_integral Exponent>
requires(
!std::same_as<std::remove_cv_t<T>, bool>
&& !std::same_as<std::remove_cv_t<Exponent>, bool>
)
constexpr T integer_pow(T base, Exponent exponent) {
return ipow(base, exponent);
}
} // namespace math
} // namespace m1une
#line 11 "math/squarefree_count.hpp"
namespace m1une {
namespace math {
namespace squarefree_count_detail {
inline std::vector<int> mobius_prefix(int limit) {
std::vector<std::int8_t> mobius(limit + 1, 0);
if (limit >= 1) mobius[1] = 1;
{
std::vector<int> primes;
std::vector<bool> composite(limit + 1, false);
for (int value = 2; value <= limit; value++) {
if (!composite[value]) {
primes.push_back(value);
mobius[value] = -1;
}
for (int prime : primes) {
if (prime > limit / value) break;
const int product = prime * value;
composite[product] = true;
if (value % prime == 0) {
mobius[product] = 0;
break;
}
mobius[product] = static_cast<std::int8_t>(-mobius[value]);
}
}
}
std::vector<int> prefix(limit + 1, 0);
for (int value = 1; value <= limit; value++) {
prefix[value] = prefix[value - 1] + mobius[value];
}
return prefix;
}
} // namespace squarefree_count_detail
// Returns the number of square-free positive integers not greater than n.
inline uint64_t count_squarefree(uint64_t n) {
if (n == 0) return 0;
const uint64_t split = std::max<uint64_t>(1, floor_kth_root(n, 5U));
const uint64_t sieve_limit_u64 = isqrt(n / split);
assert(
sieve_limit_u64 <=
static_cast<uint64_t>(std::numeric_limits<int>::max())
);
const int sieve_limit = static_cast<int>(sieve_limit_u64);
const std::vector<int> mertens =
squarefree_count_detail::mobius_prefix(sieve_limit);
__int128_t direct_sum = 0;
for (int divisor = 1; divisor <= sieve_limit; divisor++) {
direct_sum += static_cast<__int128_t>(
n / static_cast<uint64_t>(divisor) /
static_cast<uint64_t>(divisor)
) * (mertens[divisor] - mertens[divisor - 1]);
}
std::vector<std::int64_t> large_mertens;
large_mertens.reserve(split - 1);
std::int64_t large_mertens_sum = 0;
for (uint64_t index = split; index-- > 1;) {
const uint64_t argument = isqrt(n / index);
const uint64_t square_root = isqrt(argument);
std::int64_t value = 1;
const uint64_t small_quotient_limit =
argument / (square_root + 1);
for (uint64_t quotient = 1;
quotient <= small_quotient_limit;
quotient++) {
const uint64_t multiplicity =
argument / quotient - argument / (quotient + 1);
value -= static_cast<std::int64_t>(multiplicity) *
mertens[static_cast<int>(quotient)];
}
for (uint64_t divisor = 2; divisor <= square_root; divisor++) {
const uint64_t quotient = argument / divisor;
if (quotient <= sieve_limit_u64) {
value -= mertens[static_cast<int>(quotient)];
} else {
const uint64_t previous_argument =
index * divisor * divisor;
assert(previous_argument < split);
const uint64_t position = split - previous_argument - 1;
assert(position < large_mertens.size());
value -= large_mertens[position];
}
}
large_mertens.push_back(value);
large_mertens_sum += value;
}
const std::int64_t grouped_sum =
large_mertens_sum -
static_cast<std::int64_t>(split - 1) * mertens[sieve_limit];
const __int128_t answer = direct_sum + grouped_sum;
assert(answer >= 0);
assert(answer <= std::numeric_limits<uint64_t>::max());
return static_cast<uint64_t>(answer);
}
} // namespace math
} // namespace m1une