容易看出，
[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) 比较特殊，它等于1.
然后剩下的就简单了：
[G(x)=sum_{i,j<=K/g}ij[(i,j)=1]=sum_{i=1}^xi^2varphi(i)
]
这个直接预处理。
[sum_{g=1}^Kf(g)G(left lfloor frac{K}{g}
ight
floor )
]
这个整除分块。不过我们需要 (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
", (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);
}
}