Luogu

LOJ

分析

可以发现这两条规则长得很像更相减损术，推一推可以得到
$$
f(a,b)=\frac{a\times b}{\gcd(a,b)^2}f(\gcd(a,b),\gcd(a,b))
$$
于是
$$
\begin{aligned}
&\sum_{i=1}^n\sum_{j=1}^nf(i,j)\\
=&\sum_{i=1}^n\sum_{j=1}^n\frac{ij}{\gcd(i,j)^2}f(\gcd(i,j),\gcd(i,j))\\
=&\sum_{d=1}^nf(d,d)\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}ij[\gcd(i,j)=1]\\
=&\sum_{d=1}^nf(d,d)\sum_{i=1}^{n/d}i^2\varphi(i)
\end{aligned}
$$
最后那一步可以考虑实际意义，证明就不写了（

如果没有修改的话就可以预处理 $i^2\varphi(i)$ 的前缀和然后数论分块计算了。然而现在还有若干次对于 $f(d,d)$ 的修改，可以想到用一个数据结构来维护。

思考一下，我们在数论分块时需要进行 $\mathcal{O}(Q)$ 次修改和求 $\mathcal{O}(Q\sqrt{n})$ 次前缀和，因此用值域分块来做到 $\mathcal{O}(\sqrt{n})$ 修改和 $\mathcal{O}(1)$ 求前缀和即可。

代码

// ====================================
//   author: M_sea
//   website: https://m-sea-blog.com/
// ====================================
#include <bits/stdc++.h>
#define file(x) freopen(#x".in","r",stdin); freopen(#x".out","w",stdout)
using namespace std;
typedef long long ll;

ll read() {
    ll X=0; char c=getchar();
    while (c<'0'||c>'9') c=getchar();
    while (c>='0'&&c<='9') X=X*10+c-'0',c=getchar();
    return X;
}

const int N=4000000+10,M=2000+10;
const int mod=1e9+7;

int n,Q,f[N];

int p[N],v[N],cnt=0,phi[N],s[N];
void sieve(int n) {
    phi[1]=1;
    for (int i=2;i<=n;++i) {
        if (!v[i]) p[++cnt]=i,phi[i]=i-1;
        for (int j=1;j<=cnt&&i*p[j]<=n;++j) {
            v[i*p[j]]=1;
            if (i%p[j]) phi[i*p[j]]=phi[i]*(p[j]-1);
            else { phi[i*p[j]]=phi[i]*p[j]; break; }
        }
    }
    for (int i=1;i<=n;++i) s[i]=1ll*i*i%mod*phi[i]%mod;
    for (int i=2;i<=n;++i) s[i]=(s[i-1]+s[i])%mod;
}

namespace T {
    int B,bl[N],L[M],R[M],sb[M],s[N];
    void build() {
        B=sqrt(n);
        for (int i=1;i<=n;++i) bl[i]=(i-1)/B+1;
        for (int i=1;i<=bl[n];++i) L[i]=(i-1)*B+1,R[i]=min(i*B,n);
        for (int i=1;i<=n;++i) s[i]=bl[i]==bl[i-1]?(s[i-1]+f[i])%mod:f[i];
        for (int i=1;i<=bl[n];++i) sb[i]=(sb[i-1]+s[R[i]])%mod;
    }
    void modify(int p,int w) {
        for (int i=p;i<=R[bl[p]];++i) s[i]=(s[i]+w)%mod;
        for (int i=bl[p];i<=bl[n];++i) sb[i]=(sb[i]+w)%mod;
    }
    int sum(int p) { return (sb[bl[p]-1]+s[p])%mod; }
}

int main() {
    Q=read(),n=read(); sieve(n);
    for (int i=1;i<=n;++i) f[i]=1ll*i*i%mod;
    T::build();
    while (Q--) {
        int a=read(),b=read(); ll x=read(); int k=read();
        int d=__gcd(a,b);
        T::modify(d,mod-f[d]),f[d]=(x/(a/d)/(b/d))%mod,T::modify(d,f[d]);
        int ans=0;
        for (int l=1,r;l<=k;l=r+1) {
            r=k/(k/l);
            ans=(ans+1ll*(T::sum(r)-T::sum(l-1)+mod)*s[k/l])%mod;
        }
        printf("%d\n",ans);
    }
    return 0;
}