首先观察$b*f(a,a+b)=(a+b)*f(a,b)$这个东西
可以化成$\frac{f(a,a+b)}{a+b}=\frac{f(a,b)}{b}$，发现这类似辗转相除求gcd
而我们两边同乘一个a就能得到$\frac{f(a)}{a}$是个定值的这个结论
那么有$f(a,b)=\frac{a*b}{gcd(a,b)^2}*f(gcd(a,b),gcd(a,b))$
为了方便现在设$gcd(i,j)=g$，现在把这个东西放进原来的式子里
$\sum\limits_{i=1}^k\sum\limits_{j=1}^kf(i,j)$
$=\sum\limits_{i=1}^k\sum\limits_{j=1}^k\frac{i*j}{g^2}*f(g,g)$
改为枚举$g$，把$f(g,g)$前提
$=\sum\limits_{d=1}^kf(d,d)\sum_{d|i}\sum_{d|j}[gcd(i,j)==d]\frac{i*j}{d^2}$
熟悉的，后面改为枚举$d$
$=\sum\limits_{d=1}^kf(d,d)\sum\limits_{i=1}^{\left\lfloor\frac{k}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{k}{d}\right\rfloor}[gcd(i,j)==1]i*j$
开始搞后面那个玩意
$\sum_{i=1}^n\sum_{j=1}^n[gcd(i,j)=1]i*j$
$=2\sum_{i=1}^n\sum_{j=1}^{i}[gcd(i,j)==1]i*j-\sum_{i=1}^ni$
我们知道$n$以内$(n>=2)$和$n$互质的数的和是$\frac{n*φ(n)}{2}$
，然后$i=1$的时候就是$1$
$=\sum_{i=1}^ni^2φ(i)$
回到原来的式子
$=\sum\limits_{d=1}^kf(d,d)\sum\limits_{i=1}^{\left\lfloor\frac{k}{d}\right\rfloor}i^2φ(i)$
可做了，前面那个预处理，后面的分块
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int N=,Sq=,mod=1e9+;
int pri[N],npr[N],phi[N];
int blo[N],stp[Sq],edp[Sq];
int val[N],bsum[N],asum[Sq],func[N];
int a,b,k,n,m,cnt,tot,sqr; long long x,ans;
int exGCD(int a,int b,int &x,int &y)
{
if(!b) {x=,y=; return a;}
int g=exGCD(b,a%b,y,x); y-=a/b*x;
return g;
}
int GCD(int a,int b)
{
return b?GCD(b,a%b):a;
}
int Inv(int nm,int md)
{
int xx,yy;
exGCD(nm,md,xx,yy);
return (xx%md+md)%md;
}
void Prework()
{
phi[]=,npr[]=true;
sqr=sqrt(n)+,stp[cnt=]=;
for(int i=;i<=n;i++)
{
if(!npr[i])
pri[++tot]=i,phi[i]=i-;
for(int j=;j<=tot&&1ll*i*pri[j]<=n;j++)
{
npr[i*pri[j]]=true;
phi[i*pri[j]]=phi[i]*pri[j];
if(i%pri[j]) phi[i*pri[j]]-=phi[i];
else break;
}
}
for(int i=;i<=n;i++)
{
val[i]=1ll*i*i%mod,blo[i]=(i-)/sqr+;
if(i%sqr==) edp[cnt++]=i,stp[cnt]=i+;
func[i]=(func[i-]+1ll*i*i%mod*phi[i]%mod)%mod;
}
(n%sqr)?edp[cnt]=n:cnt--;
for(int i=;i<=cnt;i++)
{
bsum[stp[i]]=val[stp[i]];
for(int j=stp[i]+;j<=edp[i];j++)
bsum[j]=(bsum[j-]+val[j])%mod;
asum[i]=(asum[i-]+bsum[edp[i]])%mod;
}
}
void Change(int a,int b,int x)
{
int g=GCD(a,b),bl=blo[g],v=Inv(1ll*a*b%mod,mod);
val[g]=1ll*g*g%mod*x%mod*v%mod;
bsum[g]=(g==stp[bl])?val[g]:(bsum[g-]+val[g])%mod;
for(int i=g+;i<=edp[bl];i++)
bsum[i]=(bsum[i-]+val[i])%mod;
for(int i=bl;i<=cnt;i++)
asum[i]=(asum[i-]+bsum[edp[i]])%mod;
}
long long Query(int pos)
{
return pos?(asum[blo[pos]-]+bsum[pos])%mod:;
}
int main()
{
scanf("%d%d",&m,&n),Prework();
while(m--)
{
scanf("%d%d%lld%d",&a,&b,&x,&k);
x%=mod,Change(a,b,x),ans=;
for(int i=,j;i<=k;i=j+)
j=k/(k/i),ans+=1ll*(Query(j)-Query(i-)+mod)%mod*func[k/i]%mod,ans%=mod;
printf("%lld\n",ans);
}
return ;
}