闲话 22.12.24

闲话

今晚 ABC！
看我是不是卡死在 B 题上（

jjdw；你这太菜了 不行啊

今天的推歌是 I LOST(2022Remaster) by はるまきごはん feat. 初音ミク
好最后还是回到了饭的歌呢！
这首歌我总是无意识地哼
好听好听！

\(\text{Powerful Number}\) 筛

还差一个筛子。以前一直没法理解，现在就班门弄斧说一说吧。

记 \(*\) 为狄利克雷卷积，\(\times\) 或隐式结合为数乘，\(+\) 为数加，\(S(f, n)\) 为对数论函数 \(f\) 求和，上界为 \(n\)。

\(\text{Powerful Number}\) 筛，简称 \(\text{PN}\) 筛，应用了特殊构造的狄利克雷卷积，将点值归类至 \(\text{PN}\) 处，从而加速求解。

首先介绍一下 \(\text{Powerful Number}\)，下面简称 \(\text{PN}\)。

记 \(n\) 的质因数分解为 \(\prod_{i} p_i^{c_i}\)。\(n\) 为 \(\text{PN}\) 当且仅当 \(\forall i,\ c_i > 1\)。容易发现任意 \(\text{PN}\) 可以表示为 \(a^2b^3\) 的形式：对于偶数的 \(c_i\) 直接合并入 \(a\)，其余情况取出 \(p_i^3\) 合入 \(b\)，再将剩余部分合入 \(a\)。

可以发现，\(n\) 以内的 \(\text{PN}\) 有 \(O(\sqrt n)\) 个。
证明：

\[\int_1^{\sqrt {n}} \sqrt[3]{\frac{n}{x^2}} \text dx = O(\sqrt n) \]

下面的操作中需要找到所有 \(n\) 以内的 \(\text{PN}\)。我们可以首先筛出 \(\sqrt n\) 以内的质数，随后 dfs 各质数的指数即可。由于不会产生无用搜索，总时间复杂度为 \(O(\sqrt n)\)。

假设我们需要求解积性函数 \(f\) 的前缀和 \(S(f, n)\)。应用 PN 筛的要求是存在一个易求前缀和的积性函数函数 \(g\)，满足 \(\forall p\text{ is prime},\ f(p) = g(p)\)。

我们构造这么一个 \(g\)，再构造 \(h\) 满足 \(f = g * h\)。容易发现 \(h\) 积性。观察 \(h\) 的取值。取素数 \(p\)，\(f(p) = g(1)h(p) + g(p)h(1) = h(p) + g(p) = h(p) + f(p)\)。因此 \(h\) 在任意质数处的取值均为 \(0\)。归纳得到 \(h\) 只在 \(\text{PN}\) 处有值。

根据 \(f = g*h\)，

\[\begin{aligned} & S(f, n) \\ = \ & \sum_{i=1}^n f(i) \\ = \ & \sum_{i=1}^n \sum_{d|i} h(d) g\left(\frac id\right) \\ = \ & \sum_{d=1}^n h(d) \sum_{i = 1}^{\lfloor n / d\rfloor} g(i) \\ = \ & \sum_{d=1}^n h(d) S(g,\left\lfloor\frac nd \right\rfloor) \\ = \ & \sum_{d \text{ is PN}}^n h(d) S(g,\left\lfloor\frac nd \right\rfloor) \end{aligned}\]

随后只需要找到 \(h\) 的所有取值即可。这点可以在枚举 \(\text{PN}\) 时完成。为此我们还需要得到每个 \(h(p^c)\)。

有两种方法。第一种是简单地找到 \(h(p^c)\) 仅关于 \(p,c\) 的计算式。而第二种是由 \(f(p^c) = \sum_{i = 0}^c g(p^i) h(p^{c - i})\) 得到 \(h(p^c) = f(p^c) - \sum_{i = 1}^{c} g(p^i) h(p^{c - i})\)。还可以分段打表

构造 \(\text{PN}\) 筛的步骤如下：

1. 绞尽脑汁或者利用科技想出一个好求的 \(g\)。
2. 计算 \(h(p^c)\)。
3. dfs 过程中累计答案。

