积性函数求和问题的一种筛法

前言

不知不觉鸽了这么久了，马上也要到今年的冬令营了… 还是早点把这个去年的想法写完吧。

这篇文章主要讨论曾经常在 OI 里出现的积性函数求和问题。

基本定义和性质

记号$\lfloor x\rfloor$表示不超过$x$的最大整数，记号$\lceil x\rceil$表示不小于$x$的最小整数。

记号$d|n$表示$d$整除$n$，$\gcd (x,y)$表示$x$和$y$的最大公约数，$[x]$在$x$为真时值为$1$，否则值为$0$。

定义函数$f:{\mathbb{Z}}^{+}\to \mathbb{F}$为数论函数，其中$\mathbb{F}$是一个域，称由全体数论函数构成的集合为$\mathcal{U}$。

若数论函数$g:{\mathbb{Z}}^{+}\to \mathbb{F}$满足$\forall x,y\in {\mathbb{Z}}^{+},\gcd (x,y)=1,\exists g(xy)=g(x)g(y)$，则称$g$为积性函数，称由全体积性函数构成的集合为$\mathcal{V}$。

若数论函数$h:{\mathbb{Z}}^{+}\to \mathbb{F}$满足$\forall x,y\in {\mathbb{Z}}^{+},\exists h(xy)=h(x)h(y)$，则称$h$为完全积性函数。

定义$P$为素数集合，集合中的第$i$个元素是第$i$小的素数$p_i$。对于一个正整数$x=\prod p_i^{e_{x_i}}$，定义该正整数的坐标$e_x=[x{]}_P=(e_{x,1},e_{x,2},e_{x,3},\dots )$。

对$p\in P$，定义 V p = { 1 , p , p 2 , p 3 , … } V_p = \{1, p, p^2, p^3, \ldots\} Vp​={1,p,p2,p3,…}，表示$p$的非负次幂构成的集合。

定义$f_p(n)=f(n)[n\in V_p]$，即保留$f$中下标在$V_p$中的项后的数论函数。

定义$\mathcal{U}$上的 Dirichlet 卷积运算$\ast$，满足$(f\ast g)(n)={\sum }_{d|n}f(d)g\left(\frac{n}{d}\right)$，用$\hat{\prod }$表示 Dirichlet 卷积的连续运算，Dirichlet 卷积有以下性质：

性质一：$f\in \mathcal{V}⟺f={\hat{\prod }}_{p\in P}{f}_p$，其中$f_p$只有下标在$V_p$中的位置有值。

证明： 若$f\in \mathcal{V}$，比较等式两端系数可得$f={\hat{\prod }}_{p\in P}{f}_p$；若$f={\hat{\prod }}_{p\in P}{f}_p$，容易验证$f\in \mathcal{V}$。

性质二：$\forall f,g\in \mathcal{V},f\ast g\in \mathcal{V}$

证明：$f\ast g={\hat{\prod }}_{p\in P}{f}_p{g}_p$，由性质一有$f\ast g\in \mathcal{V}$。

性质三：$\forall f\in \mathcal{V},\exists g=f^{-1}\mathrm{s}.\mathrm{t}.f\ast g=ϵ$，其中$ϵ(x)=[x=1]$

定义$n\in {\mathbb{Z}}^{+}$的整除集合为$D_n=\left\{\lfloor \frac{n}{m}\rfloor |m\in {\mathbb{Z}}^{+}\right\}$，整除集合有如下性质。

性质一：$D_n=\left\{0,1,2,\dots ,\lfloor \sqrt{n}\rfloor ,\lfloor \frac{n}{\lfloor \sqrt{n}\rfloor }\rfloor ,\lfloor \frac{n}{\lfloor \sqrt{n}\rfloor -1}\rfloor ,\dots ,n\right\}$

