[CQOI2017] 小Q的表格（分块 + 整除分块 + 数学 + 前缀和）

problem

luogu-P3700

solution

• $f(a,b)=f(b,a)$ 意味着我们只用考虑半个棋盘的信息。

• $b\ast f(a,a+b)=(a+b)\ast f(a,b)$

  会发现修改 $f(a,b)$ 就影响 $f(a,a+b)$ 进而影响 $f(a,a+2b)\dots$

  对于这个信息，如果对数学很敏感的选手，就会直接猜跟 $\gcd$ 挂钩了。

  向我们这种凡人呢就只能听天由命了

  考虑这个性质之所以很难使用上，是因为这两个式子的系数只有一项，且是与等式对面的 $f(x,y)$ 挂钩的。

  这种跟对面挂钩的我们会想到比例式，因为我们清楚如果分母和分子相同是可以约分的。

  当写成比例式 $\frac{f(a,a+b)}{a+b}=\frac{f(a,b)}{b}$，我们才会稍微察觉到如果两边都再除以一个 $a$，整体形式就很同构了。

  $\frac{f(a,a+b)}{a\ast (a+b)}=\frac{f(a,b)}{a\ast b}$ 加减运算不影响 $\gcd$ 的值。

  递归辗转相除最后就会得到 $\frac{f(a,b)}{a\ast b}=\frac{f(\gcd (a,b),\gcd (a,b))}{\gcd (a,b{)}^2}$。

  所以我们可以得出：修改 $f(x,y)$ 只会影响 $\gcd (x,y)=\gcd (a,b)$ 的所有 $f(a,b)$ 的值。

---

继续考虑处理询问，这就要推式子了。
$$\sum _{i=1}^n\sum _{j=1}^{n}f(i,j)=\sum _{i=1}^n\sum _{j=1}^n\frac{i\times j}{\gcd (i,j{)}^2}f(i,j)=\sum _{d=1}^{n}f(d,d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }i\times j\cdot [\gcd (i,j)=1]$$

$$\begin{gathered} \sum _{k=1}^n[\gcd (k,n)=1]\cdot k \\ =\sum _{k=1}^{n}k\sum _{d\mid k}\mu (d)=\sum _{d\mid n}\mu (d)\sum _{d\mid k}k=\sum _{d\mid n}\mu (d)\ast d\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }i \\ =\sum _{d\mid n}\mu (d)\ast d\ast \frac{(\lfloor \frac{n}{d}\rfloor +1)\ast \lfloor \frac{n}{d}\rfloor }{2}=\frac{n\cdot ({\sum }_{d\mid n}\mu (d)+n{\sum }_{d\mid n}\frac{\mu (d)}{d})}{2} \\ =\frac{n\cdot (\phi (n)+[n=1])}{2} \end{gathered}$$

根据两个基础数论结论：${\sum }_{d\mid n}\mu (d)=[n=1],{\sum }_{d\mid n}\frac{\mu (d)}{d}=\frac{\phi (n)}{n}$ 得出最后一个等式。

$$\sum _{d=1}^{n}f(d,d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }i\times j\cdot [\gcd (i,j)=1]\ne \sum _{d=1}^{n}f(d,d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }i\cdot \frac{i\cdot (\phi (i)+[i=1])}{2}$$

这个式子之所以不能直接用上面的结论，是因为这里的 $j$ 可能大于 $i$，根据对称性，$\times 2$ 即可。

但这样就会把所有 $i=j$ 的情况重复算了一遍，但除了 $i=j=1$ 外其余都不满足最大公约数为 $1$ 的限制。

$i=j=1$ 的贡献是 $1$，而 $\frac{1\times \phi (1)+[i=1]}{2}\times 2=2$，我们可以直接扔掉 $[i=1]$，这样这种情况算出来的贡献就是对的。

且对于 $i>1$ 的情况去掉 $[i=1]$ 的特判没有影响。
$$=\sum _{d=1}^{n}f(d,d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }i\cdot \frac{i\ast \phi (i)}{2}\times 2=\sum _{d=1}^{n}f(d,d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }{i}^2\phi (i)$$
记 $S(n)={\sum }_{i=1}^n{i}^2\phi (i)$，则答案即为 ${\sum }_{d=1}^{n}f(d,d)S(\lfloor \frac{n}{d}\rfloor )$。

