【洛谷P3700】小Q的表格

题目

题目链接：https://www.luogu.com.cn/problem/P3700
小Q是个程序员。
作为一个年轻的程序员，小Q总是被老C欺负，老C经常把一些麻烦的任务交给小Q来处理。每当小Q不知道如何解决时，就只好向你求助。
为了完成任务，小Q需要列一个表格，表格有无穷多行，无穷多列，行和列都从1开始标号。为了完成任务，表格里面每个格子都填了一个整数，为了方便描述，小Q把第a行第b列的整数记为f(a,b)。为了完成任务，这个表格要满足一些条件：

1. 对任意的正整数a,b，都要满足f(a,b)=f(b,a)；
2. 对任意的正整数a,b，都要满足b×f(a,a+b)=(a+b)×f(a,b)。

为了完成任务，一开始表格里面的数很有规律，第a行第b列的数恰好等于a×b，显然一开始是满足上述两个条件的。为了完成任务，小Q需要不断的修改表格里面的数，每当修改了一个格子的数之后，为了让表格继续满足上述两个条件，小Q还需要把这次修改能够波及到的全部格子里都改为恰当的数。由于某种神奇的力量驱使，已经确保了每一轮修改之后所有格子里的数仍然都是整数。为了完成任务，小Q还需要随时获取前k行前k列这个有限区域内所有数的和是多少，答案可能比较大，只需要算出答案mod1,000,000,007之后的结果。
\(1 ≤ m ≤ 10^4, 1 ≤ a,b,k ≤ n ≤ 4×10^6, 0 ≤ x < 10^{18}\)。

思路

太神啦 orz。
不难发现，如果我们更改 \(f(x,y)\)，那么显然会连带着更改 \(f(x,y+kx),f(x+ky,y)\)。
发现这个东西很想更相减损法，其实不难发现，我们更改 \(f(x,y)\) 后，所有 \(\gcd(i,j)=\gcd(x,y)\) 的 \(f(i,j)\) 都会被更改。
既然如此，令 \(d=\gcd(x,y)\)，考虑寻找 \(f(d,d)\) 与 \(f(x,y)\) 之间的关系。由条件二可知，\(f(d,kd)=k\times f(d,d)\)。拓展可得

\[f(x,y)=\frac{xy}{d^2}f(d,d) \]

那么我们更新 \(f(x,y)\) 的时候只需要更新 \(f(d,d)\) 即可。因为其他所有要更改的格子的 \(\gcd\) 都是 \(d\)，都可以通过 \(f(d,d)\) 推算。
那么设 \(g(i)=f(i,i)\)，对于一次修改操作，直接把 \(g(\gcd(x,y))\) 设为 \(\frac{\gcd(x,y)^2}{xy}c\) 即可。其中 \(c\) 是修改的数值。
对于一次询问操作，答案是

\[ans=\sum^{k}_{i=1}\sum^{k}_{j=1}\frac{ij}{\gcd(i,j)^2}f(\gcd(i,j)) \]

套路性把 \(\gcd(i,j)\) 提出来

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

这里我们需要知道一个结论：

\[\sum^{n}_{i=1}[\gcd(i,n)=1]i=\left\{\begin{matrix}\frac{n\varphi(n)}{2}\ (n>1)\\1\ (n=1)\end{matrix}\right. \]

证明非常简单，如果 \(\gcd(i,n)=1\)，那么 \(\gcd(n-i,n)=1\)，也就是对于任意一个与 \(n\) 互质的数字 \(i\)，都存在另一个与 \(n\) 互质的数 \(j\) 且 \(i+j=n\)。那么一共有 \(\frac{\varphi(n)}{2}\) 对这样的数字。
回到本题，我们可以把 \(i\) 看做 \(n\)，\(j\) 看做 \(i\)，那么

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

因为我们可以把 \(j\) 的枚举范围变成 \([1,i]\)，这样其实恰好做了一半。最后再把式子乘二就得到了上式。恰好 \(1\) 的贡献不用特判了。
然后整除分块就可以了。
但是我们需要支持单点修改和求前缀和，采用树状数组的话复杂度就是 \(O(m\sqrt{n}\log n)\)，无法接受。
发现除去这一部分时，修改操作复杂度 \(O(m)\)，询问操作复杂度 \(O(m\sqrt{n})\)。我们希望一个修改 \(O(\sqrt{n})\)，询问 \(O(1)\) 的数据结构。
那么直接分块维护前缀和即可。算是一个比较经典的 trick 就不介绍了。
时间复杂度 \(O(n+m\sqrt{n})\)。

代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

const int N=4000010,M=10010,MOD=1000000007;
int n,m,T,cnt,prm[N],bel[N];
bool v[N];
ll ans,phi[N],f[N],g[N],sum[2][N];

void findprm(int n)
{
	phi[1]=1;
	for (int i=2;i<=n;i++)
	{
		if (!v[i]) prm[++cnt]=i,phi[i]=i-1;
		for (int j=1;j<=cnt;j++)
		{
			if (i>n/prm[j]) break;
			v[i*prm[j]]=1; phi[i*prm[j]]=phi[i]*(prm[j]-1);
			if (i%prm[j]==0)
			{
				phi[i*prm[j]]=phi[i]*prm[j];
				break;
			}
		}
	}
}

ll fpow(ll x,ll k)
{
	ll ans=1;
	for (;k;k>>=1,x=x*x%MOD)
		if (k&1) ans=ans*x%MOD;
	return ans;
}

int gcd(int x,int y)
{
	return y?gcd(y,x%y):x;
}

ll calc(int i)
{
	return (sum[0][bel[i]]+sum[1][i])%MOD;
}

int main()
{
	scanf("%d%d",&m,&n);
	findprm(n);
	T=sqrt(n)+1;
	for (int i=1;i<=n;i++)
	{
		bel[i]=(i-1)/T+1;
		phi[i]=(phi[i-1]+phi[i]*i%MOD*i)%MOD;
		f[i]=1LL*i*i%MOD;
		g[i]=(g[i-1]+f[i])%MOD;
		sum[1][i]=(g[i]-g[(bel[i]-1)*T])%MOD;
	}
	for (int i=1;i<=T;i++) sum[0][i]=g[(i-1)*T];
	for (int o=1;o<=m;o++)
	{
		int x,y,k; ll p;
		scanf("%d%d",&x,&y);
		scanf("%lld",&p);
		scanf("%d",&k);
		int d=gcd(x,y);
		p=p%MOD*d%MOD*d%MOD*fpow(x,MOD-2)%MOD*fpow(y,MOD-2)%MOD;
		for (int i=bel[d]+1;i<=T;i++)
			sum[0][i]=(sum[0][i]+p-f[d])%MOD;
		for (int i=d;i<=bel[d]*T;i++)
			sum[1][i]=(sum[1][i]+p-f[d])%MOD;
		f[d]=p; ans=0;
		for (int l=1,r;l<=k;l=r+1)
		{
			r=k/(k/l);
			ans=(ans+(calc(r)-calc(l-1))*phi[k/l])%MOD;
		}
		printf("%lld\n",(ans%MOD+MOD)%MOD);
	}
	return 0;
}