讨论时间复杂度。

第一部分复杂度是 \(O(1\text{ min})\sim O(\infty)\) 不等。

第二部分默认采用第二种方法，第一种方法易分析。
由于需要找到每个 \(p \le \sqrt n\) 的每个 \(c = O(\log n)\) 位置的点值，因此点值数量是 \(O(\sqrt n)\) 的。每个 \(h(p^c)\) 都需要重新枚举一遍 \(1\sim c\)，因此总时间复杂度为 \(O(\sqrt n \log n)\)。这个上界较松。同时根据 \(h\) 的性质可以再加优化。

第三部分的复杂度依赖于计算 \(S(g, n)\) 的复杂度。如果可以 \(O(1)\) 计算，则这部分的复杂度为 \(O(\sqrt n)\)。而借助杜教筛的总复杂度有 \(O(n^{2/3})\)。

\(\text{PN}\) 筛部分的空间复杂度显然为 \(O(\sqrt n)\)。

这东西被 fjzzq 称作”杜教筛的一个拓展“。确实，从构造到生成都很像。
但是感觉这东西还是在与杜教筛剥离的情况下能得到更优秀的（理论）复杂度。

放一道 command-block 博客里的例题：

设积性函数 \(f\) 满足 \(f(p^c) = p^{c\times k_1} + p^{c\times k_2}\)，求 \(S(f, n)\)。

\(1\le n\le 10^{12}, 1\le k_1, k_2\le 10\)。

考虑素数拟合，\(f(p) = p^{k_1} + p^{k_2}\)，容易发现构造 \(g_1 = \text{id}_{k_1}, g_2 = \text{id}_{k_2}\) 即可得到 \(g = g_1 * g_2\)。
考虑 \(h(p^c)\)，我们就需要得到得到 \(g\) 的前 \(\sqrt n\) 个点值。首先得到 \(g_1, g_2\)，这可以直接线性筛出来，复杂度为 \(\widetilde O(\sqrt n)\)。随后在 \(O(\sqrt n \log n)\) 的复杂度内做狄利克雷卷积可得。
\(S(g, n)\) 可以考虑 \(k\) 的取值较小，\(\text{id}_k\) 任意位置的前缀和都可以 \(O(1)\)，这样可以通过整除分块 \(O(\sqrt n)\) 算得。在取 \(\text{PN}\) 时计算，时间复杂度易证得 \(O(\sqrt n \log n)\)。

因此有总时间复杂度 \(O(\sqrt n\log n)\)。

推荐阅读：一类积性函数的特殊前缀和问题，whx1003

板子就放 \(\text{Min25}\) 板子的代码吧。贺 OIwiki 的。

需要注意的是可能乘爆 ll。Submission.

code

#include <bits/stdc++.h>
#define int long long
using namespace std; using pii = pair<int,int>; using vi = vector<int>; using vp = vector<pii>; using ll = long long; 
using ull = unsigned long long; using db = double; using ld = long double; using lll = __int128_t;
mt19937 rnd(chrono::steady_clock::now().time_since_epoch().count());
template <typename T> T rand(T l, T r) { return uniform_int_distribution<T>(l, r)(rnd); }
template <typename T> void get(T & x) {
	x = 0; char ch = getchar(); bool f = false; while (ch < '0' or ch > '9') f = f or ch == '-', ch = getchar();
	while ('0' <= ch and ch <= '9') x = (x << 1) + (x << 3) + ch - '0', ch = getchar(); f && (x = -x); 
} template <typename T, typename ... Args> void get(T & a, Args & ... b) { get(a); get(b...); }
#define iot ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
#define timer cerr << 1. * clock() / CLOCKS_PER_SEC << '\n';
template <typename T1, typename T2> T1 min(T1 a, T2 b) { return a < b ? a : b; }
template <typename T1, typename T2> T1 max(T1 a, T2 b) { return a > b ? a : b; }
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
#define ufile(x) 
#define rep(i,s,t) for (register int i = (s), i##_ = (t) + 1; i < i##_; ++ i)
#define pre(i,s,t) for (register int i = (s), i##_ = (t) - 1; i > i##_; -- i)
#define Aster(i, s) for (int i = head[s], v; i; i = e[i].next)
#define all(s) s.begin(), s.end()
#define eb emplace_back
#define pb pop_back
#define em emplace
const int N = 2e6 + 10, L = 2e6;
const int inf = 0x3f3f3f3f;
const ll infll = 0x3f3f3f3f3f3f3f3fll;

