题目描述
小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之后的结果。
题解
这题思路太神了。
乍一看需要维护一个二维的东西看起来很麻烦，但实际上是可以转化为一维上的问题的。
先来看第一个限制，关于对角线对称的数字相等，这样我们的棋盘大小变为原来的一半。
然后看第二个限制，挪一下位置，变成了
f(a,a+b)/(a+b)=f(a,b)/b
f(a,a+b)/((a+b)*a)=f(a,b)/(a*b)
我们发现了规律，f(a,b)是和a*b正相关的，而且这种关系非常像辗转相减，那么一直减到最后就会变成f(gcd,gcd)。
于是我们得出了结论，每个位置(a,b)的值之和gcd(a,b)相关。
于是它变成了一个序列上的问题，每次修改只需要修改一个数。
然后考虑答案。
ans=∑f(g,g)∑∑i*j*(gcd(i,j)==g) i,j,g<=k
后面的那个东西∑∑i*j*(gcd(i,j)==g)经过推导后可以得到∑i*i*Φ(i) i<=k/g
然后就可以除法分块了。
但我们还需要支持单点修改，区间求和，可以用树状数组，但是这题的数据范围比较特别，所以用分块搞就可以了，根号修改，O(1)查询。
代码
#include<iostream>
#include<cstdio>
#include<cmath>
#define N 4000009
using namespace std;
typedef long long ll;
const int mod=1e9+;
ll n,m,sum[N],pre[N],phi[N],prime[N],f[N],val[N],n1,be[N];
bool vis[N];
inline ll rd(){
ll x=;char c=getchar();bool f=;
while(!isdigit(c)){if(c=='-')f=;c=getchar();}
while(isdigit(c)){x=(x<<)+(x<<)+(c^);c=getchar();}
return f?-x:x;
}
inline void work(ll x,ll y){
pre[x]=(pre[x]+y)%mod;
for(int i=x+;be[i]==be[i-];++i)pre[i]=(pre[i]+y)%mod;
for(int i=be[x];i<=be[n];++i)sum[i]=(sum[i]+y)%mod;
}
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
inline ll query(int x){return (pre[x]+sum[be[x]-])%mod;}
inline void prework(){
phi[]=;
for(int i=;i<=n;++i){
if(!vis[i])prime[++prime[]]=i,phi[i]=i-;
for(int j=;j<=prime[]&&i*prime[j]<=n;++j){
vis[i*prime[j]]=;
if(i%prime[j]==){phi[i*prime[j]]=phi[i]*prime[j];break;}
phi[i*prime[j]]=phi[i]*phi[prime[j]];
}
}
for(int i=;i<=n;++i)f[i]=(f[i-]+phi[i]*val[i])%mod;
}
int main(){
m=rd();n=rd();n1=sqrt(n);
for(int i=;i<=n;++i){
be[i]=(i-)/n1+;val[i]=1ll*i*i%mod;
(sum[be[i]]+=val[i])%=mod;
if(be[i]==be[i-])pre[i]=(pre[i-]+val[i])%mod;
else pre[i]=val[i];
}
for(int i=;i<=be[n];++i)(sum[i]+=sum[i-])%=mod;
prework();ll a,b,x,k;
while(m--){
a=rd();b=rd();x=rd();k=rd();ll g=gcd(a,b);ll ans=,r;
x=x/((a/g)*(b/g));x%=mod;work(g,(x-val[g]+mod)%mod);val[g]=x;
for(ll l=;l<=k;l=r+){
r=k/(k/l);
ans+=(query(r)-query(l-))*f[k/l]%mod;
ans=(ans%mod+mod)%mod;
}
printf("%lld\n",ans);
}
return ;
}