#line 2 "multiplicative-function/count-square-free.hpp"
#include <unordered_map>
#include <vector>
using namespace std;
#line 2 "math/enumerate-quotient.hpp"
#line 2 "math/isqrt.hpp"
#include <cmath>
using namespace std;
// floor(sqrt(n)) を返す (ただし n が負の場合は 0 を返す)
long long isqrt(long long n) {
if (n <= 0) return 0;
long long x = sqrt(n);
while ((x + 1) * (x + 1) <= n) x++;
while (x * x > n) x--;
return x;
}
#line 4 "math/enumerate-quotient.hpp"
namespace EnumerateQuotientImpl {
long long fast_div(long long a, long long b) { return 1.0 * a / b; };
long long slow_div(long long a, long long b) { return a / b; };
} // namespace EnumerateQuotientImpl
// { (q, l, r) : forall x in (l,r], floor(N/x) = q }
// を引数に取る関数f(q, l, r)を渡す。範囲が左に半開なのに注意
// 商は小さい方から走査する
template <typename T, typename F>
void enumerate_quotient(T N, const F& f) {
T sq = isqrt(N);
#define FUNC(d) \
T upper = N, quo = 0; \
while (upper > sq) { \
T thres = d(N, (++quo + 1)); \
f(quo, thres, upper); \
upper = thres; \
} \
while (upper > 0) { \
f(d(N, upper), upper - 1, upper); \
upper--; \
}
if (N <= 1e12) {
FUNC(EnumerateQuotientImpl::fast_div);
} else {
FUNC(EnumerateQuotientImpl::slow_div);
}
#undef FUNC
}
/**
* @brief 商の列挙
*/
#line 2 "multiplicative-function/mf-famous-series.hpp"
#line 2 "multiplicative-function/enumerate-multiplicative-function.hpp"
#line 4 "multiplicative-function/enumerate-multiplicative-function.hpp"
using namespace std;
#line 2 "prime/prime-enumerate.hpp"
#line 5 "prime/prime-enumerate.hpp"
using namespace std;
// Prime Sieve {2, 3, 5, 7, 11, 13, 17, ...}
vector<int> prime_enumerate(int N) {
vector<bool> sieve(N / 3 + 1, 1);
for (int p = 5, d = 4, i = 1, sqn = sqrt(N); p <= sqn; p += d = 6 - d, i++) {
if (!sieve[i]) continue;
for (int q = p * p / 3, r = d * p / 3 + (d * p % 3 == 2), s = 2 * p,
qe = sieve.size();
q < qe; q += r = s - r)
sieve[q] = 0;
}
vector<int> ret{2, 3};
for (int p = 5, d = 4, i = 1; p <= N; p += d = 6 - d, i++)
if (sieve[i]) ret.push_back(p);
while (!ret.empty() && ret.back() > N) ret.pop_back();
return ret;
}
#line 7 "multiplicative-function/enumerate-multiplicative-function.hpp"
// f(n, p, c) : n = pow(p, c), f is multiplicative function
template <typename T, T (*f)(int, int, int)>
struct enumerate_multiplicative_function {
enumerate_multiplicative_function(int _n)
: ps(prime_enumerate(_n)), a(_n + 1, T()), n(_n), p(ps.size()) {}
vector<T> run() {
a[1] = 1;
dfs(-1, 1, 1);
return a;
}
private:
vector<int> ps;
vector<T> a;
int n, p;
void dfs(int i, long long x, T y) {
a[x] = y;
if (y == T()) return;
for (int j = i + 1; j < p; j++) {
long long nx = x * ps[j];
long long dx = ps[j];
if (nx > n) break;
for (int c = 1; nx <= n; nx *= ps[j], dx *= ps[j], ++c) {
dfs(j, nx, y * f(dx, ps[j], c));
}
}
}
};
template <typename T, T (*f)(int, int, int)>
using enamurate_multiplicative_function =
enumerate_multiplicative_function<T, f>;
/**
* @brief 乗法的関数の列挙
*/
#line 4 "multiplicative-function/mf-famous-series.hpp"
namespace multiplicative_function {
template <typename T>
T moebius(int, int, int c) {
return c == 0 ? 1 : c == 1 ? -1 : 0;
}
template <typename T>
T sigma0(int, int, int c) {
return c + 1;
}
template <typename T>
T sigma1(int n, int p, int) {
return (n - 1) / (p - 1) + n;
}
template <typename T>
T totient(int n, int p, int) {
return n - n / p;
}
} // namespace multiplicative_function
template <typename T>
static constexpr vector<T> mobius_function(int n) {
enumerate_multiplicative_function<T, multiplicative_function::moebius<T>> em(
n);
return em.run();
}
template <typename T>
static constexpr vector<T> sigma0(int n) {
enumerate_multiplicative_function<T, multiplicative_function::sigma0<T>> em(
n);
return em.run();
}
template <typename T>
static constexpr vector<T> sigma1(int n) {
enumerate_multiplicative_function<T, multiplicative_function::sigma1<T>> em(
n);
return em.run();
}
template <typename T>
static constexpr vector<T> totient(int n) {
enumerate_multiplicative_function<T, multiplicative_function::totient<T>> em(
n);
return em.run();
}
/**
* @brief 有名な乗法的関数
*/
#line 10 "multiplicative-function/count-square-free.hpp"
long long count_square_free(long long N) {
long long B = pow(N, 0.4);
auto pre = mobius_function<int>(B + 1);
for (int i = 1; i <= B; i++) pre[i] += pre[i - 1];
unordered_map<long long, long long> dp;
auto mu = [&](auto rc, long long n) -> long long {
if (n <= B) return pre[n];
if (dp.count(n)) return dp[n];
long long cur = 1;
enumerate_quotient(n, [&](long long q, long long l, long long r) {
if (q != n) cur -= rc(rc, q) * (r - l);
});
return dp[n] = cur;
};
long long ans = 0;
long long upper = isqrt(N);
for (long long i = 1; upper > B; i++) {
long long nxt = isqrt(N / (i + 1));
ans += i * (mu(mu, upper) - mu(mu, nxt));
upper = nxt;
}
while (upper > 0) {
ans += (pre[upper] - pre[upper - 1]) * (N / (upper * upper));
upper--;
}
return ans;
}
/**
* @brief 無平方数の数え上げ
*/
無平方数の数え上げ