P3700 [CQOI2017] 小 Q 的表格 题解

这道题是我的第九道黑题，也是我第二道完全独立做出的黑题。

首先考虑哪些格子会互相修改，也就是修改了其中的一个就会修改另一些。

观察第一个式子 \(f(a,b)=f(b,a)\)，感觉没啥用，它只告诉了我们整个表格是对称的。

所以考虑在第二个式子 \(b \times f(a,a+b) = (a+b) \times f(a,b)\)，感觉很复杂，考虑转化一下式子。

不难注意到这个式子的右边有 \(a+b\)，左边也有 \(a+b\)。所以考虑转换成一个分数等于一个分数的形式：

\[\frac{f(a,a+b)}{a+b} = \frac{f(a,b)}{b} \]

感觉还是少了什么，怎么办？观察到 \(f(a,a+b)\) 和 \(f(a,b)\) 里面都有一个 \(a\)，所以：

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

这样感觉就有些对劲了。然后怎么办？？？

考虑观察式子，我能不能把 \(a\) 换成 \(b\)，把 \(b\) 换成 \(a-b\)：

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

熟悉数论的同学就发现了，如果我们这么一直换，这不就是更相减损术吗！所以，我们最终会换出来：

\[\frac{f(a,b)}{ab}=\frac{f(g,g)}{g^2},g=\gcd(a,b) \]

再变换一下：

\[f(a,b)=\frac{ab}{g^2} f(g,g),g=\gcd(a,b) \]

再观察一下这个式子。如果我们改变了 \(f(a,b)\)，我们就会改变 \(f(g,g)\)，而 \(f(g,g)\) 又会连锁改变什么呢……

于是，我们有了一个发现：

对于 \(\gcd(a,b) = \gcd(x,y)\)，\(f(a,b)\) 和 \(f(x,y)\) 会产生连锁改变。

---

至此我们发现上面的发现和公式就已经可以完整概括题目中的两个式子了，所以我们可以忽略题目中给出的两个式子，考虑答案是什么。

首先答案是这个玩意：

\[\sum_{i=1}^k \sum_{j=1}^k \frac{ij}{(i,j)^2} f((i,j),(i,j)) \]

这里为了简化写法，将 \(f(d,d)\) 写为 \(f(d)\)。

\[\sum_{d=1}^k f(d) \sum_{i=1}^{\lfloor \frac{k}{d}\rfloor} \sum_{j=1}^{\lfloor \frac{k}{d}\rfloor} \frac{id \times jd}{d^2}[(i,j)=1]=\sum_{d=1}^k f(d) \sum_{i=1}^{\lfloor \frac{k}{d}\rfloor} \sum_{j=1}^{\lfloor \frac{k}{d}\rfloor} ij[(i,j)=1] \]

这个时候看见了 \([(i,j)=1]\)，好，\(\mu\) 函数启动，但是千万不要往莫比乌斯上面想，这会真的使问题超级复杂化的！

看见互质，但是没有莫比乌斯反演可以想，怎么办呢？？？考虑使用欧拉函数简化这个式子的后面那一部分。

---

欧拉函数只能处理互质的个数，但如果需要处理互质的数的和，该怎么办呢？先考虑原来的式子：

\[\sum_{i=1}^{x} \sum_{j=1}^{x} ij[(i,j)=1] \]

（这里使用 \(x\) 来代替 \(\lfloor \frac{k}{d}\rfloor\)，因为后者太占据篇幅了。）

因为 \(2 \le i=j\) 时 \(ij[(i,j)=1]\) 这个东西为 \(0\)，所以我们可以忽略，只处理 \(i>j\) 和 \(i<j\) 的情况，但最后还要因为 \((1,1)=1\) 再加上 \(1\)（因为 \(i=1\) 时 \(i=j\) 同样不能忽略）：

\[\sum_{i=1}^x \sum_{j=1}^{i-1} ij[(i,j)=1]+\sum_{i=1}^x \sum_{j=i+1}^x ij[(i,j)=1]+1 \]

很容易发现后面的一部分实际上和前面的一部分是一模一样的。

\[2\sum_{i=1}^x \sum_{j=1}^{i-1} ij[(i,j)=1]+1 \]

---

为了求解这个式子，我们先来看另一个式子。我们现在需要简化这个式子：

\[\begin{aligned} \sum_{i=1}^n i [(i,n)=1]&=\sum_{i=1}^n i \sum_{d|(i,n)} \mu(d) \\ &=\sum_{d|n} \mu(d) \sum_{d|i}^n i \\ &=\sum_{d|n} \mu(d) \sum_{i=1}^{\frac{n}{d}} id \\ &=\sum_{d|n} \mu(d) d \sum_{i=1}^{\frac{n}{d}} i \\ &=\sum_{d|n} \mu(d) d \frac{\frac{n}{d}(\frac{n}{d}+1)}{2} \\ &=\frac{n}{2} \sum_{d|n} \mu(d) (\frac{n}{d}+1) \\ &=\frac{n}{2} (n\sum_{d|n} \frac{\mu(d)}{d}+\sum_{d|n} \mu(d) ) \end{aligned} \]

