题目链接
题解
根据题中的式子，手玩一下发现和(gcd)很像
化一下式子：
[egin{aligned}
bf(a,a + b) &= (a + b)f(a,b) \
frac{f(a,a + b)}{a + b} &= frac{f(a,b)}{b} \
frac{f(a,a + b)}{a(a + b)} &= frac{f(a,b)}{ab} \
frac{f(a,b)}{ab} &= frac{f(d,d)}{d^2} \
end{aligned}
]
其中(d = gcd(a,b))
那么我们有
[egin{aligned}
ans &= sumlimits_{i = 1}^{k} sumlimits_{j = 1}^{k} f(i,j) \
&= sumlimits_{d = 1}^{k} f(d,d) sumlimits_{d|i} sumlimits_{d|j} frac{ij}{d^2} quad [gcd(i,j) == d] \
&= sumlimits_{d = 1}^{k} f(d,d) sumlimits_{i = 1}^{lfloor frac{k}{d}
floor} sumlimits_{j = 1}^{lfloor frac{k}{d}
floor} ij quad [i perp j] \
end{aligned}
]
令
[g(n) = sumlimits_{i = 1}^{n} sumlimits_{j = 1}^{n} ij quad [i perp j]
]
由于
[sumlimits_{i = 1}^{n} i quad [i perp n] = frac{nvarphi(n)}{2}
]
所以
[egin{aligned}
g(n) &= sumlimits_{i = 1}^{n} i imes 2 imes frac{ivarphi(i)}{2} \
&= sumlimits_{i = 1}^{n} i^2varphi(i) \
end{aligned}
]
我们可以线性筛(O(n))预处理出(g(n))
对于答案的式子，可以(O(sqrt{k}))整除分块
所以我们只需要(O(1))计算(f(d,d))的前缀和
分块即可
块外维护块的前缀和，块内维护块内前缀和
这样修改是(O(sqrt{n}))的，修改复杂度(O(msqrt{n}))
且询问时(O(1))的
总复杂度(O(msqrt{n}))
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<map>
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define mp(a,b) make_pair<int,int>(a,b)
#define cls(s) memset(s,0,sizeof(s))
#define cp pair<int,int>
#define LL long long int
using namespace std;
const int maxn = 4000005,maxm = 400005,INF = 1000000000,P = 1000000007;
inline LL read(){
LL out = 0,flag = 1; char c = getchar();
while (c < 48 || c > 57){if (c == '-') flag = -1; c = getchar();}
while (c >= 48 && c <= 57){out = (out << 3) + (out << 1) + c - 48; c = getchar();}
return out * flag;
}
int m,n,sum[maxn],S[maxn],b[maxn],L[maxn],R[maxn],bi,B;
inline void modify(int u,int x){
int v = x % P;
if (u == L[b[u]]) v = ((v - S[u]) % P + P) % P;
else v = ((v - (S[u] - S[u - 1]) % P) % P + P) % P;
v = (v + P) % P;
for (int i = u; i <= R[b[u]]; i++)
S[i] = (S[i] + v) % P;
for (int i = b[u]; i <= bi; i++)
sum[i] = (sum[i] + v) % P;
}
inline int query(int u){
if (!u) return 0;
return (S[u] + sum[b[u] - 1]) % P;
}
int p[maxn],pi,isn[maxn],phi[maxn],g[maxn],val[maxn];
void init(){
phi[1] = 1;
for (register int i = 2; i <= n; i++){
if (!isn[i]) p[++pi] = i,phi[i] = i - 1;
for (int j = 1; j <= pi && i * p[j] <= n; j++){
isn[i * p[j]] = true;
if (i % p[j] == 0){
phi[i * p[j]] = phi[i] * p[j];
break;
}
phi[i * p[j]] = phi[i] * (p[j] - 1);
}
}
for (register int i = 1; i <= n; i++)
g[i] = (g[i - 1] + 1ll * val[i] * phi[i] % P) % P;
}
LL gcd(LL a,LL b){return b ? gcd(b,a % b) : a;}
int main(){
m = read(); n = read(); B = (int)sqrt(n) + 1;
for (register int i = 1; i <= n; i++){
b[i] = i / B + 1; val[i] = 1ll * i * i % P;
if (b[i] != b[i - 1]) R[b[i - 1]] = i - 1,L[b[i]] = i;
sum[b[i]] = (sum[b[i]] + val[i]) % P;
if (i != L[b[i]]) S[i] = S[i - 1];
S[i] = (S[i] + val[i]) % P;
}
R[b[n]] = n; bi = b[n];
for (register int i = 1; i <= bi; i++)
sum[i] = (sum[i] + sum[i - 1]) % P;
init();
LL a,b,x,k,d,ans;
while (m--){
a = read(); b = read(); x = read(); k = read();
d = gcd(a,b); x /= (a / d) * (b / d); x %= P;
modify(d,x);
ans = 0;
for (int i = 1,nxt; i <= k; i = nxt + 1){
nxt = k / (k / i);
ans = (ans + 1ll * (query(nxt) - query(i - 1)) % P * g[k / i] % P) % P;
}
printf("%lld
",(ans % P + P) % P);
}
return 0;
}