乗法的関数の和の内部処理

components/multiplicative_sum_backend.hpp · GitHub で見る

共通ページ

乗法的関数の prefix sum が、大きい上限での素数和と乗法的関数の和に使う補助部品です。通常は直接使う必要はありません。

平方因子を持たない数に対応する和を作り、強力数上の補正を加えます。素数冪の補正値は共有し、同じ値の再計算を避けます。3 の逆元を使います。

素数和は、完全乗法的な関数の和から素数ごとの Dirichlet 逆変換で復元します。この復元は \(\displaystyle O(\sqrt n\log\log n)\) 時間です。

作業領域は \(\displaystyle O(\sqrt n)\) です。

依存ライブラリ

整数商の列挙・添字

検証コード

コード

#line 1 "components/multiplicative_sum_backend.hpp"

#line 1 "components/floor_quotients.hpp"

#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <iterator>
#include <utility>
#include <vector>

namespace floor_quotient_detail {
inline std::uint64_t square_root(std::uint64_t n) {
	std::uint64_t root = std::sqrt(static_cast<long double>(n));
	while (__uint128_t(root) * root > n) --root;
	while (__uint128_t(root + 1) * (root + 1) <= n) ++root;
	return root;
}
inline std::uint64_t divide(std::uint64_t n, std::uint64_t divisor) {
	if constexpr (std::numeric_limits<double>::is_iec559 && std::numeric_limits<double>::radix == 2 && std::numeric_limits<double>::digits >= 53) {
		if (n <= (std::uint64_t(1) << 53)) {
			std::uint64_t value = double(n) / double(divisor);
			return value - (value * divisor > n);
		}
	}
	return n / divisor;
}
}

class FloorQuotientRange {
	std::uint64_t n_, root_, large_;

public:
	explicit FloorQuotientRange(std::uint64_t n = 0) : n_(n), root_(floor_quotient_detail::square_root(n)), large_(root_ - (n < root_ * (root_ + 1))) {}
	std::uint64_t max_n() const { return n_; }
	std::uint64_t root() const { return root_; }
	std::uint64_t large_size() const { return large_; }
	std::uint64_t size() const { return root_ + large_; }
	bool empty() const { return !n_; }
	std::uint64_t operator[](std::uint64_t i) const {
		assert(i < size());
		return i < root_ ? i + 1 : floor_quotient_detail::divide(n_, size() - i);
	}
	class iterator {
		const FloorQuotientRange *range_ = nullptr;
		std::uint64_t index_ = 0;
		friend class FloorQuotientRange;
		iterator(const FloorQuotientRange *range, std::uint64_t i) : range_(range), index_(i) {}

public:
		using value_type = std::uint64_t;
		using difference_type = std::int64_t;
		using iterator_category = std::input_iterator_tag;
		using iterator_concept = std::forward_iterator_tag;
		using reference = value_type;
		using pointer = void;
		iterator() = default;
		value_type operator*() const { return (*range_)[index_]; }
		iterator &operator++() { ++index_; return *this; }
		iterator operator++(int) { auto result = *this; ++*this; return result; }
		bool operator==(const iterator &) const = default;
	};
	iterator begin() const { return iterator(this, 0); }
	iterator end() const { return iterator(this, size()); }
};

class FloorQuotients {
	std::uint64_t n_, root_ = 0;
	std::size_t large_ = 0;
	std::vector<std::uint64_t> values_;

public:
	explicit FloorQuotients(std::uint64_t n = 0) : n_(n) {
		if (!n) return;
		root_ = floor_quotient_detail::square_root(n);
		large_ = n / (root_ + 1);
		values_.reserve(large_ + root_);
		for (std::size_t i = 1; i <= large_; ++i) values_.push_back(floor_quotient_detail::divide(n, i));
		for (auto x = root_; x; --x) values_.push_back(x);
	}
	FloorQuotients(const FloorQuotients&) = default;
	FloorQuotients& operator=(const FloorQuotients&) = default;
	FloorQuotients(FloorQuotients&& other) noexcept : FloorQuotients() { swap(other); }
	FloorQuotients& operator=(FloorQuotients&& other) noexcept {
		if (this != &other) { FloorQuotients empty; swap(empty); swap(other); }
		return *this;
	}
	void swap(FloorQuotients& other) noexcept {
		std::swap(n_, other.n_), std::swap(root_, other.root_), std::swap(large_, other.large_);
		values_.swap(other.values_);
	}
	std::uint64_t max_n() const { return n_; }
	std::uint64_t root() const { return root_; }
	std::size_t size() const { return values_.size(); }
	std::size_t large_size() const { return large_; }
	const std::vector<std::uint64_t>& values() const { return values_; }
	bool contains(std::uint64_t x) const { return x && x <= n_ && (x <= root_ || n_ / (n_ / x) == x); }
	std::size_t index(std::uint64_t x) const {
		assert(x && x <= n_);
		return x <= root_ ? size() - x : n_ / x - 1;
	}
	std::size_t divided_index(std::size_t i, std::uint64_t divisor) const {
		assert(i < size() && divisor && divisor <= values_[i]);
		if (i >= large_) return size() - values_[i] / divisor;
		const auto d = (i + 1) * divisor;
		return d <= large_ ? d - 1 : size() - n_ / d;
	}
};

#line 6 "components/multiplicative_sum_backend.hpp"
#include <algorithm>
#line 10 "components/multiplicative_sum_backend.hpp"

