BZOJ4815 [CQOI2017]小Q的表格 【数论 + 分块】

题目链接

BZOJ4815

题解

根据题中的式子，手玩一下发现和\(gcd\)很像
化一下式子：

\[\begin{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)\)

那么我们有

\[\begin{aligned} ans &= \sum\limits_{i = 1}^{k} \sum\limits_{j = 1}^{k} f(i,j) \\ &= \sum\limits_{d = 1}^{k} f(d,d) \sum\limits_{d|i} \sum\limits_{d|j} \frac{ij}{d^2} \quad [gcd(i,j) == d] \\ &= \sum\limits_{d = 1}^{k} f(d,d) \sum\limits_{i = 1}^{\lfloor \frac{k}{d} \rfloor} \sum\limits_{j = 1}^{\lfloor \frac{k}{d} \rfloor} ij \quad [i \perp j] \\ \end{aligned} \]

令

\[g(n) = \sum\limits_{i = 1}^{n} \sum\limits_{j = 1}^{n} ij \quad [i \perp j] \]

由于

\[\sum\limits_{i = 1}^{n} i \quad [i \perp n] = \frac{n\varphi(n)}{2} \]

所以

\[\begin{aligned} g(n) &= \sum\limits_{i = 1}^{n} i \times 2 \times \frac{i\varphi(i)}{2} \\ &= \sum\limits_{i = 1}^{n} i^2\varphi(i) \\ \end{aligned} \]

我们可以线性筛\(O(n)\)预处理出\(g(n)\)
对于答案的式子，可以\(O(\sqrt{k})\)整除分块
所以我们只需要\(O(1)\)计算\(f(d,d)\)的前缀和

分块即可
块外维护块的前缀和，块内维护块内前缀和
这样修改是\(O(\sqrt{n})\)的，修改复杂度\(O(m\sqrt{n})\)
且询问时\(O(1)\)的

总复杂度\(O(m\sqrt{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\n",(ans % P + P) % P);
	}
	return 0;
}