证明： 当$m\le \sqrt{n}$时，$\lfloor n/m\rfloor$的取值已经被包含在集合中，否则$\lfloor \frac{n}{m}\rfloor \le \lfloor \sqrt{n}\rfloor$。而此时$m>\sqrt{n}$，若$\lfloor \frac{n}{m}\rfloor \ne \lfloor \frac{n}{m+1}\rfloor$，则有$\frac{n}{m+1}\le \lfloor \frac{n}{m}\rfloor \le \frac{n}{m}$，但$\frac{n}{m}-\frac{n}{m+1}<1$，因此$\lfloor \frac{n}{m}\rfloor =\lfloor \frac{n}{m+1}\rfloor +1$，于是结论成立。

性质二：$D\left(\lfloor n/m\rfloor \right)\subseteq D(n),n,m\in {\mathbb{Z}}^{+}$

证明： 令$n=kab+rb+c,r\in [0,a)\cap \mathbb{Z},c\in [0,b)\cap \mathbb{Z}$，则$\lfloor n/ab\rfloor =k=\lfloor (ka+r)/a\rfloor =\lfloor \lfloor n/b\rfloor /a\rfloor$，取$b=m$可得结论成立。

性质三： 令$U_n(x)={\min }_{y\in D_n,y\ge x}y$，其中$1\le x\le n$，则$\lfloor n/U_n(x)\rfloor =\lfloor n/x\rfloor$

证明： 不妨设$y=\lfloor n/x\rfloor$，那么$xy\le n$，于是$x\le n/y$，因此这样的$x$的上界为$\lfloor n/y\rfloor \in D_n$，因此$\lfloor n/U_n(x)\rfloor \le \lfloor n/x\rfloor$，而另一个方向是显然的，于是得证。

推论：$\forall x,y\in [1,n]\cap \mathbb{Z},xy\le n⟺U_n(x)U_n(y)\le n$

对$f\in \mathcal{U},n\in {\mathbb{Z}}^{+}$定义$f$定义在$n$的整除集合上的前缀和函数为$S_{n,f}(x)={\sum }_{y\le U_n(x)}f(y),1\le x\le n$，那么在存储$S_{n,f}$时我们只要存储$D_n$中所有下标的值即可。

定理： 给定$S_{n,f}$和$S_{n,g}$，可以在$O(n^{2/3}{\log }^{1/3}n)$的时间内求得$S_{n,f\ast g}$。

证明： 即对$m\in D_n$求$S_{n,f\ast g}(m)={\sum }_{x\in D_m}\left(S_{n,f}(\lfloor m/x\rfloor )-S_{n,f}(\lfloor m/U_m(x+1)\rfloor )\right)S_{n,g}(x)$，由此可以在 $O(\sqrt{m})$ 的时间内对特定的$m$求得该值。另一方面，满足$x,y\in D_n,xy\le S$的$(x,y)$对数是$O(S\log S)$的，因此由性质二和性质三的推论，对所有$m\le S$可以在$O(S\log S)$内求得答案。而对超过$S$的$m$消耗的总时间为$O({\int }_1^{n/S}(n/x{)}^{1/2})=O(n/S^{1/2})$，取$S=n^{2/3}{\log }^{-2/3}n$时得到$O(n^{2/3}{\log }^{1/3}n)$的时间复杂度，故命题正确。

另外注意到，若$f$和$g$是积性的，我们可以在$O(S)$的时间内求得前$S$项的值，于是复杂度可以降到$O(n^{2/3})$。

特别地，如果$S_{n,f}(x)=c$对所有$x\le \sqrt{n}$成立，合并可以用暴力在$O(\sqrt{n}\log n)$的时间内完成。

定义 Powerful Number 为所有的满足$\forall k\in {\mathbb{Z}}^{+},e_{x,k}\ne 1$的$x$，令$W$为全体 Powerful Number 组成的集合。

性质：$n$以内的 Powerful Number 的数目是$\Theta (\sqrt{n})$的。

证明： 一方面，所有的平方数均 powerful，因此数目是$\Omega (\sqrt{n})$的。另一方面，所有的 powerful number 总能表示为$a^2{b}^3$的形式，这样的数的数目为$O({\int }_1^{n^{1/3}}(n/x^3{)}^{1/2}\mathrm{d}x)+O(n^{1/3})=O(n^{1/2})$，因此原命题得证。