namespace multiplicative_sum_detail {
template<class T> struct LinearWeights {
	T count, weighted;
	template<class I> explicit LinearWeights(I value) : count(T(value)), weighted(T(value)) {}
	LinearWeights(T count_, T weighted_) : count(count_), weighted(weighted_) {}
	LinearWeights& operator+=(const LinearWeights& b) { count += b.count; weighted += b.weighted; return *this; }
	LinearWeights& operator-=(const LinearWeights& b) { count -= b.count; weighted -= b.weighted; return *this; }
	friend LinearWeights operator+(LinearWeights a, const LinearWeights& b) { return a += b; }
	friend LinearWeights operator-(LinearWeights a, const LinearWeights& b) { return a -= b; }
	friend LinearWeights operator*(const LinearWeights& a, const LinearWeights& b) { return LinearWeights(a.count * b.count, a.weighted * b.weighted); }
};

template<class T> struct LinearPrime {
	T weighted;
	explicit LinearPrime(const LinearWeights<T>& value) : weighted(value.weighted) {}
	explicit LinearPrime(T value) : weighted(value) {}
	operator LinearWeights<T>() const { return LinearWeights<T>(T(1), weighted); }
	friend LinearPrime operator*(const LinearPrime& a, const LinearPrime& b) { return LinearPrime(a.weighted * b.weighted); }
	friend LinearWeights<T> operator*(const LinearPrime& a, const LinearWeights<T>& b) { return LinearWeights<T>(b.count, a.weighted * b.weighted); }
	friend LinearWeights<T> operator*(const LinearWeights<T>& a, const LinearPrime& b) { return LinearWeights<T>(a.count, a.weighted * b.weighted); }
};
template<class T> struct PrimeValue { using type = T; };
template<class T> struct PrimeValue<LinearWeights<T>> { using type = LinearPrime<T>; };

template<class T> class Backend {
	using U = std::uint64_t;
	U n, root, large;
	std::size_t length;
	std::vector<U> primes;
	std::vector<typename PrimeValue<T>::type> first;
	std::vector<T> corrections;
	std::vector<std::size_t> offsets;
	static U square_root(U value) {
		U result = std::sqrt(static_cast<long double>(value));
		while (__uint128_t(result) * result > value) --result;
		while (__uint128_t(result + 1) * (result + 1) <= value) ++result;
		return result;
	}
	static U cube_root(U value) {
		U result = std::cbrt(static_cast<long double>(value));
		while (__uint128_t(result) * result * result > value) --result;
		while (__uint128_t(result + 1) * (result + 1) * (result + 1) <= value) ++result;
		return result;
	}
	static U fifth_root(U value) {
		U result = std::pow(static_cast<long double>(value), 0.2L);
		auto power = [](U x) { const __uint128_t square = __uint128_t(x) * x; return square * square * x; };
		while (power(result) > value) --result;
		while (power(result + 1) <= value) ++result;
		return result;
	}
	std::size_t small(U x) const { return x - 1; }
	std::size_t big(U d) const { return length - d; }
	std::size_t pi(U x) const { return std::upper_bound(primes.begin(), primes.end() - 1, x) - primes.begin(); }
	struct Tree {
		std::vector<T> data;
		explicit Tree(const std::vector<T>& a, std::size_t size) : data(size + 1, T(0)) {
			T previous(0);
			for (std::size_t i = 1; i <= size; ++i) { data[i] = a[i - 1] - previous; previous = a[i - 1]; }
			for (std::size_t i = 1; i <= size; ++i) if (i + (i & -i) <= size) data[i + (i & -i)] += data[i];
		}
		void add(std::size_t i, T value) { for (++i; i < data.size(); i += i & -i) data[i] += value; }
		T sum(std::size_t i) const { T result(0); for (; i; i -= i & -i) result += data[i]; return result; }
		void write(std::vector<T>& a) {
			for (std::size_t i = data.size() - 1; i; --i) if (i + (i & -i) < data.size()) data[i + (i & -i)] -= data[i];
			T result(0); for (std::size_t i = 1; i < data.size(); ++i) { result += data[i]; a[i - 1] = result; }
		}
	};
	T powerful(const std::vector<T>& base, std::size_t begin, U product, U limit) const {
		T answer = base[product <= large ? big(product) : small(limit)];
		for (std::size_t i = begin; i < first.size(); ++i) {
			const U p = primes[i];
			U power = p * p;
			if (power > limit) break;
			U next_limit = limit / power;
			std::size_t correction = offsets[i];
			for (;;) {
				const T weight = corrections[correction++];
				if constexpr (requires { weight == T(0); }) if (weight == T(0)) {
					if (next_limit < p) break;
					power *= p;
					next_limit /= p;
					continue;
				}
				if (i + 1 < first.size() && primes[i + 1] * primes[i + 1] <= next_limit)
					answer += weight * powerful(base, i + 1, product * power, next_limit);
				else answer += weight * base[product * power <= large ? big(product * power) : small(next_limit)];
				if (next_limit < p) break;
				power *= p;
				next_limit /= p;
			}
		}
		return answer;
	}
	template<bool SquareFree> void enumerate(Tree& tree, std::size_t begin, U product, T coefficient, U bound, unsigned depth = 0) const {
		if (depth == 3) return;
		for (std::size_t i = begin; i < first.size() && primes[i] <= bound / product; ++i) {
			const U next = product * primes[i];
			const T value = coefficient * first[i];
			tree.add(next <= root ? small(next) : big(n / next), value);
			enumerate<SquareFree>(tree, i + SquareFree, next, value, bound, depth + 1);
		}
	}

template<bool SquareFree> std::vector<T> build(const std::vector<T>* prefix, T inv3) const {
		const T before = SquareFree ? (*prefix)[large] : T(0);
		auto prime_large = [&](U d) {
			if constexpr (SquareFree) return (*prefix)[d - 1] - before;
			else return T(0);
		};
		const U tree_large = std::min<U>(std::floor(std::pow(n, 3.0 / 8)), large), tree_bound = n / (tree_large + 1);
		const U min_cut = std::max<U>(tree_bound / (root + 1), fifth_root(tree_bound));
		const U cut = std::max<U>(square_root(square_root(root)), min_cut);
		const auto pi_cut = pi(cut), pi_fourth = pi(square_root(root));
		std::vector<T> out(length, T(0));
		out[0] = T(1);
		for (std::size_t i = pi_cut; i < first.size(); ++i) out[small(primes[i])] = first[i];
		for (U i = 1; i < root; ++i) out[i] += out[i - 1];
		for (std::size_t i = pi_fourth, r = first.size(); i < first.size(); ++i) {
			const U p = primes[i], limit = n / p;
			while (r && U(r - 1) * primes[r - 1] > limit) --r;
			for (std::size_t j = i + SquareFree; j < r; ++j) out[big(limit / primes[j])] += first[i] * first[j];
			const U bound = std::min(primes[std::max(r, i + SquareFree)], root + 1);
			if (p <= large) out[big(large / p)] += first[i] * out[small(root)];
			if (bound <= limit) out[big(limit / bound)] -= first[i] * out[small(bound - 1)];
		}
		for (U i = large; i-- > 1;) out[big(i)] += out[big(i + 1)];
		for (std::size_t i = pi_fourth, r = first.size(); i < first.size(); ++i) {
			const U p = primes[i], limit = n / p;
			while (r && U(r - 1) * primes[r - 1] > limit) --r;
			const U bound = std::min(primes[std::max(r, i + SquareFree)], root + 1), split = large / p, end = limit / bound;
			for (U j = end; j > split; --j) out[big(j)] += first[i] * out[small(limit / j)];
			for (U j = split; j; --j) out[big(j)] += first[i] * prime_large(j * p);
		}
		const U bound_prime = primes[pi_fourth], bound_z = n / bound_prime / bound_prime / bound_prime;
		for (U i = 1; i <= bound_z; ++i) {
			const U limit = n / i;
			const auto bound1 = pi(large / i), bound2 = pi(square_root(limit / bound_prime)), bound3 = pi(cube_root(limit));
			T value(0);
			for (std::size_t j = pi_fourth; j < bound1; ++j) value += first[j] * out[big(i * primes[j])];
			for (std::size_t j = pi_fourth; j < bound2; ++j) {
				const T correction = first[j] * first[j] * (out[small(limit / primes[j] / primes[j])] - out[small(bound_prime - 1)]);
				if constexpr (SquareFree) value -= correction;
				else value += correction;
			}
			for (std::size_t j = pi_fourth; j < bound3; ++j) value += first[j] * first[j] * first[j];
			out[big(i)] += value * inv3;
		}
		for (U i = 1; i <= large; ++i) out[big(i)] += prime_large(i) + out[small(root)];
		const std::size_t tree_length = length - tree_large;
		Tree tree(out, tree_length);
		for (std::size_t i = pi_fourth; i-- > pi_cut;) {
			const U p = primes[i], limit = n / p, limit2 = limit / p;
			const U split = std::min(large / p, tree_large), end = std::min(limit2, tree_large);
			const U index2 = p * p <= root ? small(p * p) : std::min(big(limit2), tree_length);
			auto access = [&](U index) { return index2 <= index && index < tree_length ? tree.sum(index + 1) : out[index]; };
			if constexpr (SquareFree) {
				const T previous = out[small(p)];
				for (U j = 1; j <= split; ++j) out[big(j)] += first[i] * (access(big(j * p)) - previous);
				for (U j = split + 1; j <= end; ++j) out[big(j)] += first[i] * (access(small(limit / j)) - previous);
				enumerate<true>(tree, i + 1, p, first[i], tree_bound);
			} else {
				const T previous = out[small(p - 1)];
				enumerate<false>(tree, i, p, first[i], tree_bound);
				for (U j = end; j > split; --j) out[big(j)] += first[i] * (access(small(limit / j)) - previous);
				for (U j = split; j; --j) out[big(j)] += first[i] * (access(big(j * p)) - previous);
			}
		}
		tree.write(out);
		for (std::size_t i = pi_cut; i--;) {
			const U p = primes[i], limit = n / p, split = large / p;
			if constexpr (SquareFree) {
				for (U j = 1; j <= split; ++j) out[big(j)] += first[i] * out[big(j * p)];
				for (U j = split + 1; j <= large; ++j) out[big(j)] += first[i] * out[small(limit / j)];
				U j = root;
				for (U q = root / p; q; --q) {
					const T delta = first[i] * out[small(q)];
					for (const U stop = q * p; j >= stop; --j) out[small(j)] += delta;
				}
			} else {
				U j = p;
				for (U q = 1; j <= root; ++q) {
					const T delta = first[i] * out[small(q)];
					for (const U stop = std::min(root + 1, (q + 1) * p); j < stop; ++j) out[small(j)] += delta;
				}
				for (U j = large; j > split; --j) out[big(j)] += first[i] * out[small(limit / j)];
				for (U j = split; j; --j) out[big(j)] += first[i] * out[big(j * p)];
			}
		}
		return out;
	}
  public:
	Backend(const FloorQuotients& quotients, const std::vector<std::size_t>& prime_list) : n(quotients.max_n()), root(quotients.root()), large(quotients.large_size()), length(quotients.size()), primes(prime_list.begin(), prime_list.end()) { primes.push_back(root + 1); }
	template<class F> T sum(const std::vector<T>& prefix, F& f, T inv3) {
		if (!n) return T(0);
		first.clear();
		first.reserve(primes.size() - 1);
		corrections.clear();
		corrections.reserve(primes.size());
		offsets.clear();
		offsets.reserve(primes.size() - 1);
		for (std::size_t i = 0; i + 1 < primes.size(); ++i) {
			const U p = primes[i];
			first.emplace_back(f(p, 1, p));
			offsets.push_back(corrections.size());
			T previous(0);
			U power = p * p, limit = n / power;
			for (unsigned e = 2;; ++e) {
				previous = f(p, e, power) - first[i] * previous;
				corrections.push_back(previous);
				if (limit < p) break;
				power *= p;
				limit /= p;
			}
		}
		auto out = build<true>(&prefix, inv3);
		return powerful(out, 0, 1, n);
	}
	template<class Prefix> std::vector<T> prime_sum(Prefix& integer_prefix, T inv3) {
		if (!n) return {};
		first.clear();
		for (std::size_t i = 0; i + 1 < primes.size(); ++i) first.emplace_back(integer_prefix(primes[i]) - integer_prefix(primes[i] - 1));
		const auto base = build<false>(nullptr, inv3);
		std::vector<T> result(length, T(0));
		for (U i = 1; i <= large; ++i) result[i - 1] = integer_prefix(n / i) - base[big(i)];
		for (std::size_t j = 0; j < first.size() && primes[j] <= large; ++j) {
			const U p = primes[j], end = large / p;
			for (U i = 1; i <= end; ++i) result[i - 1] -= first[j] * result[i * p - 1];
		}
		T cumulative(0);
		std::size_t j = 0;
		for (U i = 1; i <= root; ++i) {
			if (i == primes[j]) cumulative += first[j++];
			result[length - i] = cumulative;
		}
		for (U i = 0; i < large; ++i) result[i] += cumulative;
		return result;
	}
};
}

