【题解】CQOI2017 - 小 Q 的表格
https://www.luogu.com.cn/problem/P3700
首先考虑题面给的两个式子。由式二可以得到：
\[\dfrac{f(a,a+b)}{a(a+b)}=\dfrac{f(a,b)}{ab} \]发现这个很像辗转相除，可得
\[\dfrac{f(a,b)}{ab}=\dfrac{f(a,a\bmod b)}{a(a\bmod b)} \]然后由式一转换，最终可得
\[\dfrac{f(a,b)}{ab}=\dfrac{f(\gcd(a,b),\gcd(a,b))}{\gcd(a,b)^2} \]设 \(F(d)\) 表示 \(f(d,d)\)，则每次修改 \((a,b,k)\) 就会使 \(F(d)\) 变为 \(\dfrac{k\gcd(a,b)^2}{ab}\)。
再来用 \(F\) 表示答案：
\[\begin{aligned} \sum_{i=1}^n\sum_{j=1}^nf(i,j)&=\sum_{k=1}^nF(k)\sum_{i=1}^n\sum_{j=1}^n\dfrac{ij}{\gcd(i,j)^2}[\gcd(i,j)=k]\\ &=\sum_{k=1}^nF(k)\sum_{i=1}^{\frac nk}\sum_{j=1}^{\frac nk}ij[\gcd(i,j)=1]\\ &=\sum_{k=1}^nF(k)\sum_{i=1}^{\frac nk}i^2\varphi(i) \end{aligned}\]最后一步是一个结论，来证明一下：
首先有 \(\dfrac{i\varphi(i)}2 = \sum_{j=1}^{i}j[\gcd(i,j)=1]~(i\neq 1)\)，因为若 \(k\) 与 \(i\) 互素那么 \(i-k\) 也与 \(i\) 互素，所以 \(i\) 以下所有与 \(i\) 互素的数的平均数为 \(\dfrac{i}2\)，乘上总数 \(\varphi(i)\) 即得证。
那么上式中 \(i,j\) 没有大小关系，所以要 \(*2\)，就和分母上的 \(/2\) 消掉了。
但问题是 \(i=j\) 得情况算重复了一次：
- 若 \(i=j\neq 1\)，则不会对答案产生贡献，故不用管。
- 若 \(i=j=1\)，此时正好 \(\dfrac{i\varphi(i)}2\) 其实等于 \(\dfrac{1}2\)，两个一加正好对了。
所以我们就证明了这个式子。
//P3700
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 4e6 + 10;
const ll P = 1e9 + 7;
int m, n;
ll qp(ll x, ll y){
ll ans = 1;
while(y){
if(y & 1){
ans = ans * x % P;
}
x = x * x % P;
y >>= 1;
}
return ans;
}
int pr[N], cnt, vis[N], phi[N];
ll g[N];
void euler(){
phi[1] = 1;
for(int i = 2; i <= n; ++ i){
if(!vis[i]){
pr[++cnt] = i;
phi[i] = i-1;
}
for(int j = 1; j <= cnt && i * pr[j] <= n; ++ j){
vis[i*pr[j]] = 1;
if(i % pr[j]){
phi[i*pr[j]] = phi[i] * phi[pr[j]];
} else {
phi[i*pr[j]] = phi[i] * pr[j];
break;
}
}
}
for(int i = 1; i <= n; ++ i){
g[i] = (ll)phi[i] * i % P * i % P;
g[i] = (g[i] + g[i-1]) % P;
}
}
int inb[N], bl, L[N], R[N];
ll vl[N], bv[N], fr[N];
void blc(){
bl = sqrt(n);
for(int i = 1; i <= n; ++ i){
fr[i] = (ll)i * i % P;
vl[i] = ((ll)i * i % P + vl[i-1]) % P;
inb[i] = (i - 1) / bl + 1;
if(inb[i] != inb[i-1]){
R[inb[i-1]] = i-1;
L[inb[i]] = i;
}
}
R[inb[n]] = n;
}
void add(int x, ll v){
for(int i = x; i <= R[inb[x]]; ++ i){
vl[i] = (vl[i] + v) % P;
}
for(int i = inb[x]+1; i <= inb[n]; ++ i){
bv[i] = (bv[i] + v) % P;
}
fr[x] = (fr[x] + v) % P;
}
ll qry(ll x){
return (vl[x] + bv[inb[x]]) % P;
}
ll qry(ll l, ll r){
return (qry(r) - qry(l-1) + P) % P;
}
int main(){
scanf("%d%d", &m, &n);
euler();
blc();
while(m--){
int a, b, k;
ll x;
scanf("%d%d%lld%d", &a, &b, &x, &k);
int d = __gcd(a, b);
x = x % P * qp(a/d, P-2) % P * qp(b/d, P-2) % P;
add(d, (x - fr[d] + P) % P);
ll ans = 0;
for(int l = 1, r; l <= k;){
r = k / (k / l);
ans = (ans + qry(l, r) * g[k/l]) % P;
l = r + 1;
}
printf("%lld\n", ans);
}
return 0;
}