推到一半被一个函数卡住真是太屑了的说（菜是原罪）
首先对于(b imes f(a,a+b)=(a+b) imes f(a,b))我们把它化成一个整齐一点的式子：(frac{f(a,a+b)}{a(a+b)}=frac{f(a,b)}{ab})
这样是不是看出点什么了，然后考虑(f)里的变化，是不是相当于一个数减去另一个数
如果我们减到不能减，然后再换过来继续减，直到两数相等，容易发现此时我们相当于在进行辗转相减的过程，实际上容易发现会递归到(gcd(a,b))，因此我们有：
那么就意为着我们只需要考虑(f(d,d))对棋盘的影响即可，考虑对于前(n)行前(n)列的棋盘求和：
其中(S(n)=frac{n(n+1)}{2})，前面的部分可以求出前缀和然后除法分块，那么问题是后面的这部分怎么做（到这里就不会了）
看了题解发现方法确实是妙，我们令(G(n)=sum_{i=1}^n mu(i) imes i^2 imes S(lfloorfrac{n}{i} floor)^2)，发现主要是里面(S(lfloorfrac{n}{i} floor))有个下取整很难搞
考虑运用关于下取整的技巧——差分：
然后就有：
第二步的转化自已化开一下就出来了，最后一步的原因是(varphiast I=idLeftrightarrow varphiast Iast mu=idast muLeftrightarrow varphiast epsilon=idast muLeftrightarrowvarphi=idast mu)，就进行了这一步化(mu)为(varphi)
于是(G(n))就可以(O(n))预处理了
考虑我们还需要修改，显然每次修改都只需要修改一个对应的(f(d,d))即可，而我们同时也要维护(f(d,d))的前缀和
如果直接用树状数组维护单次询问的复杂度达到(O(sqrt nlog n))，这样可能会被卡掉
但是我们套路地发现这个修改的复杂度可以大一些，因为前面不需要乘(sqrt n)的询问次数，因此我们容易想到用分块(O(sqrt n))修改(O(1))询问，总复杂度就是(O(msqrt n))的了
PS：这题不知道为什么BZOJ上一直RE，Luogu和LOJ上都没问题，下载了数据也是对的233
#include<cstdio>
#include<cmath>
#define RI register int
#define CI const int&
using namespace std;
const int N=4e6,mod=1e9+7;
int n,m,a,b,k,prime[N+5],cnt,phi[N+5],g[N+5]; bool vis[N+5]; long long x;
inline int gcd(CI x,CI y)
{
return y?gcd(y,x%y):x;
}
inline int sum(CI x,CI y)
{
int t=x+y; return t>=mod?t-mod:t;
}
inline int sub(CI x,CI y)
{
int t=x-y; return t<0?t+mod:t;
}
#define P(x) prime[x]
inline void init(CI n)
{
RI i,j; for (vis[1]=phi[1]=1,i=2;i<=n;++i)
{
if (!vis[i]) prime[++cnt]=i,phi[i]=i-1;
int tp; for (j=1;j<=cnt&&(tp=i*P(j))<=n;++j)
{
vis[tp]=1; if (i%P(j)) phi[tp]=phi[i]*(P(j)-1);
else { phi[tp]=phi[i]*P(j); break; }
}
}
for (i=1;i<=n;++i) g[i]=sum(g[i-1],1LL*i*i%mod*phi[i]%mod);
}
#undef P
class Value_Blocks
{
private:
static const int S=2005;
int a[N+5],pfx[S][S],tag[S],bel[N+5],L[S],R[S],size;
public:
inline void init(CI n)
{
RI i,j; for (size=(int)sqrt(n),i=1;i<=n;++i) bel[i]=(i-1)/size+1;
for (i=1;i<=bel[n];++i) L[i]=(i-1)*size+1,R[i]=i*size; R[n]=n;
for (i=1;i<=bel[n];++i) for (j=L[i];j<=R[i];++j)
tag[i]=sum(tag[i],a[j]=1LL*j*j%mod),pfx[i][j-L[i]+1]=sum(pfx[i][j-L[i]],a[j]);
for (i=1;i<=bel[n];++i) tag[i]=sum(tag[i],tag[i-1]);
}
inline void modify(CI x,CI y)
{
RI i; int mv=sub(y,a[x]); a[x]=y;
for (i=x;i<=R[bel[x]];++i) pfx[bel[x]][i-L[bel[x]]+1]=sum(pfx[bel[x]][i-L[bel[x]]+1],mv);
for (i=bel[x];i<=bel[n];++i) tag[i]=sum(tag[i],mv);
}
inline int query(CI x)
{
//int ret=0; for (RI i=1;i<=x;++i) ret=sum(ret,a[i]); return ret;
return x?sum(tag[bel[x]-1],pfx[bel[x]][x-L[bel[x]]+1]):0;
}
}B;
int main()
{
for (scanf("%d%d",&m,&n),init(n),B.init(n);m;--m)
{
scanf("%d%d%lld%d",&a,&b,&x,&k); int d=gcd(a,b),ans=0;
B.modify(d,(x/(1LL*(a/d)*(b/d)))%mod); for (RI l=1,r;l<=k;l=r+1)
r=k/(k/l),ans=sum(ans,1LL*sub(B.query(r),B.query(l-1))*g[k/l]%mod);
printf("%d
",ans);
}
return 0;
}