#ifndef MULTIPLICATIVE_SUM_BACKEND_HPP
#define MULTIPLICATIVE_SUM_BACKEND_HPP

#include "floor_quotients.hpp"
#include <cstddef>
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <vector>

namespace multiplicative_sum_detail {
template<class T> struct LinearWeights {
	T count, weighted;
	template<class I> explicit LinearWeights(I value) : count(T(value)), weighted(T(value)) {}
	LinearWeights(T count_, T weighted_) : count(count_), weighted(weighted_) {}
	LinearWeights& operator+=(const LinearWeights& b) { count += b.count; weighted += b.weighted; return *this; }
	LinearWeights& operator-=(const LinearWeights& b) { count -= b.count; weighted -= b.weighted; return *this; }
	friend LinearWeights operator+(LinearWeights a, const LinearWeights& b) { return a += b; }
	friend LinearWeights operator-(LinearWeights a, const LinearWeights& b) { return a -= b; }
	friend LinearWeights operator*(const LinearWeights& a, const LinearWeights& b) { return LinearWeights(a.count * b.count, a.weighted * b.weighted); }
};

template<class T> struct LinearPrime {
	T weighted;
	explicit LinearPrime(const LinearWeights<T>& value) : weighted(value.weighted) {}
	explicit LinearPrime(T value) : weighted(value) {}
	operator LinearWeights<T>() const { return LinearWeights<T>(T(1), weighted); }
	friend LinearPrime operator*(const LinearPrime& a, const LinearPrime& b) { return LinearPrime(a.weighted * b.weighted); }
	friend LinearWeights<T> operator*(const LinearPrime& a, const LinearWeights<T>& b) { return LinearWeights<T>(b.count, a.weighted * b.weighted); }
	friend LinearWeights<T> operator*(const LinearWeights<T>& a, const LinearPrime& b) { return LinearWeights<T>(a.count, a.weighted * b.weighted); }
};
template<class T> struct PrimeValue { using type = T; };
template<class T> struct PrimeValue<LinearWeights<T>> { using type = LinearPrime<T>; };

template<class T> class Backend {
	using U = std::uint64_t;
	U n, root, large;
	std::size_t length;
	std::vector<U> primes;
	std::vector<typename PrimeValue<T>::type> first;
	std::vector<T> corrections;
	std::vector<std::size_t> offsets;
	static U square_root(U value) {
		U result = std::sqrt(static_cast<long double>(value));
		while (__uint128_t(result) * result > value) --result;
		while (__uint128_t(result + 1) * (result + 1) <= value) ++result;
		return result;
	}
	static U cube_root(U value) {
		U result = std::cbrt(static_cast<long double>(value));
		while (__uint128_t(result) * result * result > value) --result;
		while (__uint128_t(result + 1) * (result + 1) * (result + 1) <= value) ++result;
		return result;
	}
	static U fifth_root(U value) {
		U result = std::pow(static_cast<long double>(value), 0.2L);
		auto power = [](U x) { const __uint128_t square = __uint128_t(x) * x; return square * square * x; };
		while (power(result) > value) --result;
		while (power(result + 1) <= value) ++result;
		return result;
	}
	std::size_t small(U x) const { return x - 1; }
	std::size_t big(U d) const { return length - d; }
	std::size_t pi(U x) const { return std::upper_bound(primes.begin(), primes.end() - 1, x) - primes.begin(); }
	struct Tree {
		std::vector<T> data;
		explicit Tree(const std::vector<T>& a, std::size_t size) : data(size + 1, T(0)) {
			T previous(0);
			for (std::size_t i = 1; i <= size; ++i) { data[i] = a[i - 1] - previous; previous = a[i - 1]; }
			for (std::size_t i = 1; i <= size; ++i) if (i + (i & -i) <= size) data[i + (i & -i)] += data[i];
		}
		void add(std::size_t i, T value) { for (++i; i < data.size(); i += i & -i) data[i] += value; }
		T sum(std::size_t i) const { T result(0); for (; i; i -= i & -i) result += data[i]; return result; }
		void write(std::vector<T>& a) {
			for (std::size_t i = data.size() - 1; i; --i) if (i + (i & -i) < data.size()) data[i + (i & -i)] -= data[i];
			T result(0); for (std::size_t i = 1; i < data.size(); ++i) { result += data[i]; a[i - 1] = result; }
		}
	};
	T powerful(const std::vector<T>& base, std::size_t begin, U product, U limit) const {
		T answer = base[product <= large ? big(product) : small(limit)];
		for (std::size_t i = begin; i < first.size(); ++i) {
			const U p = primes[i];
			U power = p * p;
			if (power > limit) break;
			U next_limit = limit / power;
			std::size_t correction = offsets[i];
			for (;;) {
				const T weight = corrections[correction++];
				if constexpr (requires { weight == T(0); }) if (weight == T(0)) {
					if (next_limit < p) break;
					power *= p;
					next_limit /= p;
					continue;
				}
				if (i + 1 < first.size() && primes[i + 1] * primes[i + 1] <= next_limit)
					answer += weight * powerful(base, i + 1, product * power, next_limit);
				else answer += weight * base[product * power <= large ? big(product * power) : small(next_limit)];
				if (next_limit < p) break;
				power *= p;
				next_limit /= p;
			}
		}
		return answer;
	}
	template<bool SquareFree> void enumerate(Tree& tree, std::size_t begin, U product, T coefficient, U bound, unsigned depth = 0) const {
		if (depth == 3) return;
		for (std::size_t i = begin; i < first.size() && primes[i] <= bound / product; ++i) {
			const U next = product * primes[i];
			const T value = coefficient * first[i];
			tree.add(next <= root ? small(next) : big(n / next), value);
			enumerate<SquareFree>(tree, i + SquareFree, next, value, bound, depth + 1);
		}
	}

	template<bool SquareFree> std::vector<T> build(const std::vector<T>* prefix, T inv3) const {
		const T before = SquareFree ? (*prefix)[large] : T(0);
		auto prime_large = [&](U d) {
			if constexpr (SquareFree) return (*prefix)[d - 1] - before;
			else return T(0);
		};
		const U tree_large = std::min<U>(std::floor(std::pow(n, 3.0 / 8)), large), tree_bound = n / (tree_large + 1);
		const U min_cut = std::max<U>(tree_bound / (root + 1), fifth_root(tree_bound));
		const U cut = std::max<U>(square_root(square_root(root)), min_cut);
		const auto pi_cut = pi(cut), pi_fourth = pi(square_root(root));
		std::vector<T> out(length, T(0));
		out[0] = T(1);
		for (std::size_t i = pi_cut; i < first.size(); ++i) out[small(primes[i])] = first[i];
		for (U i = 1; i < root; ++i) out[i] += out[i - 1];
		for (std::size_t i = pi_fourth, r = first.size(); i < first.size(); ++i) {
			const U p = primes[i], limit = n / p;
			while (r && U(r - 1) * primes[r - 1] > limit) --r;
			for (std::size_t j = i + SquareFree; j < r; ++j) out[big(limit / primes[j])] += first[i] * first[j];
			const U bound = std::min(primes[std::max(r, i + SquareFree)], root + 1);
			if (p <= large) out[big(large / p)] += first[i] * out[small(root)];
			if (bound <= limit) out[big(limit / bound)] -= first[i] * out[small(bound - 1)];
		}
		for (U i = large; i-- > 1;) out[big(i)] += out[big(i + 1)];
		for (std::size_t i = pi_fourth, r = first.size(); i < first.size(); ++i) {
			const U p = primes[i], limit = n / p;
			while (r && U(r - 1) * primes[r - 1] > limit) --r;
			const U bound = std::min(primes[std::max(r, i + SquareFree)], root + 1), split = large / p, end = limit / bound;
			for (U j = end; j > split; --j) out[big(j)] += first[i] * out[small(limit / j)];
			for (U j = split; j; --j) out[big(j)] += first[i] * prime_large(j * p);
		}
		const U bound_prime = primes[pi_fourth], bound_z = n / bound_prime / bound_prime / bound_prime;
		for (U i = 1; i <= bound_z; ++i) {
			const U limit = n / i;
			const auto bound1 = pi(large / i), bound2 = pi(square_root(limit / bound_prime)), bound3 = pi(cube_root(limit));
			T value(0);
			for (std::size_t j = pi_fourth; j < bound1; ++j) value += first[j] * out[big(i * primes[j])];
			for (std::size_t j = pi_fourth; j < bound2; ++j) {
				const T correction = first[j] * first[j] * (out[small(limit / primes[j] / primes[j])] - out[small(bound_prime - 1)]);
				if constexpr (SquareFree) value -= correction;
				else value += correction;
			}
			for (std::size_t j = pi_fourth; j < bound3; ++j) value += first[j] * first[j] * first[j];
			out[big(i)] += value * inv3;
		}
		for (U i = 1; i <= large; ++i) out[big(i)] += prime_large(i) + out[small(root)];
		const std::size_t tree_length = length - tree_large;
		Tree tree(out, tree_length);
		for (std::size_t i = pi_fourth; i-- > pi_cut;) {
			const U p = primes[i], limit = n / p, limit2 = limit / p;
			const U split = std::min(large / p, tree_large), end = std::min(limit2, tree_large);
			const U index2 = p * p <= root ? small(p * p) : std::min(big(limit2), tree_length);
			auto access = [&](U index) { return index2 <= index && index < tree_length ? tree.sum(index + 1) : out[index]; };
			if constexpr (SquareFree) {
				const T previous = out[small(p)];
				for (U j = 1; j <= split; ++j) out[big(j)] += first[i] * (access(big(j * p)) - previous);
				for (U j = split + 1; j <= end; ++j) out[big(j)] += first[i] * (access(small(limit / j)) - previous);
				enumerate<true>(tree, i + 1, p, first[i], tree_bound);
			} else {
				const T previous = out[small(p - 1)];
				enumerate<false>(tree, i, p, first[i], tree_bound);
				for (U j = end; j > split; --j) out[big(j)] += first[i] * (access(small(limit / j)) - previous);
				for (U j = split; j; --j) out[big(j)] += first[i] * (access(big(j * p)) - previous);
			}
		}
		tree.write(out);
		for (std::size_t i = pi_cut; i--;) {
			const U p = primes[i], limit = n / p, split = large / p;
			if constexpr (SquareFree) {
				for (U j = 1; j <= split; ++j) out[big(j)] += first[i] * out[big(j * p)];
				for (U j = split + 1; j <= large; ++j) out[big(j)] += first[i] * out[small(limit / j)];
				U j = root;
				for (U q = root / p; q; --q) {
					const T delta = first[i] * out[small(q)];
					for (const U stop = q * p; j >= stop; --j) out[small(j)] += delta;
				}
			} else {
				U j = p;
				for (U q = 1; j <= root; ++q) {
					const T delta = first[i] * out[small(q)];
					for (const U stop = std::min(root + 1, (q + 1) * p); j < stop; ++j) out[small(j)] += delta;
				}
				for (U j = large; j > split; --j) out[big(j)] += first[i] * out[small(limit / j)];
				for (U j = split; j; --j) out[big(j)] += first[i] * out[big(j * p)];
			}
		}
		return out;
	}
  public:
	Backend(const FloorQuotients& quotients, const std::vector<std::size_t>& prime_list) : n(quotients.max_n()), root(quotients.root()), large(quotients.large_size()), length(quotients.size()), primes(prime_list.begin(), prime_list.end()) { primes.push_back(root + 1); }
	template<class F> T sum(const std::vector<T>& prefix, F& f, T inv3) {
		if (!n) return T(0);
		first.clear();
		first.reserve(primes.size() - 1);
		corrections.clear();
		corrections.reserve(primes.size());
		offsets.clear();
		offsets.reserve(primes.size() - 1);
		for (std::size_t i = 0; i + 1 < primes.size(); ++i) {
			const U p = primes[i];
			first.emplace_back(f(p, 1, p));
			offsets.push_back(corrections.size());
			T previous(0);
			U power = p * p, limit = n / power;
			for (unsigned e = 2;; ++e) {
				previous = f(p, e, power) - first[i] * previous;
				corrections.push_back(previous);
				if (limit < p) break;
				power *= p;
				limit /= p;
			}
		}
		auto out = build<true>(&prefix, inv3);
		return powerful(out, 0, 1, n);
	}
	template<class Prefix> std::vector<T> prime_sum(Prefix& integer_prefix, T inv3) {
		if (!n) return {};
		first.clear();
		for (std::size_t i = 0; i + 1 < primes.size(); ++i) first.emplace_back(integer_prefix(primes[i]) - integer_prefix(primes[i] - 1));
		const auto base = build<false>(nullptr, inv3);
		std::vector<T> result(length, T(0));
		for (U i = 1; i <= large; ++i) result[i - 1] = integer_prefix(n / i) - base[big(i)];
		for (std::size_t j = 0; j < first.size() && primes[j] <= large; ++j) {
			const U p = primes[j], end = large / p;
			for (U i = 1; i <= end; ++i) result[i - 1] -= first[j] * result[i * p - 1];
		}
		T cumulative(0);
		std::size_t j = 0;
		for (U i = 1; i <= root; ++i) {
			if (i == primes[j]) cumulative += first[j++];
			result[length - i] = cumulative;
		}
		for (U i = 0; i < large; ++i) result[i] += cumulative;
		return result;
	}
};
}