// int mod;
// const int mod = 10007;
// const int mod = 469762049, g = 3;
// const int mod = 998244353; // const int g = 3;
// const int mod = 1004535809, g = 3;
const int mod = 1e9 + 7;
// const int mod = 1e9 + 9;
// const int mod = 1e9 + 3579, bse = 131;
/* add / sub */ 			template <typename T1, typename T2> T1 add(T1 a, T2 b) { return (a += b) >= mod ? a - mod : a; } template <typename T1, typename ...Args> T1 add(T1 a, Args ... b) { return add(a, add(b...)); } template <typename T1, typename T2> T1 sub(T1 a, T2 b) { return (a -= b) < 0 ? a + mod : a; } template <typename T1, typename ...Args> T1 sub(T1 a, Args ... b) { return sub(a, add(b...)); } template <typename T1, typename T2> void addi(T1 & a, T2 b) { (a += b) >= mod ? (a -= mod) : true; } template <typename T1, typename ...Args> void addi(T1 & a, Args ...b) { addi(a, add(b...)); } template <typename T1, typename T2> void subi(T1 & a, T2 b) { (a -= b) < 0 ? (a += mod) : true; } template <typename T1, typename ...Args> void subi(T1 & a, Args ...b) { subi(a, add(b...)); }
/* Fastmod / mul */ 		struct FastMod { int m; ll b; void init(int _m) { m = _m; if (m == 0) m = 1; b = ((lll)1<<64) / m; } FastMod(int _m) { init(_m); } int operator() (ll a) {ll q = ((lll)a * b) >> 64; a -= q * m; if (a >= m) a -= m; return a; } } Mod(mod); int mul(int a, int b) { return Mod(1ll * a * b); } template <typename T1, typename T2> int mul(T1 a, T2 b) { return Mod((long long)(1ll * a * b)); } template <typename T, typename ...Args> int mul(T a, Args ...b) { return mul(a, mul(b...)); }  template <typename T1, typename T2> void muli(T1 & a, T2 b) { a = Mod(1ll * a * b); } template <typename T1, typename ...Args> void muli(T1 & a, Args ...b) { muli(a, mul(b ...)); } // /* trivial multiple function(idk whether its useful*/ int mul(int a, int b) { return 1ll * a * b % mod; } template <typename T1, typename T2> int mul(T1 a, T2 b) { return (long long)(1ll * a * b) % mod; } template <typename T, typename ...Args> int mul(T a, Args ...b) { return mul(a, mul(b...)); }
/* qp init C */ 			template <typename T1, typename T2> T1 qp(T1 a, T2 b) { T1 ret = 1; for (; b > 0; a = mul(a, a), b >>= 1) if (b & 1) ret = mul(ret, a); return ret; } int fac[N], inv[N], ifc[N]; template<typename T1> void init_fac(T1 bnd, int mask = 0b111) { if ((mask >> 2) & 1) { static int __var_fac_bnd; if (__var_fac_bnd == 0) fac[0] = fac[1] = __var_fac_bnd = 1; rep(i, __var_fac_bnd + 1, bnd) fac[i] = mul(fac[i - 1], i); __var_fac_bnd = bnd; } if ((mask >> 1) & 1) { static int __var_inv_bnd; if (__var_inv_bnd == 0) inv[0] = inv[1] = __var_inv_bnd = 1; rep(i, __var_inv_bnd + 1, bnd) inv[i] = mul(mod - mod / i, inv[mod % i]); __var_inv_bnd = bnd; } if ((mask >> 0) & 1) { static int __var_ifac_bnd, __var_inv_bnd; if (__var_ifac_bnd == 0) ifc[0] = ifc[1] = __var_ifac_bnd = 1; if (__var_inv_bnd < bnd) init_fac(bnd, 0b010); rep(i, __var_ifac_bnd + 1, bnd) ifc[i] = mul(ifc[i - 1], inv[i]); __var_ifac_bnd = bnd;  } } template<typename T1, typename T2> int C(T1 n, T2 m) { if (n < m) return 0; return mul(fac[n], ifc[m], ifc[n - m]); } template<typename T1, typename T2> int invC(T1 n, T2 m) { if (n < m) return 0; return mul(ifc[n], fac[m], fac[n - m]); }
// /* inv 2/3/6 / phi */ 		constexpr int __var_pow(int a, int b) { int ret = 1; for (; b > 0; b >>= 1, a = 1ll * a * a % mod) if (b & 1) ret = 1ll * ret * a % mod; return ret; } constexpr int __var_phi(int x) { if (x <= 2) return 1; int ret = x; for (int i = 2; 1ll * i * i <= x; ++i)  if (x % i == 0) { ret = ret / i * (i - 1); while (x % i == 0) x /= i; } if (x > 1) ret = ret / x * (x - 1); return ret; } const int __phi_mod = __var_phi(mod); const int inv2 = __var_pow(2, __phi_mod - 1), inv3 = __var_pow(3, __phi_mod - 1), inv6 = __var_pow(6, __phi_mod - 1);

