P3700 [CQOI2017]小Q的表格

P3700 [CQOI2017]小Q的表格

有一个表格，满足：

$$\begin{aligned} &f(a,b)=f(b,a)\\ &b \times f(a,a+b)=(a+b)\times f(a,b) \end{aligned}$$

且一开始 $f(a,b)=a\times b$，然后带上一个单点修改操作。

每当修改了一个格子的数之后，为了让表格继续满足上述两个条件，你还需要把这次修改能够波及到的全部格子里都改为恰当的数。

你还需要随时输出前 $k$ 行前 $k$ 列这个有限区域内所有数的和 $\pmod {10^9+7}$。

考虑从下面那个式子得出一些可以反演的东西，有：

$$\begin{aligned} b \times f(a,a+b)&=(a+b)\times f(a,b)\\ ab \times f(a,a+b)&=a(a+b)\times f(a,b)\\ \frac{f(a,a+b)}{a(a+b)}&=\frac{f(a,b)}{ab}\\ \frac{f(a,b)}{ab}&=\frac{f(\gcd(a,b),\gcd(a,b))}{\gcd(a,b)^2} \end{aligned}$$

设 $\gcd(a,b)=d$，则有：

$$\begin{aligned} \frac{f(a,b)}{ab}&=\frac{f(d,d)}{d^2}\\ f(a,b)&=\frac{f(d,d)\times ab}{d^2} \end{aligned}$$

考虑到这个矩阵要变换，因为修改一个位置 $(a,b)$ 的值后有且仅有 $\gcd(x,y)=d$ 的 $(x,y)$ 位置的 $f$ 值会跟着变化，那么用一个树状数组维护 $f(d,d)$ 即可。

设 $f(d,d)=g(d)$，再考虑答案：

$$\begin{aligned} ans&=\sum_{d=1}^{n}g(d)\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{ij}{\gcd(i,j)^2}[\gcd(i,j)=d]\\ &=\sum_{d=1}^{n}g(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}ij[\gcd(i,j)=1]\\ &=\sum_{d=1}^{n}g(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}ij\sum_{p|\gcd(i,j)}\mu(p)\\ &=\sum_{d=1}^{n}g(d)\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p) \sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}i[p|i]\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}j[p|j]\\ &=\sum_{d=1}^{n}g(d)\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)p^2\left(\sum_{i=1}^{\lfloor\frac{\lfloor\frac{n}{d}\rfloor}{p}\rfloor}i\right)^2\\ &=\sum_{d=1}^{n}g(d)\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)p^2S(\lfloor\frac{\lfloor\frac{n}{d}\rfloor}{p}\rfloor) ^2 \end{aligned}$$

设 $G(n)=\sum\limits_{i=1}^{n}i^2\mu(i)S(\lfloor\frac{n}{i}\rfloor)^2$，则原式等于：

$$\sum_{d=1}^{n}g(d)G(\lfloor\frac{n}{d}\rfloor)$$

根据 $\lfloor\frac{n}{i}\rfloor-\lfloor\frac{n-1}{i}\rfloor=[i|n]$，有：

$$\begin{aligned} G(n)-G(n-1)&=\sum_{i|n}i^2\mu(i)(S(\lfloor\frac{n}{i}\rfloor)^2-S(\lfloor\frac{n}{i}\rfloor-1)^2\\ &=\sum_{i|n}i^2\mu(i)(\frac{n}{i})^3\\ &=n^2\sum_{i|n}\mu(i)\frac{n}{i}\\ &=n^2\varphi(n) \end{aligned}$$

那么有：

$$G(n)=\sum_{i=1}^{n}i^2\varphi(i)$$

那么原式等于：

$$\sum_{d=1}^{n}g(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}i^2\varphi(i)$$

线性筛，维护前缀和，再数论分块即可。

#include <bits/stdc++.h>

using namespace std;

#define int long long

const int mod = 1e9 + 7;

inline int read()
{
	int x = 0, f = 1;
	char c = getchar();
	while(c < '0' || c > '9')
	{
		if(c == '-') f = -1;
		c = getchar();
	}
	while(c >= '0' && c <= '9')
	{
		x = ((x << 1) + (x << 3) + (c ^ 48)) % mod;
		c = getchar();
	}
	return x * f;
}

inline void write(int x)
{
	if(x < 0)
	{
		putchar('-');
		x = -x;
	}
	if(x > 9)
		write(x / 10);
	putchar(x % 10 + '0');
}

const int _ = 4e6 + 10;

int n, m;

int cnt, pr[_], vis[_], phi[_];

int f[_];

int c[_];

int val[_];

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

void inc(int &x, int y)
{
	x += y;
	if(x >= mod) x -= mod;
	else if(x < 0) x += mod;
}

int add(int x, int y)
{
	int temp = x + y;
    if (temp >= mod)
        temp -= mod;
    else if (temp < 0)
        temp += mod;
    return temp;
}

int qpow(int x, int y)
{
	int res = 1;
	while(y)
	{
		if(y & 1) res = res * x % mod;
		x = x * x % mod;
		y >>= 1;
	}
	return res;
}

void init()
{
	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] == 0)
			{
				phi[i * pr[j]] = phi[i] * pr[j];
				break;
			}
			phi[i * pr[j]] = phi[i] * (pr[j] - 1);
		}
	}
	for(int i = 1; i <= n; ++i)
		f[i] = (f[i - 1] + i * i % mod * phi[i] % mod) % mod;
}

#define lowbit(x) (x & -x)

void update(int x, int val)
{
	while(x <= n)
	{
		inc(c[x], val);
		x += (x & -x);
	}
}

int query(int x)
{
	int res = 0;
	while(x)
	{
		inc(res, c[x]);
		x -= (x & -x);
	}
	return res;
}

int Query(int l, int r)
{
	return add(query(r), -query(l - 1));
}

signed main()
{
	m = read(), n = read();
	init();
	for(int i = 1; i <= n; ++i)
	{
		val[i] = i * i % mod;
		update(i, val[i]);
	}
	while(m--)
	{
		int a = read(), b = read(), x = read(), k = read(), d = Gcd(a, b);
        int upt = x * d % mod * d % mod * qpow(a * b % mod, mod - 2ll) % mod;
        update(d, add(upt, -val[d]));
        val[d] = upt;
        int ans = 0;
        int S = sqrt(k);
        for(int i = 1; i < S; ++i) inc(ans, f[k / i] * Query(i, i) % mod);
        for (int l = S, r; l <= k; l = r + 1) 
		{
            r = k / (k / l);
            inc(ans, f[k / l] * Query(l, r) % mod);
        }
        write(ans);
        putchar('\n');
	}
	return 0;
}