#endif

#line 1 "components/multiplicative_sum_backend.hpp"



#line 1 "components/floor_quotients.hpp"



#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <iterator>
#include <utility>
#include <vector>

namespace floor_quotient_detail {
inline std::uint64_t square_root(std::uint64_t n) {
	std::uint64_t root = std::sqrt(static_cast<long double>(n));
	while (__uint128_t(root) * root > n) --root;
	while (__uint128_t(root + 1) * (root + 1) <= n) ++root;
	return root;
}
inline std::uint64_t divide(std::uint64_t n, std::uint64_t divisor) {
	if constexpr (std::numeric_limits<double>::is_iec559 && std::numeric_limits<double>::radix == 2 && std::numeric_limits<double>::digits >= 53) {
		if (n <= (std::uint64_t(1) << 53)) {
			std::uint64_t value = double(n) / double(divisor);
			return value - (value * divisor > n);
		}
	}
	return n / divisor;
}
}

class FloorQuotientRange {
	std::uint64_t n_, root_, large_;

  public:
	explicit FloorQuotientRange(std::uint64_t n = 0) : n_(n), root_(floor_quotient_detail::square_root(n)), large_(root_ - (n < root_ * (root_ + 1))) {}
	std::uint64_t max_n() const { return n_; }
	std::uint64_t root() const { return root_; }
	std::uint64_t large_size() const { return large_; }
	std::uint64_t size() const { return root_ + large_; }
	bool empty() const { return !n_; }
	std::uint64_t operator[](std::uint64_t i) const {
		assert(i < size());
		return i < root_ ? i + 1 : floor_quotient_detail::divide(n_, size() - i);
	}
	class iterator {
		const FloorQuotientRange *range_ = nullptr;
		std::uint64_t index_ = 0;
		friend class FloorQuotientRange;
		iterator(const FloorQuotientRange *range, std::uint64_t i) : range_(range), index_(i) {}