ll n;
int ans, inv2, inv6, h[N][40];
bool vis[N][40];
int mp[N];

int g[N], sg[N], prime[N], cnt, phi[N];
void sieve(int n = L) {
	phi[1] = 1;
	rep(i,2,n) {
		if (!phi[i]) phi[i] = i - 1, prime[++ cnt] = i;
		rep(j,1,cnt) {
			int t = i * prime[j];
			if (t > n) break;
			if (i % prime[j] == 0) {
				phi[t] = 1ll * phi[i] * prime[j];
				break;
			} else phi[t] = 1ll * phi[i] * phi[prime[j]];
		}
	} 
	rep(i,1,n) g[i] = mul(i, phi[i]);
	rep(i,1,n) sg[i] = add(sg[i - 1], g[i]);
	rep(i,1,cnt) h[i][0] = 1, h[i][1] = 0, vis[i][0] = vis[i][1] = 1;
	inv2 = (mod + 1) / 2;
	inv6 = (mod + 1) / 6;
	memset(mp, -1, sizeof mp);
}

int S1(ll x) { x = Mod(x); return Mod(x * (x + 1) >> 1); }
int S2(ll x) { x = Mod(x); return mul(x, x + 1, 2 * x + 1, inv6); }

int getg(ll x) {
	if (x <= L) return sg[x];
	if (~mp[n / x]) return mp[n / x];
	int ret = S2(x);
	for (ll l = 2, r, t; l <= x; l = r + 1) {
		t = x / l, r = x / t;
		ret = sub(ret, mul(sub(S1(r), S1(l - 1)), getg(t)));
	} 
	return mp[n / x] = ret;
}

void dfs(ll d, int hd, int pid) {
	ans = add(ans, mul(hd, getg(n / d)));
	rep(i,pid,cnt) {
		if (i > 1 and d * prime[i] * prime[i] > n) break;
		for (ll x = d * prime[i] * prime[i], c = 2; 0 < x and x <= n; x *= prime[i], ++ c) {
			if (!vis[i][c]) {
				int f = qp(prime[i], c), g = mul(prime[i], prime[i] - 1), t = mul(prime[i], prime[i]);
				f = mul(f, sub(f, 1));
				rep(j,1,c) {
					f = sub(f, mul(g, h[i][c - j]));
					g = mul(g, t);
				} h[i][c] = f;
				vis[i][c] = 1;
			}
			if (h[i][c]) dfs(x, mul(hd, h[i][c]), i + 1);
		}
	}
}

signed main() {
	cin >> n;
	sieve();
	dfs(1, 1, 1);
	cout << ans << '\n';
}