divisor_summatory.hpp¶
Return sum_{i=1}^n tau(i) in O(sqrt(n)) quotient blocks.
Verified by enumerate_quotients.
\[
\displaystyle S(n)=\sum_{i=1}^n \tau(i)
\]
Implementation¶
#ifndef NOYA_DIVISOR_SUMMATORY_HPP
#define NOYA_DIVISOR_SUMMATORY_HPP 1
/// @complexity Time: O(sqrt(n)).
/// Space: O(1) for summatory functions; O(sqrt(n)) for quotient enumeration.
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <vector>
namespace noya {
/// @brief Return sum_{i=1}^n tau(i) in O(sqrt(n)) quotient blocks.
inline unsigned __int128 summatory_divisor_count(std::uint64_t n) {
unsigned __int128 result = 0;
for (std::uint64_t left = 1, right; left <= n; left = right + 1) {
std::uint64_t quotient = n / left;
right = n / quotient;
result += static_cast<unsigned __int128>(right - left + 1) * quotient;
}
return result;
}
/// @brief Return sum_{i=1}^n sigma(i) in O(sqrt(n)); n must be at most 1e18 so
/// the exact result fits unsigned 128-bit arithmetic.
inline unsigned __int128 summatory_divisor_sum(std::uint64_t n) {
assert(n <= 1000000000000000000ULL);
unsigned __int128 result = 0;
for (std::uint64_t left = 1, right; left <= n; left = right + 1) {
std::uint64_t quotient = n / left;
right = n / quotient;
unsigned __int128 count = right - left + 1;
unsigned __int128 endpoints = static_cast<unsigned __int128>(left) + right;
unsigned __int128 interval_sum =
(count % 2 == 0) ? count / 2 * endpoints : count * (endpoints / 2);
result += interval_sum * quotient;
}
return result;
}
/// @brief Return all distinct positive values floor(n / x) in increasing
/// order. Returns an empty vector when n is zero.
inline std::vector<std::uint64_t> distinct_quotients(std::uint64_t n) {
std::vector<std::uint64_t> result;
for (std::uint64_t left = 1, right; left <= n; left = right + 1) {
std::uint64_t quotient = n / left;
result.push_back(quotient);
right = n / quotient;
}
std::reverse(result.begin(), result.end());
return result;
}
} // namespace noya
#endif // NOYA_DIVISOR_SUMMATORY_HPP
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <vector>
/// @complexity Time: O(sqrt(n)).
/// Space: O(1) for summatory functions; O(sqrt(n)) for quotient enumeration.
namespace noya {
/// @brief Return sum_{i=1}^n tau(i) in O(sqrt(n)) quotient blocks.
inline unsigned __int128 summatory_divisor_count(std::uint64_t n) {
unsigned __int128 result = 0;
for (std::uint64_t left = 1, right; left <= n; left = right + 1) {
std::uint64_t quotient = n / left;
right = n / quotient;
result += static_cast<unsigned __int128>(right - left + 1) * quotient;
}
return result;
}
/// @brief Return sum_{i=1}^n sigma(i) in O(sqrt(n)); n must be at most 1e18 so
/// the exact result fits unsigned 128-bit arithmetic.
inline unsigned __int128 summatory_divisor_sum(std::uint64_t n) {
assert(n <= 1000000000000000000ULL);
unsigned __int128 result = 0;
for (std::uint64_t left = 1, right; left <= n; left = right + 1) {
std::uint64_t quotient = n / left;
right = n / quotient;
unsigned __int128 count = right - left + 1;
unsigned __int128 endpoints = static_cast<unsigned __int128>(left) + right;
unsigned __int128 interval_sum =
(count % 2 == 0) ? count / 2 * endpoints : count * (endpoints / 2);
result += interval_sum * quotient;
}
return result;
}
/// @brief Return all distinct positive values floor(n / x) in increasing
/// order. Returns an empty vector when n is zero.
inline std::vector<std::uint64_t> distinct_quotients(std::uint64_t n) {
std::vector<std::uint64_t> result;
for (std::uint64_t left = 1, right; left <= n; left = right + 1) {
std::uint64_t quotient = n / left;
result.push_back(quotient);
right = n / quotient;
}
std::reverse(result.begin(), result.end());
return result;
}
} // namespace noya