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;}elsephi[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;
}