题目链接：http://www.lydsy.com/JudgeOnline/problem.php?id=4815
题目中所给条件中的$(a,a+b)$和$(a,b)$的关系很瞩目。
然后大家都知道$(a,b)=(a,a-b)=(a,a+b)$，于是观察（猜）一下这个表格与gcd的关系。
可以发现每次修改$(a,b)$会影响到所有$(i,j)=(a,b)$的点，并且关系为$$f(i,j)=frac{i}{a}*frac{j}{b}*f(a,b)$$
所以只需要知道$f(d,d)$的值记为$f(d)$，就能推出其他的值。
然后慢慢推推推大概可以推到这一步$$ans=sum_{d=1}^nf(d)sum_{i=1}^{frac{n}{d}}sum_{j=1}^{frac{n}{d}}(i,j)[(i,j)==1]$$
可以发现这个式子中$i$和$j$是对称的$$S(frac{n}{d})=sum_{i=1}^{frac{n}{d}}sum_{j=1}^{frac{n}{d}}(i,j)[(i,j)==1]$$
不妨先设$i>j$，于是我们有$$S′(n)=sum_{i=1}^nfrac{φ(i)*i^{2}}{2}$$
由于$i$与$j$对称，所以有$$S(n)=2*S′(n)=sum_{i=1}^nφ(i)*i^{2}$$
所以最终的答案就变成了$$ans=sum_{d=1}^nf(d)S(frac{n}{d})$$
我们记录$f$的前缀和，并且分块维护这个数列，而$S$很明显是可以预处理出来的。
询问了$m$次，于是总体复杂度应该是$O(msqrt{n})$
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #include<cmath> 5 using namespace std; 6 typedef long long ll; 7 const int mod=1e9+7; 8 int inline readint(){ 9 int Num;char ch; 10 while((ch=getchar())<'0'||ch>'9');Num=ch-'0'; 11 while((ch=getchar())>='0'&&ch<='9') Num=Num*10+ch-'0'; 12 return Num; 13 } 14 ll inline readll(){ 15 ll Num;char ch; 16 while((ch=getchar())<'0'||ch>'9');Num=ch-'0'; 17 while((ch=getchar())>='0'&&ch<='9') Num=Num*10+ch-'0'; 18 return Num; 19 } 20 void outint(int x){ 21 if(x>=10) outint(x/10); 22 putchar(x%10+'0'); 23 } 24 int inline gcd(int x,int y){ 25 return !y?x:gcd(y,x%y); 26 } 27 int n,m; 28 int phi[4000010],p[4000010],cnt=0; 29 int la,blk,add[2010]; 30 int f[4000010]; 31 bool vis[4000010]; 32 void sieve(int n){ 33 for(int i=2;i<=n;i++){ 34 if(!vis[i]){ 35 p[++cnt]=i; 36 phi[i]=i-1; 37 } 38 for(int j=1;p[j]*i<=n;j++){ 39 vis[p[j]*i]=true; 40 if(i%p[j]==0){ 41 phi[p[j]*i]=phi[i]*p[j]; 42 break; 43 } 44 phi[p[j]*i]=phi[i]*(p[j]-1); 45 } 46 phi[i]=(1LL*i*i%mod*phi[i]+phi[i-1])%mod; 47 f[i]=(1LL*i*i+f[i-1])%mod; 48 } 49 } 50 void modify(int x,int ad){ 51 int l=(x-1)/blk+1, 52 r=min(n,l*blk); 53 for(int i=l+1;i<=la;i++) add[i]=(add[i]+ad)%mod; 54 for(int i=x;i<=r;i++) f[i]=(f[i]+ad)%mod; 55 } 56 int inline qry(int x){ 57 return x?(f[x]+add[(x-1)/blk+1])%mod:0; 58 } 59 int main(){ 60 m=readint(); 61 n=readint(); 62 f[1]=phi[1]=1; 63 blk=(int)sqrt(n); 64 la=(n-1)/blk+1; 65 sieve(n); 66 for(int i=1;i<=m;i++){ 67 int a=readint(), 68 b=readint(), 69 g=gcd(a,b), 70 ans=0; 71 ll x=readll(); 72 int k=readint(); 73 x=x/(1LL*(a/g)*(b/g))%mod; 74 modify(g,((x-qry(g)+qry(g-1))%mod+mod)%mod); 75 for(int j=1,now;j<=k;j=now+1){ 76 now=k/(k/j); 77 ans=(ans+1LL*(qry(now)-qry(j-1)+mod)%mod*phi[k/j])%mod; 78 } 79 outint(ans); 80 putchar(' '); 81 } 82 return 0; 83 }