定理： 若$f$中只有 Powerful Number 处有值，且我们有$S_{n,g}$，我们可以在$\tilde{O}(n^{0.6})$的时间内求得$S_{n,f\ast g}$。

证明： 由 Powerful Number 的性质，按照$m^{2/3}$分段可知，$S_{m,f}$中值不同的项的数目是$O(m^{1/3})$的。同样采用类似的方法，平衡较小部分的复杂度$\tilde{O}(S)$和较大部分的复杂度$O(n/S^{2/3})$可在$S=\tilde{O}(n^{0.6})$时取得平衡。于是命题成立。

一种筛法

OI 中出现的积性函数求和问题往往分为两种：给定$f_p$的生成规则，求$S_{n,f}(n)$或是求出$S_{n,f}$。

由于对于$\sqrt{n}<p\le n$的$p$是$O(n{\log }^{-1}n)$的，为了给出生成规则我们一般对于这一部分的$p$一定有$f(p)=y(p)$，其中$y(x)={\sum }_{i\in T}{a}_i{x}^i,T\subset {\mathbb{Z}}^{+}$，并且$|T|$是一个常数，且其中元素的大小不太大。

只求单一的$S_{n,f}(n)$的值可以用 Min-25 筛，虽然是$\Omega (n^c)\forall c<1$的但在现在的范围$(\sim 10^{11})$里还是非常 work 的，求$S_{n,f}$的传统方法有：

1. 杜教筛：已知$S_{n,f}$和$S_{n,g}$，可以在$O(n^{2/3})$内求$S_{n,f\ast g^{-1}}$。
2. 洲阁筛（扩展埃筛）：若$|S|$是常数，可以在$O(n^{3/4}{\log }^{-1}n)$内求$S_{n,f}$。

以上资料可以在《2016/2018年候选队论文集》以及 Min-25的博客 中找到。在博客中还有另一个比本文方法复杂度稍劣的方法，并且和本文方法一样没有很大的效率提升。

下面我们将介绍一种$O(n^{2/3}{\log }^{-1}n)$的求$S_{n,f}$的方法。设$f={\hat{\prod }}_{p\in P}{f}_p$，并定义$f_{\mathcal{D}}={\hat{\prod }}_{p\in \mathcal{D}}{f}_p$，那么：

定理一： 给定$f_p$，可以在$O(\sqrt{n}{\log }_{p}n)$的时间内由$S_{n,g}$推得$S_{n,g\ast f_p}$。

证明： 暴力枚举$f_p$中的每一项并更新即可。

因此，若我们对于$p\le n^{1/\alpha }$使用该方法更新，其复杂度为$O(\alpha n^{1/\alpha +1/2})$。

接下来我们将$f=\hat{\prod }{f}_p$按照$p$分成三个部分：${\mathcal{D}}_1=P\cap [1,n^{1/\alpha }]$,${\mathcal{D}}_2=P\cap (n^{1/\alpha },n^{1/2}]$,${\mathcal{D}}_3=P\cap (n^{1/2},n]$。若我们求出了$S_{n,f_{{\mathcal{D}}_2}}$和$S_{n,f_{{\mathcal{D}}_3}}$，由之前的性质，我们可以在$O(\alpha n^{1/\alpha +1/2})$的时间内将其合并为$S_{n,f}$。首先考虑如何求$S_{n,f_{{\mathcal{D}}_2}}$。

定理二：$\forall 1\le k\le {\log }_{2}n$，恰含$k$个素因子的素因子大小至少为$n^{1/\alpha }$的数的个数是$O((\alpha +1{)}^{k}n{\log }^{-1}n)$的。