	  public:
		using value_type = std::uint64_t;
		using difference_type = std::int64_t;
		using iterator_category = std::input_iterator_tag;
		using iterator_concept = std::forward_iterator_tag;
		using reference = value_type;
		using pointer = void;
		iterator() = default;
		value_type operator*() const { return (*range_)[index_]; }
		iterator &operator++() { ++index_; return *this; }
		iterator operator++(int) { auto result = *this; ++*this; return result; }
		bool operator==(const iterator &) const = default;
	};
	iterator begin() const { return iterator(this, 0); }
	iterator end() const { return iterator(this, size()); }
};

class FloorQuotients {
	std::uint64_t n_, root_ = 0;
	std::size_t large_ = 0;
	std::vector<std::uint64_t> values_;

  public:
	explicit FloorQuotients(std::uint64_t n = 0) : n_(n) {
		if (!n) return;
		root_ = floor_quotient_detail::square_root(n);
		large_ = n / (root_ + 1);
		values_.reserve(large_ + root_);
		for (std::size_t i = 1; i <= large_; ++i) values_.push_back(floor_quotient_detail::divide(n, i));
		for (auto x = root_; x; --x) values_.push_back(x);
	}
	FloorQuotients(const FloorQuotients&) = default;
	FloorQuotients& operator=(const FloorQuotients&) = default;
	FloorQuotients(FloorQuotients&& other) noexcept : FloorQuotients() { swap(other); }
	FloorQuotients& operator=(FloorQuotients&& other) noexcept {
		if (this != &other) { FloorQuotients empty; swap(empty); swap(other); }
		return *this;
	}
	void swap(FloorQuotients& other) noexcept {
		std::swap(n_, other.n_), std::swap(root_, other.root_), std::swap(large_, other.large_);
		values_.swap(other.values_);
	}
	std::uint64_t max_n() const { return n_; }
	std::uint64_t root() const { return root_; }
	std::size_t size() const { return values_.size(); }
	std::size_t large_size() const { return large_; }
	const std::vector<std::uint64_t>& values() const { return values_; }
	bool contains(std::uint64_t x) const { return x && x <= n_ && (x <= root_ || n_ / (n_ / x) == x); }
	std::size_t index(std::uint64_t x) const {
		assert(x && x <= n_);
		return x <= root_ ? size() - x : n_ / x - 1;
	}
	std::size_t divided_index(std::size_t i, std::uint64_t divisor) const {
		assert(i < size() && divisor && divisor <= values_[i]);
		if (i >= large_) return size() - values_[i] / divisor;
		const auto d = (i + 1) * divisor;
		return d <= large_ ? d - 1 : size() - n_ / d;
	}
};