然后似乎就可以整除分块了，然而整除分块还要用到 $f$ 的前缀和。

很遗憾，这个 $f$ 还要是支持修改的，所以需要一个数据结构。

树状数组查询 $O(m\sqrt{n}\log n)$ 不太行，但修改却只有一个 $\log$ 很快。当然这也能过。

我们需要平衡一下查询和修改的时间复杂度，那就是分块啦！

直接把 $f$ 变成 $f$ 的前缀和，单点修改变成区间修改 $\gcd (a,b)\sim n$，之后查询前缀和就变成了单点询问了。

具体而言 $f(a,b)\leftarrow x$，即为 $f(\gcd (a,b),\gcd (a,b))\leftarrow \frac{x\cdot \gcd (a,b{)}^2}{ab}$（这个是由上面最开始递归辗转相除的式子变形得来的）

code

#include <bits/stdc++.h>
using namespace std;
#define int long long 
#define mod 1000000007
#define maxn 4000005
#define maxB 2005
int m, n, cnt;
int a[maxn], block[maxn], L[maxB], R[maxB], f[maxn], g[maxB], phi[maxn], s[maxn], prime[maxn], vis[maxn];

int gcd( int x, int y ) {
	if( ! y ) return x;
	else return gcd( y, x % y );
}

void init() {
	int B = sqrt( n );
	for( int i = 1;i <= n;i ++ ) {
		a[i] = i * i % mod;
		f[i] = ( f[i - 1] + a[i] ) % mod;
		block[i] = ( i - 1 ) / B + 1;
	}
	for( int i = n;i >= 1;i -- ) L[block[i]] = i; 
	for( int i = 1;i <= n;i ++ ) R[block[i]] = i;
}
void modify( int l, int r, int x ) {
	int k = x - a[l]; a[l] = x;
	if( block[l] == block[r] )
		for( int i = l;i <= r;i ++ ) f[i] = ( f[i] + k ) % mod;
	else {
		for( int i = l;i <= R[block[l]];i ++ ) f[i] = ( f[i] + k ) % mod;
		for( int i = L[block[r]];i <= r;i ++ ) f[i] = ( f[i] + k ) % mod;
		for( int i = block[l] + 1;i < block[r];i ++ ) g[i] = ( g[i] + k ) % mod;
	}
}
int query( int x ) { return ( f[x] + g[block[x]] ) % mod; }

void sieve() {
	phi[1] = 1;
	for( int i = 2;i <= n;i ++ ) {
		if( ! vis[i] ) prime[++ cnt] = i, phi[i] = i - 1;
		for( int j = 1;j <= cnt and i * prime[j] <= n;j ++ ) {
			vis[i * prime[j]] = 1;
			if( i % prime[j] == 0 ) {
				phi[i * prime[j]] = phi[i] * prime[j] % mod;
				break;
			}
			else
				phi[i * prime[j]] = phi[i] * phi[prime[j]] % mod;
		}
	}
	for( int i = 1;i <= n;i ++ ) s[i] = ( s[i - 1] + i * i % mod * phi[i] ) % mod;
}
int solve( int n ) {
	int ans = 0;
	for( int l = 1, r;l <= n;l = r + 1 ) {
		r = n / ( n / l );
		ans = ( ans + ( query( r ) - query( l - 1 ) ) * s[n / l] ) % mod;
	}
	return ans;
}

signed main() {
	scanf( "%lld %lld", &m, &n );
	init();
	sieve();
	for( int i = 1, a, b, x, k;i <= m;i ++ ) {
		scanf( "%lld %lld %lld %lld", &a, &b, &x, &k );
		int d = gcd( a, b );
		x = x / ( a / d ) / ( b / d ) % mod;
		modify( d, n, x );
		printf( "%lld\n", ( solve( k ) + mod ) % mod );
	}
	return 0;
}