题目
题目链接：https://www.luogu.com.cn/problem/P3700
小Q是个程序员。
作为一个年轻的程序员，小Q总是被老C欺负，老C经常把一些麻烦的任务交给小Q来处理。每当小Q不知道如何解决时，就只好向你求助。
为了完成任务，小Q需要列一个表格，表格有无穷多行，无穷多列，行和列都从1开始标号。为了完成任务，表格里面每个格子都填了一个整数，为了方便描述，小Q把第a行第b列的整数记为f(a,b)。为了完成任务，这个表格要满足一些条件：
- 对任意的正整数a,b，都要满足f(a,b)=f(b,a)；
- 对任意的正整数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 imes 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) 是修改的数值。
对于一次询问操作，答案是
套路性把 (gcd(i,j)) 提出来
这里我们需要知道一个结论：
证明非常简单，如果 (gcd(i,n)=1)，那么 (gcd(n-i,n)=1)，也就是对于任意一个与 (n) 互质的数字 (i)，都存在另一个与 (n) 互质的数 (j) 且 (i+j=n)。那么一共有 (frac{varphi(n)}{2}) 对这样的数字。
回到本题，我们可以把 (i) 看做 (n)，(j) 看做 (i)，那么
因为我们可以把 (j) 的枚举范围变成 ([1,i])，这样其实恰好做了一半。最后再把式子乘二就得到了上式。恰好 (1) 的贡献不用特判了。
然后整除分块就可以了。
但是我们需要支持单点修改和求前缀和，采用树状数组的话复杂度就是 (O(msqrt{n}log n))，无法接受。
发现除去这一部分时，修改操作复杂度 (O(m))，询问操作复杂度 (O(msqrt{n}))。我们希望一个修改 (O(sqrt{n}))，询问 (O(1)) 的数据结构。
那么直接分块维护前缀和即可。算是一个比较经典的 trick 就不介绍了。
时间复杂度 (O(n+msqrt{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
",(ans%MOD+MOD)%MOD);
}
return 0;
}