因为 \(\sum_{d|n} \mu(d) =\epsilon(n),\sum_{d|n} \frac{\mu(d)}{d}=\frac{\varphi(n)}{n}\)：

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

回到原来的式子。

\[2\sum_{i=1}^x \sum_{j=1}^{i-1} ij[(i,j)=1]+1 \]

可以把上面的式子套到这个式子里面，得到：

\[\sum_{i=1}^x i^2 \varphi(i) \]

至此，我们已经把式子尽可能地简化了。

---

我们从原来的式子：

\[\sum_{d=1}^k f(d) \sum_{i=1}^{\lfloor \frac{k}{d}\rfloor} \sum_{j=1}^{\lfloor \frac{k}{d}\rfloor} ij[(i,j)=1] \]

变成了这个式子：

\[\sum_{d=1}^k f(d) \sum_{i=1}^{\lfloor \frac{k}{d}\rfloor} i^2 \varphi(i) \]

是不是很神奇？

而且，观察到后面的一部分是个定值，这意味着我们可以预处理它！这里，不妨设 \(S(x)=\sum_{i=1}^x i^2 \varphi(i)\)。

所以这个式子又可以进一步缩短：

\[\sum_{d=1}^k f(d) S(\lfloor \frac{k}{d}\rfloor) \]

我们的答案就是这个式子，感觉已经简单了许多，这都是推式子的功劳。

---

考虑计算上面的式子，观察发现可以使用整除分块。 那么我们单次询问的复杂度就是 \(O(\sqrt k) \approx O(\sqrt n)\) 的，而 \(n\) 很大，所以我们不能让复杂度更大了。

那么我们就需要很快的求出 \(f\) 的一段区间的区间和，秉承着不能增高复杂度的原则，求一段区间和的时间复杂度必须是 \(O(1)\)，这样有些严格了。

而我们在修改的时候默认是 \(O(1)\) 的，因为我们要修改 \((a,b)\) 只需要修改 \((d,d),d=\gcd(a,b)\) 这一个位置。而因为我们计算答案的时候复杂度已经至少是 \(O(\sqrt n)\) 了，所以我们让修改变成 \(O(\sqrt n)\) 的复杂度也没有关系。

既然我们区间查询必须是 \(O(1)\) 的，而单点修改却绰绰有余，那么我们可不可以均摊一下，将单点修改变成 \(O(\sqrt n)\) 呢？

于是，我们现在需要的就是一个单点修改 \(O(\sqrt n)\)，区间查询 \(O(1)\) 的数据结构，每一次查询一个区间的和。显然会想到分块。

分块维护前缀和数组即可。

#include <bits/stdc++.h>
#define int long long
using namespace std;
int m, n;
const int N = 4e6+10, mod = 1e9+7;
int len;
int pos[N];
int t[N], pre[N];
int nw[N];

int back(int x) {
	if (x > n)
		return n + 1;
	return min(n, pos[x] * len);
}

int front(int x) {
	return (pos[x] - 1) * len + 1;
}//计算块的左右断电

void upd(int x, int v) {
	if (x == front(x))
		t[pos[x]] = (t[pos[x]] + v + mod) % mod;
	else
		for (int i = x; i <= back(x); i++)
			pre[i] = (pre[i] + v + mod) % mod;
	for (int i = pos[x] + 1; i <= pos[n]; i++)
		t[i] = (t[i] + v + mod) % mod;
}

int query(int x) {
	return (t[pos[x]] + pre[x]) % mod;
}

void init() {
	len = sqrt(n);
	for (int i = 1; i <= n; i++)
		pos[i] = (i - 1) / len + 1, nw[i] = i * i % mod, pre[i] = (pre[i - 1] + nw[i]) % mod;
}
int phi[N], ps[N], p[N];
bool f[N];

void sieve() {
	phi[1] = ps[1] = 1;
	for (int i = 2; i <= n; i++) {
		if (!f[i])
			p[++p[0]] = i, phi[i] = i - 1;
		for (int j = 1; j <= p[0] && i * p[j] <= n; j++) {
			f[i * p[j]] = 1;
			if (i % p[j])
				phi[i * p[j]] = (p[j] - 1) * phi[i];
			else {
				phi[i * p[j]] = p[j] * phi[i];
				break;
			}
		}
		ps[i] = (ps[i - 1] + phi[i] % mod *i % mod *i % mod) % mod;
	}
}//线性筛

signed main() {
	cin >> m >> n;
	sieve(), init();
	while (m--) {
		int a, b, x, k;
		cin >> a >> b >> x >> k;
		int y = __gcd(a, b), bl = x / (a / y) / (b / y);
		bl %= mod, upd(y, -nw[y]), upd(y, bl), nw[y] = bl;
		int ans = 0;
		for (int l = 1, r = 0; l <= k && r <= k; l = r + 1) {
			r = k / (k / l);
			ans = (ans + ps[k / l] % mod*(query(r) - query(l - 1) + mod) % mod) % mod;
		}
		cout << ans << endl;
	}
	return 0;
}