LOJ561「LibreOJ Round #9」CommonAnts 的调和数

给定长为 \(n\) 的初始为 \(0\) 的序列 \(\{a\}_{i=1}^n\)。

\(m\) 次操作，给定正整数 \(p,x\)，表示令 \(\forall i\le\lfloor\frac np\rfloor,a_{pi}:=a_{pi}+xi\)。

\(q\) 次询问，给定正整数 \(k\)，表示求 \((\sum_{i=1}^{\lfloor n/k\rfloor}i\cdot a_{ki})\bmod 998244353\)。

\(p_i,k_i\le n\le 10^9,m,q\le 2\cdot 10^5,x_i\le 10^9,\omega(\prod_{i=1}^m p_i\prod_{i=1}^qk_i)\le 10\)。

---

你真的会埃氏筛么？（雾

设 \(z\) 是 \(\prod_{i=1}^m p_i\prod_{i=1}^q k_i\) 的所有质因子之积。搜搜可以得到 \(\le n\) 的数中 \(z^{\infty}\) 的因数个数 \(d\approx 2\cdot 10^5\)。

将所有位置的数在其与 \(z^{\infty}\) 的 \(\gcd\) 处统计贡献，系数大概是

\[f(k)=\sum_{i=1}^ki^2[i\perp z] \]

其中 \(k\) 是 \(\lfloor n/i\rfloor\) 的形式，只有 \(2\sqrt n\) 个数，反演一下就可以算出来。

求和的时候类似 Dirichlet 前/后缀和，时间复杂度 \(O((d+\sqrt n)\omega)\)。

#include<bits/stdc++.h>
#define PB emplace_back
#define MP make_pair
#define fi first
#define se second
using namespace std;
typedef long long LL;
typedef pair<int, int> pii;
const int N = 200003, mod = 998244353, inv3 = (mod+1)/3;
template<typename T>
void read(T &x){
	int ch = getchar(); x = 0; bool f = false;
	for(;ch < '0' || ch > '9';ch = getchar()) f |= ch == '-';
	for(;ch >= '0' && ch <= '9';ch = getchar()) x = x * 10 + ch - '0';
	if(f) x = -x;
} void qmo(int &x){x += x >> 31 & mod;}
int n, m, q, pc, cnt, a[N], tot, b[N], pri[10], que[N];
unordered_map<int, int> f, dp;
int SUP(int x){return (x*(x+1ll)>>1)%mod*(2*x+1)%mod*inv3%mod;}
void fact(int x){
	for(int i = 0;i < pc;++ i)
		while(!(x % pri[i])) x /= pri[i];
	for(int i = 2;i * i <= x;++ i) if(!(x % i)){
		pri[pc++] = i; x /= i; while(!(x % i)) x /= i;
	} if(x > 1) pri[pc++] = x;
}
void dfs(int d, LL prd){
	if(d == pc){a[cnt++] = prd; return;}
	while(prd <= n){dfs(d+1, prd); prd *= pri[d];}
}
int main(){
	read(n); read(m); read(q);
	for(int i = 1, p, v;i <= m;++ i){
		read(p); read(v); qmo(dp[p] += v - mod); fact(p);
	} for(int i = 1;i <= q;++ i){read(que[i]); fact(que[i]);}
	dfs(0, 1); sort(a, a + cnt);
	for(int i = 1;i * i <= n;++ i){
		b[tot++] = i; if(i * (i+1) <= n) b[tot++] = n / i;
	} sort(b, b + tot);
	for(int i = 0;i < tot;++ i) f[b[i]] = SUP(b[i]);
	for(int i = 0;i < pc;++ i)
		for(int j = tot-1;~j;-- j)
			qmo(f[b[j]] -= (LL)f[b[j] / pri[i]] * pri[i] % mod * pri[i] % mod);
	for(int i = 0;i < pc;++ i)
		for(int j = 0;j < cnt;++ j) if(!(a[j] % pri[i]))
			dp[a[j]] = (dp[a[j]] + (LL)dp[a[j] / pri[i]] * pri[i]) % mod;
	for(int i = 0;i < cnt;++ i) dp[a[i]] = (LL)dp[a[i]] * f[n / a[i]] % mod;
	for(int i = 0;i < pc;++ i)
		for(int j = cnt-1;~j;-- j) if((LL)a[j] * pri[i] <= n)
			dp[a[j]] = (dp[a[j]] + (LL)dp[a[j] * pri[i]] * pri[i]) % mod;
	for(int i = 1;i <= q;++ i) printf("%d\n", dp[que[i]]);
}