题目链接

做法

首先打表找规律可以得知
令$d=gcd(x,y)$

$f_{x,y}=f_{d,d}*\frac{x*y}{d^2}$

那么ans就可以换一种表现形式

$ans=\sum_{i=1}^{k}\sum_{j=1}^{k}f_{i,j}\\ ans=\sum_{i=1}^{k}\sum_{j=1}^{k}f_{d,d}*\frac{x*y}{d^2}$

改变枚举顺序，先枚举d

$ans=\sum_{d=1}^{k}*f_{d,d}*\sum_{i=1}^{\lfloor\frac{k}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{k}{d}\rfloor}i*j*[gcd(i,j)=1]$

我们令

$A(x)=\sum_{i=1}^{x}\sum_{j=1}^{x}i*j*[gcd(i,j)=1]$

考虑$A(x)$比$A(x-1)$多了什么
多的这部分为

$\sum_{i=1}^{x}i*x*[gcd(i,x)==1]*2-x*x*[gcd(x,x)=1]$

显然$A(1)=1$
所以我们这里只需要考虑$x>1$的情况，后面减的那项为0
而由于$i<x$且$gcd(i,x)=1$，所以$gcd(i,x-i)=1$也成立。
这也就意味着所有与$x$互质且小于$x$的数都可以两两组成和为$x$的组。
所以

$A(x)=A(x-1)+\phi(x)*x^2$

也就是说

$A(x)=\sum_{i=1}^{x}\phi(i)*i^2$

我们再回到开头$ans$的式子

$ans=\sum_{d=1}^{k}*f_{d,d}*A(\lfloor\frac{k}{d}\rfloor)$

$\lfloor\frac{k}{d}\rfloor$只有$\sqrt{k}=2000$种取值，我们数论分块计算。
而对于$f$就变成了一个单点修改，区间求值的问题。
再注意到实际上修改只用进行$10^4$次，而查询需要$2000*10^4$次。
我们就将序列分块，维护块的前缀和，再维护每个块中的前缀和，就可以做到$O(\sqrt{n})$修改，$O(1)$查询了。

代码

#include<cstdio>
#include<algorithm>
#include<cctype>
#include<cstring>
#include<iostream>
#include<cmath>
#define LL long long
#define N (4000005)
using namespace std;
int m,n,pri,ans,x,y,k,g,unit,tot,r;
LL d;
int A[N],prime[N],phi[N],f[N],be[N],sum[2005],pr[2005][2005];
bool nprime[N];
const int P=1000000007; 
template <typename T> void read(T&t) {
    t=0;
    bool fl=true;
    char p=getchar();
    while (!isdigit(p)) {
        if (p=='-') fl=false;
        p=getchar();
    }
    do {
        (t*=10)+=p-48;p=getchar();
    }while (isdigit(p));
    if (!fl) t=-t;
}
inline int Inc(int a,int b){
    return (a+b>=P)?(a+b-P):(a+b);
}
inline int gcd(int a,int b){
    while (a%b!=0){
        int tmp=b;
        b=a%b;
        a=tmp;
    }
    return b;
}
void change(int x,int d){
    d=Inc(d,P);
    for (int i=be[x];i<=tot;i++) sum[i]=Inc(sum[i],d);
    for (int i=x-(be[x]-1)*unit;i<=unit;i++) pr[be[x]][i]=Inc(pr[be[x]][i],d);
}
int query(int x){
    return Inc(sum[be[x]-1],pr[be[x]][x-(be[x]-1)*unit]);
}
signed main(){
    read(m),read(n);
    phi[1]=1;
    for (int i=2;i<=n;i++){
        if (!nprime[i]) prime[++pri]=i,phi[i]=i-1;
        for (int j=1;j<=pri&&i*prime[j]<=n;j++){
            nprime[i*prime[j]]=1;
            if (i%prime[j]==0){
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            else phi[i*prime[j]]=phi[i]*(prime[j]-1);
        }
    }
    for (int i=1;i<=n;i++) A[i]=Inc(A[i-1],1ll*phi[i]*i%P*i%P);
    unit=sqrt(n);
    for (int i=1;i<=n;i++){
        be[i]=(i-1)/unit+1;
    }
    tot=be[n];
    for (int i=1;i<=n;i++){
        f[i]=1ll*i*i%P;
        if (be[i]!=be[i-1]) sum[be[i]]=sum[be[i-1]];
        sum[be[i]]=Inc(sum[be[i]],f[i]);
        int num=i-(be[i]-1)*unit;
        pr[be[i]][num]=Inc(pr[be[i]][num-1],f[i]);
    }
    for (int i=1;i<=m;i++){
        read(x),read(y),read(d),read(k);
        g=gcd(x,y);
        d=d/(x/g)/(y/g)%P;
        change(g,d-f[g]);
        f[g]=d;
        ans=0;
        for (int l=1;l<=k;l=r+1){
            r=k/(k/l);
            ans=Inc(ans,1ll*Inc(query(r)-query(l-1),P)*A[k/l]%P);
        } 
        printf("%d\n",ans);
    }
    return 0;
}