神仙题，神仙题
这是一道很适合盯着发呆的题目
看到这个规律
这也没什么规律啊
于是自闭了
盯着发呆一个小时之后发现，这个(f(a,a+b))和(f(a,b))有关系
因为修改((a,b))就一定会影响((a,a+b))，同时也会影响((a,a+2b)...)
停
这不是更相减损术吗
于是我们得出了第一个结论
修改(a,b)这个值只会影响(gcd(x,y)=gcd(a,b))的(f(x,y))的值
但是这样我们还是没有什么办法来维护啊，毕竟矩阵那么大，我们修改一次影响的数那么多
我们大胆猜想格子的值存在某种关系，如果(gcd(a,b)=d)，那么(f(a,b))肯定和(f(d,d))存在某种关系
尝试去求一下这个关系
显然(f(d,kd)=frac{k}{(k-1)}f(d,(k-1)d)=k imes f(d,kd))
显然纵坐标也会有这样的性质
于是(f(k_1d,k_2d),k_1perp k_2)，就会有(k_1 imes k_2 imes f(d,d)=f(k_1d,k_2d))
其实也就是这样
考虑把(f(a,b))写成(frac{ab}{d^2}f(d,d))
于是我们只需要记录(f(d,d))的值了，这样就可以处理修改操作了
接下来把(f(d,d))简记做(f(d))
现在我们要求的柿子是
考虑一下枚举(gcd)
那个除以((i,j)^2)消失了是因为我们后面乘上的是(i,j)，本来就是都除以(d)了的
之后只要记住一条，千万别反演就好了
我们能通过欧拉函数把上面的柿子写成这个样子
至于为什么，我们需要这个柿子
至于这个柿子这么来的，我们证明一下
交换一下求和符号
我们现在需要两条很基础的结论
这里就不再证明了
根据上面那条结论我们有
我们设(S(n)=sum_{i=1}^ni^2 imes varphi(i))
答案就是
我们现在就可以尽情的整除分块了
但是由于(f)需要支持修改我们还要查询前缀和，于是看起来有点自闭，因为树状数组的复杂度高达(O(msqrt{n}logn))，好像不是很科学
但是修改却快的一批，低到(O(mlogn))，考虑一个神奇的数据结构，可以做到(O(1))单点求和
自然是神奇的分块了，我们直接把(f)做成前缀和，单点修改我们直接搞成区间修改，之后我们单点查询前缀和就可以很快了
代码
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define re register
#define LL long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
const int maxn=4e6+5;
const LL mod=1e9+7;
inline int gcd(int a,int b) {return !b?a:gcd(b,a%b);}
inline LL read() {
char c=getchar();LL x=0;while(c<'0'||c>'9') c=getchar();
while(c>='0'&&c<='9') x=(x<<3)+(x<<1)+c-48,c=getchar();return x;
}
int f[maxn],p[maxn>>1];
LL phi[maxn],pre[maxn],F[maxn];
int n,m;
struct Block {
int sz,tot;
int l[3005],r[3005],id[maxn];
LL tag[3005];
inline void Build() {
sz=std::sqrt(n);
for(re int i=1;i<=n;i+=sz) {
l[++tot]=i,r[tot]=min(i+sz-1,n);
for(re int j=l[tot];j<=r[tot];j++) id[j]=tot;
}
}
inline LL ask(int x) {return (pre[x]+tag[id[x]])%mod;}
inline void change(int x,LL val) {
int j=id[x];
if(x==l[j]) tag[j]+=val,tag[j]=(tag[j]+mod)%mod;
else for(re int i=x;i<=r[j];i++) pre[i]=(pre[i]+val+mod)%mod;
j++;while(j<=tot) tag[j]+=val,tag[j]=(tag[j]+mod)%mod,j++;
}
}B;
int main() {
m=read(),n=read();
phi[1]=1;
for(re int i=2;i<=n;i++) {
if(!f[i]) p[++p[0]]=i,phi[i]=i-1;
for(re int j=1;j<=p[0]&&p[j]*i<=n;j++) {
f[p[j]*i]=1;
if(i%p[j]==0) {phi[p[j]*i]=p[j]*phi[i];break;}
phi[p[j]*i]=phi[p[j]]*phi[i];
}
}
for(re LL i=1;i<=n;i++) phi[i]=phi[i-1]+i*i%mod*phi[i]%mod,phi[i]%=mod;
for(re LL i=1;i<=n;i++) F[i]=(i*i)%mod;
for(re int i=1;i<=n;i++) pre[i]=pre[i-1]+F[i],pre[i]%=mod;
B.Build();
int a,b,k;LL v;
while(m--) {
a=read(),b=read(),v=read(),k=read();
int t=gcd(a,b);LL ans=0;
B.change(t,-1ll*F[t]);
F[t]=v/((LL)(a/t)*(LL)(b/t));
B.change(t,F[t]);
for(re int l=1,r;l<=k;l=r+1) {
r=k/(k/l);
ans+=phi[k/l]*(B.ask(r)-B.ask(l-1)+mod)%mod,ans%=mod;
}
printf("%lld
",ans);
}
return 0;
}