参考：http://blog.csdn.net/qq_33229466/article/details/70174227
看这个等式的形式就像高精gcd嘛…所以随便算一下就发现每次修改(a,b)影响到的都是横纵坐标gcd为gcd(a,b)的，进而发现可以把gcd(i,j)==d的一部分都归到d上，f(a,b)=f(d,d)ab/d/d ，这样二维就变成一维了，设为f。
然后答案就是：
[ans=sum_{d=1}^{k}f(d)sum_{i=1}^{k}sum_{j=1}^{k}[gcd(i,j)==d]frac{i*j}{d^2}
]
[ans=sum_{d=1}^{k}f(d)sum_{i=1}^{frac{k}{d}}sum_{j=1}^{frac{k}{d}}[gcd(i,j)==1]i*j
]
[s(n)=sum_{i=1}^{n}sum_{j=1}^{n}[gcd(i,j)==1]i*j=sum_{i=1}^{n}varphi(i)i^2
]
[ans=sum_{d=1}^{k}f(d)s(frac{n}{d})
]
（f应该先按照原始状态写出来
因为m比n小很多，所以用分块维护即可
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int N=4000005,mod=1000000007;
int m,n,t[N],phi[N],q[N],tot,inv[N],f[N],w[N],c[N],kuai,lz[N];
bool v[N];
int read()
{
int r=0,f=1;
char p=getchar();
while(p>'9'||p<'0')
{
if(p=='-')
f=-1;
p=getchar();
}
while(p>='0'&&p<='9')
{
r=r*10+p-48;
p=getchar();
}
return r*f;
}
int gcd(int a,int b)
{//cout<<"gcd"<<endl;
return b==0?a:gcd(b,a%b);
}
void add(int x,int v)
{//cout<<"add"<<endl;
for(int i=x;i<=(x+kuai-1)/kuai*kuai;i++)
w[i]=((w[i]+v)%mod+mod)%mod;
for(int i=(x+kuai-1)/kuai+1;i<=(n+kuai-1)/kuai;i++)
lz[i]=((lz[i]+v)%mod+mod)%mod;
}
int ques(int x)
{//cout<<"ques"<<endl;
int re=0;
for(int i=1,la;i<=x;i=la+1)
{
la=x/(x/i);
re=(re+1ll*(w[la]+lz[(la+kuai-1)/kuai]-w[i-1]-lz[(i+kuai-2)/kuai])%mod*f[x/i]%mod)%mod;
}
return (re+mod)%mod;
}
int main()
{
m=read(),n=read();
kuai=sqrt(n);
v[1]=1,phi[1]=1;
for(int i=2;i<=n;i++)
{
if(!v[i])
{
q[++tot]=i;
phi[i]=i-1;
}
for(int j=1;j<=tot&&i*q[j]<=n;j++)
{
int k=i*q[j];
v[k]=1;
if(i%q[j]==0)
{
phi[k]=phi[i]*q[j];
break;
}
phi[k]=phi[i]*(q[j]-1);
}
}
for(int i=1;i<=n;i++)
f[i]=(f[i-1]+1ll*i*i%mod*phi[i]%mod)%mod;
for(int i=1;i<=n;i++)
w[i]=(w[i-1]+1ll*i*i%mod)%mod,c[i]=1ll*i*i%mod;
while(m--)
{
int a=read(),b=read();
long long x;
scanf("%lld",&x);
int k=read(),g=gcd(a,b),gai=1ll*x/(a/g)/(b/g)%mod;
add(g,((gai-c[g])%mod+mod)%mod);
c[g]=gai;
printf("%d
",ques(k));
}
return 0;
}