证明： 当$k=1$时该式成立。若当$k=m$时该式成立，考虑$k=m+1$时这样的数的个数显然为
$$O\left((\alpha +1{)}^m{\int }_{n^{1/\alpha }}^n\frac{\mathrm{d}x}{\ln x}\frac{n}{x\ln x}+(\alpha +1{)}^m{\int }_{n^{1/\alpha }}^n\frac{\mathrm{d}x}{\ln x}\right)=O((\alpha +1{)}^{m+1}n{\log }^{-1}n)$$
传统的计算$S_{n,f_{{\mathcal{D}}_2}}$的方法总要考虑素因子的顺序，这导致了状态的两个维度的产生。即便每个状态的转移已经最优，复杂度也难以优化。因此我们尝试减少这一部分的状态数：如果我们不考虑素因子的顺序会发生什么呢？也就是说，令$g_k(x)=f_{{\mathcal{D}}_2}(x)[\sum e_{x,i}=k]$，我们希望通过计算$g_{k-1}$与${\sum }_{p\in {\mathcal{D}}_2}{f}_p$的卷积来推得$S_{n,g_k}$。

定理三： 若$\alpha$为常数，则计算一次这样的卷积的时间复杂度为$O(n^{2/3}{\log }^{-1}n)$。

证明： 由于只在素数及其幂次处有值，可以注意到$S_{n,{\sum }_{p\in {\mathcal{D}}_2}{f}_p}$中不同的值的期望个数是$O(n^{1/2}{\log }^{-1}n)$的。因次在计算$S_{n,g_{k-1}\ast {\sum }_{p\in {\mathcal{D}}_2}{f}_p}$时，对一个$m$单次枚举的复杂度可以做到$O(m^{1/2}{\log }^{-1}m)$：我们仍然按照 $\sqrt{m}$ 分段，较大部分时我们的素数个数是除了一个 log 的，而较小时值的个数也有限。而由定理二，若$\alpha$为常数，则计算较小部分的复杂度也可以做到$O(S{\log }^{-1}n)$，于是取$S=O(n^{2/3})$可以达到定理中的时间复杂度。

由定理三，当$\alpha$为常数时，我们总的计算时间复杂度为$O(n^{2/3}{\log }^{-1}n)$。

但我们究竟算了什么呢？容易发现，对一个数$x$来说，其对答案的贡献为$(\sum e_{x,i})!\prod f_{{\mathcal{D}}_2}(p_i{)}^{e_{x,i}}(e_{x,i}){!}^{-1}$，因此若$f_{{\mathcal{D}}_2}(p_i^e)=f_{{\mathcal{D}}_2}(p_i{)}^e(e!{)}^{-1}$，我们算出的值恰好就是$S_{n,g_k}$的$k!$倍，只要在最后除去即可。

注意到若我们构造另一个积性函数$f_0$满足$f_0(p_i^e)=f_{{\mathcal{D}}_2}(p_i{)}^e(e!{)}^{-1}$，那么$f_{{\mathcal{D}}_2}\ast f_0^{-1}$不为$0$的位置一定是 Powerful Number 集合$W$的子集。于是我们可以先求出$f_0$的前缀和，再和一个 Powerful 的函数做卷积，这样我们就在$O(n^{2/3}{\log }^{-1}n)$的时间内求出了$S_{n,f_{{\mathcal{D}}_2}}$。

那么如何求$S_{n,f_{{\mathcal{D}}_3}}$呢，我们将问题转化为$|T|$个求所有最小素因子超过$\sqrt{n}$的数的$k$次方之和的问题，令$u_k(x)=x^k[x\in P][x>\sqrt{n}]$，即求$S_{n,u_k}$。注意到超过$\sqrt{n}$的因子至多只有一个，令$v_k(x)=x^k[\forall p_i>\sqrt{n},\exists e_{x,i}=0]$，则有：$S_{n,u_k\ast (v_k+ϵ)}=S_{n,i{d}^k}$，而$S_{n,u_k}$也只有$\sqrt{n}$以上的项有值，因此有$S_{n,v_k+ϵ}$和$S_{n,i{d}^k}$后可以用一个简单的容斥做到$O(\sqrt{n}\log n)$求$S_{n,u_k}$。我们可以利用之前的方法求出$S_{n,v_k}$，这样我们这一部分的时间复杂度也是$O(n^{2/3}{\log }^{-1}n)$。

