乗法的関数のprefix sum の列挙
(multiplicative-function/enumerate-sum-of-multiplicative-function.hpp)
Depends on
Verified with
Code
#pragma once
#include <functional>
#include <unordered_map>
using namespace std;
#include "../math/enumerate-quotient.hpp"
#include "../math/isqrt.hpp"
/**
* S(f, n) = f(1) + f(2) + ... + f(n) とする
* f と g のディリクレ積を h とする
* S(h, n) と S(g, n) が高速に計算できる, かつ g(1) = 1 のとき
* S(f, N/i) を O(N^{3/4}) で列挙できる
*
* うまくやると O~(N^{2/3}) に落ちたり g(1) != 1 に対応できる
* https://codeforces.com/blog/entry/54150
*/
template <typename T, typename SG, typename SH>
struct enumerate_mf_prefix_sum {
long long N, sq;
const SG sg;
const SH sh;
vector<T> small, large;
T& ref(long long x) {
if (x <= sq) {
return small[x];
} else if (N <= 1000000000000LL) {
return large[1.0 * N / x];
} else {
return large[N / x];
}
}
enumerate_mf_prefix_sum(long long _n, SG _sg, SH _sh)
: N(_n), sq(isqrt(N)), sg(_sg), sh(_sh) {
small.resize(sq + 1);
large.resize(sq + 1);
enumerate_quotient(N, [&](long long n, long long, long long) {
T& cur = (ref(n) = sh(n));
enumerate_quotient(n, [&](long long q, long long l, long long r) {
if (q != n) cur -= ref(q) * (sg(r) - sg(l));
});
});
}
T get(long long n) { return ref(n); }
T operator()(long long n) { return get(n); }
};
/**
* @brief 乗法的関数のprefix sum の列挙
*/
#line 2 "multiplicative-function/enumerate-sum-of-multiplicative-function.hpp"
#include <functional>
#include <unordered_map>
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 9 "multiplicative-function/enumerate-sum-of-multiplicative-function.hpp"
/**
* S(f, n) = f(1) + f(2) + ... + f(n) とする
* f と g のディリクレ積を h とする
* S(h, n) と S(g, n) が高速に計算できる, かつ g(1) = 1 のとき
* S(f, N/i) を O(N^{3/4}) で列挙できる
*
* うまくやると O~(N^{2/3}) に落ちたり g(1) != 1 に対応できる
* https://codeforces.com/blog/entry/54150
*/
template <typename T, typename SG, typename SH>
struct enumerate_mf_prefix_sum {
long long N, sq;
const SG sg;
const SH sh;
vector<T> small, large;
T& ref(long long x) {
if (x <= sq) {
return small[x];
} else if (N <= 1000000000000LL) {
return large[1.0 * N / x];
} else {
return large[N / x];
}
}
enumerate_mf_prefix_sum(long long _n, SG _sg, SH _sh)
: N(_n), sq(isqrt(N)), sg(_sg), sh(_sh) {
small.resize(sq + 1);
large.resize(sq + 1);
enumerate_quotient(N, [&](long long n, long long, long long) {
T& cur = (ref(n) = sh(n));
enumerate_quotient(n, [&](long long q, long long l, long long r) {
if (q != n) cur -= ref(q) * (sg(r) - sg(l));
});
});
}
T get(long long n) { return ref(n); }
T operator()(long long n) { return get(n); }
};
/**
* @brief 乗法的関数のprefix sum の列挙
*/
Back to top page