題目：https://www.lydsy.com/JudgeOnline/problem.php?id=4815

大概就是推式子的時候注意有兩個邊界都是 n ，考慮變成 2*... 之類的。

分塊維護 f[ ] 的字首和。很好的思路是修改一個位置後字首和陣列需要區間加，整塊地打上加法標記就行了。

自己本來想維護整塊之間的字首和，還有塊內的字首和；卻WA得不行。之後再探究為什麼WA吧。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define ll long long
using
 
 namespace std;
const int N=4e6+5,M=2005,mod=1e9+7;
int n,g[N],phi[N],f[N],pri[N];bool vis[N];
int base,bh[N],s[M],si[M],fl[N];
void upd(int &x){while(x>=mod)x-=mod;while(x<0)x+=mod;}
int gcd(int a,int b){return b?gcd(b,a%b):a;}
int pw(int x,int k)
{int ret=1;while(k){if(k&1)ret=(ll)ret*x%mod;x=(ll)x*x%mod;k>>=1
 
;}return ret;}
void init()
{
  phi[1]=g[1]=1; int cnt=0;
  for(int i=2;i<=n;i++)
    {
      if(!vis[i])pri[++cnt]=i,phi[i]=i-1;
      for(int j=1;j<=cnt&&(ll)i*pri[j]<=n;j++)
    {
      vis[i*pri[j]]=1;
      if(i%pri[j]==0){phi[i*pri[j]]=(ll)phi[i]*pri[j]%mod;break;}
      else phi[i*pri[j]]=(ll)phi[i]*phi[pri[j]]%mod;
    }
      g[i]
 
=(g[i-1]+(ll)i*i%mod*phi[i])%mod;//presum!!
    }
  base=sqrt(n);
  for(int i=1;i<=n;i++)f[i]=(ll)i*i%mod,fl[i]=(fl[i-1]+f[i])%mod;
  for(int i=1,j=1,k=1;i<=n;i++,k++)
    {
      bh[i]=j;if(k==base)k=0,j++;
    }
  /*
  for(int i=1,j=1,k=base;i<=n;i++)
    {
      f[i]=(ll)i*i%mod;bh[i]=j;
      si[j]+=f[i]; upd(si[j]); fl[i]=si[j];
      if(i==k)s[j]=s[j-1]+si[j],j++,k+=base;
    }
  */
}
int calc(int x){int ret=fl[x]+s[bh[x]];upd(ret);return ret;}
int main()
{
  int T;scanf("%d%d",&T,&n);init();
  int x,y,tn; ll w;
  while(T--)
    {
      scanf("%d%d%lld%d",&x,&y,&w,&tn);//w not %mod!!!
      int u=gcd(x,y),d=bh[u];
      int tf=w/(x/u)/(y/u)%mod;//not inv
 
      /*//also ok
      int chg=tf+calc(u-1)-calc(u);upd(chg);
      for(int i=u;bh[i]==bh[u];i++)
    fl[i]+=chg,upd(fl[i]);
      for(int i=bh[u]+1;i<=bh[n];i++)//n not tn!!!
    s[i]+=chg,upd(s[i]);
      */
      int pl=tf-f[u];upd(pl);f[u]=tf;
      for(int i=u;bh[i]==bh[u];i++)
    fl[i]+=pl,upd(fl[i]);
      for(int i=bh[u]+1;i<=bh[n];i++)
    s[i]+=pl,upd(s[i]);
      /*
      si[d]=si[d]-f[u]+tf;upd(si[d]);
      for(int i=d;i<=bh[n];i++)
    s[i]=s[i-1]+si[i],upd(s[i]);
      f[u]=tf;
 
      fl[u]=(bh[u-1]==bh[u]?fl[u-1]:0)+f[u];
      for(int i=u+1,j=d*base;i<=j;i++)
    fl[i]=fl[i-1]+f[i],upd(fl[i]);
      */
      int ans=0;
      for(int i=1,j;i<=tn;i=j+1)
    {
      int d=tn/i,sm=0; j=tn/d;
      /*
      if(bh[j]-bh[i]<=1)
        for(int l=i;l<=j;l++)
          sm+=f[l],upd(sm);
      else
        {
          sm=s[bh[j]-1]-s[bh[i]-1]+fl[j];
          upd(sm);
          if(bh[i-1]==bh[i])sm-=fl[i-1],upd(sm);
        }
      ans=(ans+(ll)sm*g[d])%mod;
      */
      ans=(ans+(ll)(calc(j)-calc(i-1))*g[d])%mod;upd(ans);
    }
      printf("%d\n",ans);
    }
  return 0;
}