不妨取$\alpha =6$，总的时间复杂度即为$O(n^{2/3}{\log }^{-1}n)$。

下面是一份示例代码，不保证没有错误、也希望有人能写一个常数更小一点的版本。

const int mod = 1e9 + 7;
inline void inc(int &a, int b){ (a += b) < mod ?: a -= mod; }
inline void inc(int &a, i64 b){ a = (a + b) % mod; }
inline void dec(int &a, int b){ (a -= b) >= 0 ?: a += mod; }
inline void dec(int &a, i64 b){ inc(a, mod - b % mod); }
int fpow(int a, int t){
	static int r;
	for (r = 1; t; t >>= 1, a = (i64)a * a % mod) if (t & 1) r = (i64)r * a % mod;
	return r;
}
const int unit = 7;
i64 x;
int B, S, RS, size, count, fac[unit], ifac[unit], inv[unit];
vector<int> primes, pi;
struct array_t {
	int f[unit];
	array_t(){ memset(f, 0, unit * 4); }
	int &operator[](int i){ return f[i]; }
	int to_val(){
		int sum = 0;
		for (int i = 0; i < unit; ++i) inc(sum, (i64)f[i] * ifac[i]);
		return sum;
	}
} ;
inline void inc(array_t &a, array_t &b){
	for(int i = 0; i < unit; ++i) if (b[i]) inc(a[i], b[i]);
}
inline void inc(array_t &a, array_t &b, int v){
	for(int i = 0; i + 1 < unit; ++i) if (b[i]) inc(a[i + 1], (i64)b[i] * v);
}
template<class value_t>
struct data_t {
	vector<value_t> lsum, ssum;
	data_t(){ lsum.resize(size), ssum.resize(size); }
	void update(int p, vector<int> g){
		int i, j, e, l = S / p;
		for (i = 1; i <= S; ++i) dec(ssum[i], lsum[S]);
		for (i = S; i >= 2; --i) dec(lsum[i], lsum[i - 1]);
		for (i = 1; i <= S; ++i) for (j = i, e = 0; j <= l; inc(ssum[i], (i64)ssum[j *= p] * g[++e]));
		for (i = 2; i <= S; ++i) dec(ssum[i - 1], ssum[i]);
		for (i = S; i >= 1; --i) if (lsum[i]) {
			for (j = i, e = 0; j <= l; inc(lsum[j *= p], (i64)lsum[i] * g[++e]));
			for (j = x / p / j; j; j /= p) inc(ssum[j], (i64)lsum[i] * g[++e]);
		}
		for (i = 2; i <= S; ++i) inc(lsum[i], lsum[i - 1]);
		inc(ssum[S], lsum[S]);
		for (i = S; i >= 2; --i) inc(ssum[i - 1], ssum[i]);
	}
	void gsum(){
		int i;
		for (i = 2; i <= S; ++i) inc(lsum[i], lsum[i - 1]);
		inc(ssum[S], lsum[S]); if (S == RS) lsum[S] = ssum[S];
		for (i = S; i >= 2; --i) inc(ssum[i - 1], ssum[i]);
	}
	void gdif(){
		int i;
		for (i = 2; i <= S; ++i) dec(ssum[i - 1], ssum[i]);
		dec(ssum[S], lsum[S]);
		for (i = S; i >= 2; --i) dec(lsum[i], lsum[i - 1]);
	}
	void add(i64 a, value_t b){ inc(a <= S ? lsum[a] : ssum[x / a], b); }
	void print(){
		debug("function address %llu: \n", u64(this));
		for (int i = 1; i <= S; ++i) debug("S(%d) = %d\n", i, lsum[i]);
		for (int i = S; i >= 1; --i) debug("S(%llu) = %d\n", x / i, ssum[i]);
	}
} ;
void init(){
	int i, j;
	S = sqrtl(x), RS = x / S, size = S + 1, B = powl(x, 1.0 / unit) + 1e-7, pi.resize(size);
	vector<char> sieve(size);
	for (fac[0] = i = 1; i < unit; ++i) fac[i] = (i64)fac[i - 1] * i;
	for (inv[1] = 1, i = 2; i < unit; ++i) inv[i] = mod - (i64)(mod / i) * inv[mod % i] % mod;
	for (ifac[0] = i = 1; i < unit; ++i) ifac[i] = (i64)ifac[i - 1] * inv[i] % mod;
	for (i = 4; i <= S; i += 2) sieve[i] = true;
	for (i = 9; i <= S; i += 6) sieve[i] = true;
	for (i = 5; i * i <= S; i += 2) if (!sieve[i]) {
		for (j = i * i; j <= S; j += 6 * i) sieve[j] = true;
		for (j = i * (i + 6 - 2 * (i % 3)); j <= S; j += 6 * i) sieve[j] = true;
	}
	primes.pb(0);
	for (i = 2; i <= S; ++i) {
		pi[i] = pi[i - 1] + !sieve[i];
		if (!sieve[i]) primes.pb(i), ++count;
	}
	primes.pb(size);
}
data_t<int> solve(vector<vector<int>> f){
	vector<vector<int>> g(count + 1), h(count + 1);
	int i, j, k, l, v, m;
	for (i = 1; i <= count; ++i) {
		m = f[i].size(), g[i].resize(m), h[i].resize(m);
		if (i > pi[B]) {
			for (g[i][0] = j = 1; j < m; ++j) g[i][j] = (i64)g[i][j - 1] * f[i][1] % mod;
			for (j = 0; j < m; ++j) g[i][j] = (i64)g[i][j] * ifac[j] % mod;
		} else g[i] = f[i];
		for (j = 0; j < m; ++j) for (h[i][j] = f[i][j], k = 0; k < j; ++k) dec(h[i][j], (i64)h[i][k] * g[i][j - k]);
	}
	data_t<int> F, G;
	vector<int> fval(count + 1);
	for (i = 1; i <= count; ++i) fval[i] = g[i][1];
	for (i = 1; i <= count; ++i) inc(fval[i], fval[i - 1]);
	{
		data_t<array_t> ts;
		i64 block = std::max((i64)S, std::min(x, i64(powl(x, 2.0 / 3) * 9)));
		block = x / (x / block);
		vector<i64> nxt(size);
		function<void(int, i64, int, int)> dfs = [&](int pos, i64 cur, int cnt, int now){
			int i, u, l, r;
			for (i = pos; cur * primes[i] * primes[i] <= block; ++i) {
				i64 t = cur, a;
				int e = 0, b;
				while (t * primes[i] * primes[i] <= block) {
					t *= primes[i], ++e;
					dfs(i + 1, t, cnt + e, (i64)now * g[i][e] % mod);
					a = t * primes[i], b = (i64)now * g[i][e + 1] % mod * fac[cnt + e + 1] % mod;
					if (a <= S) nxt[a] = a, inc(ts.lsum[a][cnt + e + 1], b);
					else inc(ts.ssum[x / a][cnt + e + 1], b);
				}
			}
			now = (i64)now * fac[cnt + 1] % mod;
			for (i = pos; cur * primes[i] <= S; ++i) {
				u = cur * primes[i], nxt[u] = u;
				inc(ts.lsum[u][cnt + 1], (i64)now * (fval[i] + mod - fval[i - 1]));
			}
			i64 m = x / cur;
			l = i, r = count;
			if (cur * S > block) r = pi[block / cur];
			for (u = l; u <= r; ++u) inc(ts.ssum[m / primes[u]][cnt + 1], (i64)now * (fval[u] + mod - fval[u - 1]));
		} ;
		dfs(pi[B] + 1, 1, 0, 1);
		++ts.lsum[nxt[1] = 1][0];
		ts.gsum();
		for (i = 1; i <= S; ++i) if (!nxt[i]) nxt[i] = nxt[i - 1];
		l = x / block - 1;
		for (i = l; i >= 1; --i) {
			i64 m = x / i;
			array_t sum;
			sum[0] = 1;
			for (j = pi[B] + 1; (i64)i * primes[j] <= S && j <= count; ++j) inc(sum, ts.ssum[i * primes[j]], fval[j] + mod - fval[j - 1]);
			for (; (i64)primes[j] * primes[j] <= m && j <= count; ++j) inc(sum, ts.lsum[m / primes[j]], fval[j] + mod - fval[j - 1]);
			for (; j <= count; j = v + 1) {
				k = nxt[m / primes[j]];
				v = (i64)k * S > m ? pi[m / k] : count;
				inc(sum, ts.lsum[k], fval[v] + mod - fval[j - 1]);
			}
			ts.ssum[i] = sum;
		}
		for (i = 1; i <= S; ++i) G.lsum[i] = ts.lsum[i].to_val(), G.ssum[i] = ts.ssum[i].to_val();
	}
	for (i = pi[B]; i >= 1; --i) G.update(primes[i], g[i]);
	{
		vector<pair<i64, int>> vh;
		function<void(int, i64, int)> search = [&](int pos, i64 rem, int now){
			vh.emplace_back(x / rem, now);
			while ((i64)primes[pos] * primes[pos] <= rem) {
				i64 t = rem / primes[pos] / primes[pos];
				int e = 2;
				for (; t; t /= primes[pos], ++e) if (h[pos][e]) search(pos + 1, t, (i64)now * h[pos][e] % mod);
				++pos;
			}
		} ;
		search(1, x, 1);
		sort(vh.begin(), vh.end());
		for (i = 0; i + 1 < vh.size(); ++i) if (vh[i].first == vh[i + 1].first) inc(vh[i + 1].second, vh[i].second), vh[i].second = 0;
		vector<pair<i64, int>> vv;
		for (i = 0; i < vh.size(); ++i) if (vh[i].second) vv.pb(vh[i]);
		vh.swap(vv), vv.clear();
		vh.emplace_back(x + 1, 0);
		vector<int> hsum(RS + 1), hcnt(size), nxt(size);
		for (auto it : vh) inc(it.first > S ? hsum[x / it.first] : hcnt[it.first], it.second);
		for (i = RS - 1; i >= 1; --i) inc(hsum[i], hsum[i + 1]);
		for (i = 1; i <= S; ++i) inc(hcnt[i], hcnt[i - 1]); 
		nxt[S] = S;
		for (i = S - 1; i >= 1; --i) nxt[i] = hcnt[i] == hcnt[i + 1] ? nxt[i + 1] : i;
		i64 low = pow(x, 0.6);
		low = x / (x / low);
		for (j = 0; vh[j].first <= low; ++j) {
			for (i = 1; i * vh[j].first <= S; ++i) inc(F.lsum[i * vh[j].first], (i64)vh[j].second * (G.lsum[i] + mod - G.lsum[i - 1]));
			for (i64 it = x / vh[j].first; i * vh[j].first <= low && i <= S; ++i) inc(F.ssum[it / i], (i64)vh[j].second * (G.lsum[i] + mod - G.lsum[i - 1]));
		}
		F.gsum(), low = std::min(x / low, S + 1LL);
		for (i = S; i >= low; --i) for (j = 0; i * vh[j].first <= S; ++j) inc(F.ssum[i], (i64)vh[j].second * (G.ssum[i * vh[j].first] + mod - G.lsum[S]));
		for (i = low - 1; i >= 1; --i) {
			int w = 0;
			i64 m = x / i;
			for (j = 0; i * vh[j].first <= S; ++j) inc(w, (i64)vh[j].second * G.ssum[i * vh[j].first]);
			for (; vh[j].first <= S; ++j) inc(w, (i64)vh[j].second * G.lsum[m / vh[j].first]);
			for (j = 1; i * (j + 1) <= RS && j <= S; ++j) inc(w, (i64)G.lsum[j] * (hsum[i * j] + mod - hsum[i * (j + 1)]));
			if (j <= S) inc(w, (i64)G.lsum[j] * hsum[i * j]);
			F.ssum[i] = w;
		}
	}
	return F;
}
struct lagrange{
	vector<int> init, coef;
	int n;
	lagrange(vector<int> init): init(init){
		int i;
		coef.resize(n = init.size());
		vector<int> ifac(n), inv(n);
		for (ifac[0] = i = 1; i < n; ++i) {
			inv[i] = i > 1 ? mod - (i64)inv[mod % i] * (mod / i) % mod : 1;
			ifac[i] = (i64)ifac[i - 1] * inv[i] % mod;
		}
		for (i = 0; i < n; ++i) {
			coef[i] = (i64)init[i] * ifac[n - i - 1] % mod * ifac[i] % mod;
			if ((n & 1) == (i & 1)) coef[i] = mod - coef[i];
		}
	}
	int get(i64 w){
		int x = w % mod; 
		if (x < n) return init[x];
		int res = 0, a = 1, i;
		for (i = 0; i < n; ++i) {
			res = ((i64)res * (x - i) + (i64)a * coef[i]) % mod;
			a = (i64)a * (x - i) % mod;
		}
		return (res + mod) % mod;
	}
};
vector<vector<int>> generate(function<int(int, int)> func){
	int i, e;
	i64 t;
	vector<vector<int>> res(count + 1);
	for (i = 1; i <= count; ++i) for (res[i].pb(t = e = 1); t * primes[i] <= x; t *= primes[i], ++e) res[i].pb(func(primes[i], e));
	return res;
}
vector<vector<int>> conv(vector<vector<int>> a, vector<vector<int>> b){
	int i, j, k, l;
	decltype(a) c(count + 1);
	for (i = 1; i <= count; ++i) for (c[i].resize(l = a[i].size()), j = 0; j < l; ++j) for (k = 0; j + k < l; ++k) inc(c[i][j + k], (i64)a[i][j] * b[i][k]);
	return c;
}
data_t<int> primesum(int s){
	int i, j, e, L = std::max(s + 2, S);
	i64 t;
	vector<int> coef(L + 1);
	for (i = 1; i <= L; ++i) coef[i] = 1;
	for (i = 2; i <= L; ++i) if (coef[i] == 1) for (coef[i] = fpow(i, s), j = 2; i * j <= L; ++j) coef[i * j] = (i64)coef[i] * coef[j] % mod;
	vector<int> psum(coef.begin(), coef.begin() + s + 3);
	for (i = 1; i < psum.size(); ++i) inc(psum[i], psum[i - 1]);
	lagrange sol(psum);
	data_t<int> fw = solve(generate([=](int p, int e){ return fpow(p, e * s); }));
	for (i = 1; i <= S; ++i) fw.ssum[i] = (sol.get(x / i) + mod - fw.ssum[i]) % mod, fw.lsum[i] = 0;
	for (i = S; i >= 1; --i) for (j = 2; i * j <= S; ++j) dec(fw.ssum[i], (i64)coef[j] * fw.ssum[i * j]);
	return fw;
};
template<class value_t>
data_t<value_t> work(data_t<value_t> a, const function<int(int)> &op){
	for (int i = 1; i <= S; ++i) a.ssum[i] = op(a.ssum[i]), a.lsum[i] = op(a.lsum[i]);
	return a;
}
template<class value_t>
data_t<value_t> work(data_t<value_t> a, data_t<value_t> b, const function<int(int, int)> &op){
	data_t<value_t> c;
	for (int i = 1; i <= S; ++i) {
		c.ssum[i] = op(a.ssum[i], b.ssum[i]);
		c.lsum[i] = op(a.lsum[i], b.lsum[i]);
	}
	return c;
}
data_t<int> polysum(vector<int> poly){
	int n = poly.size(), i, j;
	vector<int> num(n + 2);
	for (i = 1; i <= n + 1; ++i) {
		for (j = n - 1; j >= 0; --j) num[i] = ((i64)num[i] * i + poly[j]) % mod;
		num[i] = (num[i] + num[i - 1]) % mod;
	}
	lagrange sol(num);
	data_t<int> res;
	for (i = 1; i <= S; ++i) res.lsum[i] = sol.get(i), res.ssum[i] = sol.get(x / i);
	return res;
}