#line 6 "components/multiplicative_sum_backend.hpp"
#include <algorithm>
#line 10 "components/multiplicative_sum_backend.hpp"

namespace multiplicative_sum_detail {
template<class T> struct LinearWeights {
	T count, weighted;
	template<class I> explicit LinearWeights(I value) : count(T(value)), weighted(T(value)) {}
	LinearWeights(T count_, T weighted_) : count(count_), weighted(weighted_) {}
	LinearWeights& operator+=(const LinearWeights& b) { count += b.count; weighted += b.weighted; return *this; }
	LinearWeights& operator-=(const LinearWeights& b) { count -= b.count; weighted -= b.weighted; return *this; }
	friend LinearWeights operator+(LinearWeights a, const LinearWeights& b) { return a += b; }
	friend LinearWeights operator-(LinearWeights a, const LinearWeights& b) { return a -= b; }
	friend LinearWeights operator*(const LinearWeights& a, const LinearWeights& b) { return LinearWeights(a.count * b.count, a.weighted * b.weighted); }
};

template<class T> struct LinearPrime {
	T weighted;
	explicit LinearPrime(const LinearWeights<T>& value) : weighted(value.weighted) {}
	explicit LinearPrime(T value) : weighted(value) {}
	operator LinearWeights<T>() const { return LinearWeights<T>(T(1), weighted); }
	friend LinearPrime operator*(const LinearPrime& a, const LinearPrime& b) { return LinearPrime(a.weighted * b.weighted); }
	friend LinearWeights<T> operator*(const LinearPrime& a, const LinearWeights<T>& b) { return LinearWeights<T>(b.count, a.weighted * b.weighted); }
	friend LinearWeights<T> operator*(const LinearWeights<T>& a, const LinearPrime& b) { return LinearWeights<T>(a.count, a.weighted * b.weighted); }
};
template<class T> struct PrimeValue { using type = T; };
template<class T> struct PrimeValue<LinearWeights<T>> { using type = LinearPrime<T>; };

template<class T> class Backend {
	using U = std::uint64_t;
	U n, root, large;
	std::size_t length;
	std::vector<U> primes;
	std::vector<typename PrimeValue<T>::type> first;
	std::vector<T> corrections;
	std::vector<std::size_t> offsets;
	static U square_root(U value) {
		U result = std::sqrt(static_cast<long double>(value));
		while (__uint128_t(result) * result > value) --result;
		while (__uint128_t(result + 1) * (result + 1) <= value) ++result;
		return result;
	}
	static U cube_root(U value) {
		U result = std::cbrt(static_cast<long double>(value));
		while (__uint128_t(result) * result * result > value) --result;
		while (__uint128_t(result + 1) * (result + 1) * (result + 1) <= value) ++result;
		return result;
	}
	static U fifth_root(U value) {
		U result = std::pow(static_cast<long double>(value), 0.2L);
		auto power = [](U x) { const __uint128_t square = __uint128_t(x) * x; return square * square * x; };
		while (power(result) > value) --result;
		while (power(result + 1) <= value) ++result;
		return result;
	}
	std::size_t small(U x) const { return x - 1; }
	std::size_t big(U d) const { return length - d; }
	std::size_t pi(U x) const { return std::upper_bound(primes.begin(), primes.end() - 1, x) - primes.begin(); }
	struct Tree {
		std::vector<T> data;
		explicit Tree(const std::vector<T>& a, std::size_t size) : data(size + 1, T(0)) {
			T previous(0);
			for (std::size_t i = 1; i <= size; ++i) { data[i] = a[i - 1] - previous; previous = a[i - 1]; }
			for (std::size_t i = 1; i <= size; ++i) if (i + (i & -i) <= size) data[i + (i & -i)] += data[i];
		}
		void add(std::size_t i, T value) { for (++i; i < data.size(); i += i & -i) data[i] += value; }
		T sum(std::size_t i) const { T result(0); for (; i; i -= i & -i) result += data[i]; return result; }
		void write(std::vector<T>& a) {
			for (std::size_t i = data.size() - 1; i; --i) if (i + (i & -i) < data.size()) data[i + (i & -i)] -= data[i];
			T result(0); for (std::size_t i = 1; i < data.size(); ++i) { result += data[i]; a[i - 1] = result; }
		}
	};
	T powerful(const std::vector<T>& base, std::size_t begin, U product, U limit) const {
		T answer = base[product <= large ? big(product) : small(limit)];
		for (std::size_t i = begin; i < first.size(); ++i) {
			const U p = primes[i];
			U power = p * p;
			if (power > limit) break;
			U next_limit = limit / power;
			std::size_t correction = offsets[i];
			for (;;) {
				const T weight = corrections[correction++];
				if constexpr (requires { weight == T(0); }) if (weight == T(0)) {
					if (next_limit < p) break;
					power *= p;
					next_limit /= p;
					continue;
				}
				if (i + 1 < first.size() && primes[i + 1] * primes[i + 1] <= next_limit)
					answer += weight * powerful(base, i + 1, product * power, next_limit);
				else answer += weight * base[product * power <= large ? big(product * power) : small(next_limit)];
				if (next_limit < p) break;
				power *= p;
				next_limit /= p;
			}
		}
		return answer;
	}
	template<bool SquareFree> void enumerate(Tree& tree, std::size_t begin, U product, T coefficient, U bound, unsigned depth = 0) const {
		if (depth == 3) return;
		for (std::size_t i = begin; i < first.size() && primes[i] <= bound / product; ++i) {
			const U next = product * primes[i];
			const T value = coefficient * first[i];
			tree.add(next <= root ? small(next) : big(n / next), value);
			enumerate<SquareFree>(tree, i + SquareFree, next, value, bound, depth + 1);
		}
	}

	template<bool SquareFree> std::vector<T> build(const std::vector<T>* prefix, T inv3) const {
		const T before = SquareFree ? (*prefix)[large] : T(0);
		auto prime_large = [&](U d) {
			if constexpr (SquareFree) return (*prefix)[d - 1] - before;
			else return T(0);
		};
		const U tree_large = std::min<U>(std::floor(std::pow(n, 3.0 / 8)), large), tree_bound = n / (tree_large + 1);
		const U min_cut = std::max<U>(tree_bound / (root + 1), fifth_root(tree_bound));
		const U cut = std::max<U>(square_root(square_root(root)), min_cut);
		const auto pi_cut = pi(cut), pi_fourth = pi(square_root(root));
		std::vector<T> out(length, T(0));
		out[0] = T(1);
		for (std::size_t i = pi_cut; i < first.size(); ++i) out[small(primes[i])] = first[i];
		for (U i = 1; i < root; ++i) out[i] += out[i - 1];
		for (std::size_t i = pi_fourth, r = first.size(); i < first.size(); ++i) {
			const U p = primes[i], limit = n / p;
			while (r && U(r - 1) * primes[r - 1] > limit) --r;
			for (std::size_t j = i + SquareFree; j < r; ++j) out[big(limit / primes[j])] += first[i] * first[j];
			const U bound = std::min(primes[std::max(r, i + SquareFree)], root + 1);
			if (p <= large) out[big(large / p)] += first[i] * out[small(root)];
			if (bound <= limit) out[big(limit / bound)] -= first[i] * out[small(bound - 1)];
		}
		for (U i = large; i-- > 1;) out[big(i)] += out[big(i + 1)];
		for (std::size_t i = pi_fourth, r = first.size(); i < first.size(); ++i) {
			const U p = primes[i], limit = n / p;
			while (r && U(r - 1) * primes[r - 1] > limit) --r;
			const U bound = std::min(primes[std::max(r, i + SquareFree)], root + 1), split = large / p, end = limit / bound;
			for (U j = end; j > split; --j) out[big(j)] += first[i] * out[small(limit / j)];
			for (U j = split; j; --j) out[big(j)] += first[i] * prime_large(j * p);
		}
		const U bound_prime = primes[pi_fourth], bound_z = n / bound_prime / bound_prime / bound_prime;
		for (U i = 1; i <= bound_z; ++i) {
			const U limit = n / i;
			const auto bound1 = pi(large / i), bound2 = pi(square_root(limit / bound_prime)), bound3 = pi(cube_root(limit));
			T value(0);
			for (std::size_t j = pi_fourth; j < bound1; ++j) value += first[j] * out[big(i * primes[j])];
			for (std::size_t j = pi_fourth; j < bound2; ++j) {
				const T correction = first[j] * first[j] * (out[small(limit / primes[j] / primes[j])] - out[small(bound_prime - 1)]);
				if constexpr (SquareFree) value -= correction;
				else value += correction;
			}
			for (std::size_t j = pi_fourth; j < bound3; ++j) value += first[j] * first[j] * first[j];
			out[big(i)] += value * inv3;
		}
		for (U i = 1; i <= large; ++i) out[big(i)] += prime_large(i) + out[small(root)];
		const std::size_t tree_length = length - tree_large;
		Tree tree(out, tree_length);
		for (std::size_t i = pi_fourth; i-- > pi_cut;) {
			const U p = primes[i], limit = n / p, limit2 = limit / p;
			const U split = std::min(large / p, tree_large), end = std::min(limit2, tree_large);
			const U index2 = p * p <= root ? small(p * p) : std::min(big(limit2), tree_length);
			auto access = [&](U index) { return index2 <= index && index < tree_length ? tree.sum(index + 1) : out[index]; };
			if constexpr (SquareFree) {
				const T previous = out[small(p)];
				for (U j = 1; j <= split; ++j) out[big(j)] += first[i] * (access(big(j * p)) - previous);
				for (U j = split + 1; j <= end; ++j) out[big(j)] += first[i] * (access(small(limit / j)) - previous);
				enumerate<true>(tree, i + 1, p, first[i], tree_bound);
			} else {
				const T previous = out[small(p - 1)];
				enumerate<false>(tree, i, p, first[i], tree_bound);
				for (U j = end; j > split; --j) out[big(j)] += first[i] * (access(small(limit / j)) - previous);
				for (U j = split; j; --j) out[big(j)] += first[i] * (access(big(j * p)) - previous);
			}
		}
		tree.write(out);
		for (std::size_t i = pi_cut; i--;) {
			const U p = primes[i], limit = n / p, split = large / p;
			if constexpr (SquareFree) {
				for (U j = 1; j <= split; ++j) out[big(j)] += first[i] * out[big(j * p)];
				for (U j = split + 1; j <= large; ++j) out[big(j)] += first[i] * out[small(limit / j)];
				U j = root;
				for (U q = root / p; q; --q) {
					const T delta = first[i] * out[small(q)];
					for (const U stop = q * p; j >= stop; --j) out[small(j)] += delta;
				}
			} else {
				U j = p;
				for (U q = 1; j <= root; ++q) {
					const T delta = first[i] * out[small(q)];
					for (const U stop = std::min(root + 1, (q + 1) * p); j < stop; ++j) out[small(j)] += delta;
				}
				for (U j = large; j > split; --j) out[big(j)] += first[i] * out[small(limit / j)];
				for (U j = split; j; --j) out[big(j)] += first[i] * out[big(j * p)];
			}
		}
		return out;
	}
  public:
	Backend(const FloorQuotients& quotients, const std::vector<std::size_t>& prime_list) : n(quotients.max_n()), root(quotients.root()), large(quotients.large_size()), length(quotients.size()), primes(prime_list.begin(), prime_list.end()) { primes.push_back(root + 1); }
	template<class F> T sum(const std::vector<T>& prefix, F& f, T inv3) {
		if (!n) return T(0);
		first.clear();
		first.reserve(primes.size() - 1);
		corrections.clear();
		corrections.reserve(primes.size());
		offsets.clear();
		offsets.reserve(primes.size() - 1);
		for (std::size_t i = 0; i + 1 < primes.size(); ++i) {
			const U p = primes[i];
			first.emplace_back(f(p, 1, p));
			offsets.push_back(corrections.size());
			T previous(0);
			U power = p * p, limit = n / power;
			for (unsigned e = 2;; ++e) {
				previous = f(p, e, power) - first[i] * previous;
				corrections.push_back(previous);
				if (limit < p) break;
				power *= p;
				limit /= p;
			}
		}
		auto out = build<true>(&prefix, inv3);
		return powerful(out, 0, 1, n);
	}
	template<class Prefix> std::vector<T> prime_sum(Prefix& integer_prefix, T inv3) {
		if (!n) return {};
		first.clear();
		for (std::size_t i = 0; i + 1 < primes.size(); ++i) first.emplace_back(integer_prefix(primes[i]) - integer_prefix(primes[i] - 1));
		const auto base = build<false>(nullptr, inv3);
		std::vector<T> result(length, T(0));
		for (U i = 1; i <= large; ++i) result[i - 1] = integer_prefix(n / i) - base[big(i)];
		for (std::size_t j = 0; j < first.size() && primes[j] <= large; ++j) {
			const U p = primes[j], end = large / p;
			for (U i = 1; i <= end; ++i) result[i - 1] -= first[j] * result[i * p - 1];
		}
		T cumulative(0);
		std::size_t j = 0;
		for (U i = 1; i <= root; ++i) {
			if (i == primes[j]) cumulative += first[j++];
			result[length - i] = cumulative;
		}
		for (U i = 0; i < large; ++i) result[i] += cumulative;
		return result;
	}
};
}