題目連結

容易看出，

\[\frac{f(a,b)}{ab} = \frac{f(g,g)}{gg} \]

轉化為維護 \(f(g,g)\)。現在我們要求和：

\[\sum_{i<=K,j<=K}f(i,j) \]

\[=\sum_{g=1}^Kf(g)\sum_{i,j<=K/g}ij[(i,j)=1] \]

然後我們需要一個技巧

\[\large \sum_{i=1}^xi[(i,x)=1]=\frac{x(\varphi(x)+\epsilon(x))}{2} \]

感性地理解一下。顯然，如果 \((x,i)=1\)，那麼必定有 \((x-i,i)=1\)，因此和 \(x\) 互質的數是“對稱的”。不過 \(x=1\)

然後剩下的就簡單了：

\[G(x)=\sum_{i,j<=K/g}ij[(i,j)=1]=\sum_{i=1}^xi^2\varphi(i) \]

這個直接預處理。

\[\sum_{g=1}^Kf(g)G(\left \lfloor \frac{K}{g} \right \rfloor ) \]

這個整除分塊。不過我們需要 \(O(1)\) 算出一個 \(f\) 的字首和。這個分塊維護一下字首和即可。

關鍵程式碼：

int m, n;
ll f[N], sum[N], G[N];
int block, blong[N], st[NN], ed[NN], tag[NN];
int pri[N], pcnt, phi[N];
bool depri[N];
inline void init() {
	for (register int i = 1; i <= n; ++i)	G[i] = (G[i - 1] + 1ll * i * i % P * phi[i]) % P;
	for (register int i = 1; i <= n; ++i)
		f[i] = 1ll * i * i % P, sum[i] = (f[i] + sum[i - 1]) % P;
}
inline void modify(int p, int x) {
	int tmp = x;//Attention!!!!!
	x = x - f[p];
	f[p] = tmp;
	for (register int i = p; i <= ed[blong[p]]; ++i)	ADD(sum[i], x);
	for (register int i = blong[p] + 1; i <= blong[n]; ++i)	ADD(tag[i], x);
}
inline ll get_sum(int p) {
	return sum[p] + tag[blong[p]];
}
inline void query(int K) {
	ll res = 0;
	for (register int l = 1, r; l <= K; l = r + 1) {
		r = (K / (K / l));
		res = (res + (get_sum(r) - get_sum(l - 1)) * G[K / l]) % P;
	}
	printf("%lld\n", (res % P + P) % P);
}
inline void work() {
	for (register int i = 1; i <= m; ++i) {
		int a, b, x, k; read(a), read(b), read(x), read(k);
		int g = get_gcd(a, b);
		x %= P;//Attention!!!!!!!!!
		modify(g, 1ll * x * g % P * g  % P * quickpow(a, P - 2) % P * quickpow(b, P - 2) % P);
		query(k);
	}
}