根据辗转相减可以知道，相互影响的就是 gcd 相同的那些数。又根据条件 2 可以知道，所有gcd相同的数的比例是不会改变的，又因为最开始 a(x,y)=xy 是一组解，因此每次操作都相当于把所有 gcd 相同的数修改为 a(x,y)=kxy 。因此，对于修改我们只需要维护 f(d) 表示对于所有 gcd(x,y)=d 的 (x,y) 有 a(x,y)=f(d)∗xy 。 接下来考虑询问，也就是求
====∑i=1n∑j=1nf(gcd(i,j))∗ij∑d=1nf(d)∑i=1n∑j=1nij∗[gcd(i,j)=d]∑d=1nf(d)∗d2∑i=1⌊nd⌋∑j=1⌊nd⌋ij∗e(gcd(i,j))∑d=1nf(d)∗d2∑i=1⌊nd⌋∑j=1⌊nd⌋ij∑x|i,x|jμ(x)∑d=1nf(d)∗d2∑x=1⌊nd⌋μ(x)x2∑i=1⌊ndx⌋∑j=1⌊ndx⌋ij 记 s(n)=∑i=1ni=n(n+1)2 g(n)=∑i=1nμ(i)∗i2∗s2(⌊ni⌋) 如果能够预处理出 g ，就可以下底函数分块来求解原问题。但是如果直接依次求是O(nn√)的，无法接受。考虑递推。注意到 ⌊nd⌋ 比 ⌊n−1d⌋ 大 1 ，当且仅当d|n，也就是 n 除以d余数为零。因此可以得到 g(n)=g(n−1)+∑i|nμ(i)∗i2(s2(ni)−s2(ni−1))=g(n−1)+n3∑i|nμ(i)∗1i 其中 h(n)=∑i|nμ(i)∗1i 为积性函数可以线性筛，于是我们就可以 O(n) 预处理出 g 。 还有一个问题，在对答案进行下底函数分块的时候需要维护f(i)∗i2的前缀和，如果直接用树状数组的话，修改是 O(logn) 的，询问是 O(n√logn) 的，无法接受。因此可以分块，牺牲修改的复杂度来降低询问的复杂度，单次修改 O(n√) ，单次询问 O(1) 。这样最终的复杂度是 O(n+mn√) 。 好像有一点卡常，考虑到比较小的 d=gcd(x,y) 被调用的次数比较多，可以在分块的时候维护后缀和，这样修改次数可以减少。 #include<cstdio> #include<cmath> #include<algorithm> using namespace std; #define LL long long const int p=1000000007,maxn=4000010,maxT=2010; LL f[maxn],g[maxn],h[maxn],inv[maxn], sum[maxT],sumin[maxT][maxT]; int prm[maxn],have[maxn],bel[maxn],L[maxT],R[maxT], q,n,T,tot; int gcd(int x,int y) { return y?gcd(y,x%y):x; } void pause() { int x; x=1; } void init() { inv[1]=g[1]=h[1]=1; for (int i=2;i<=n;i++) { inv[i]=(p-inv[p%i])*(p/i)%p; if (!have[i]) { prm[++tot]=i; h[i]=(1-inv[i]+p)%p; } for (int j=1;j<=tot&&(LL)i*prm[j]<=n;j++) { have[i*prm[j]]=1; if (i%prm[j]==0) { h[i*prm[j]]=h[i]; break; } else h[i*prm[j]]=h[i]*h[prm[j]]%p; } g[i]=(g[i-1]+(LL)i*i%p*i%p*h[i])%p; } T=sqrt(n); for (int i=1;i<=T;i++) { L[i]=R[i-1]+1; R[i]=i==T?n:L[i]+T-1; for (int j=L[i];j<=R[i];j++) { f[j]=1; bel[j]=i; sumin[i][j-L[i]+1]=(sumin[i][j-L[i]]+(LL)j*j)%p; } } for (int i=T;i;i--) sum[i]=(sum[i+1]+sumin[i][R[i]-L[i]+1])%p; } LL getsum(int l,int r) { //if (l==10&&r==13) pause(); //LL ret1; int x=bel[l],y=bel[r]; if (x==y) return (sumin[x][r-L[x]+1]-sumin[x][l-L[x]]+p)%p; return (sum[x+1]-sum[y]+p+sumin[y][r-L[y]+1]+sumin[x][R[x]-L[x]+1]-sumin[x][l-L[x]]+p)%p; /*LL ret=0; for (int i=l;i<=r;i++) ret=(ret+f[i]*i%p*i%p)%p; if (ret!=ret1) pause(); return ret;*/ } void update(int d) { for (int j=d;j<=R[bel[d]];j++) sumin[bel[d]][j-L[bel[d]]+1]=(sumin[bel[d]][j-L[bel[d]]]+f[j]*j%p*j)%p; for (int j=bel[d];j;j--) sum[j]=(sum[j+1]+sumin[j][R[j]-L[j]+1])%p; } void solve(int m) { LL ans=0; int u; for (int j=1;j<=m;j=u+1) { u=m/(m/j); ans=(ans+getsum(j,u)*g[m/j])%p; } printf("%lld\n",ans); } int main() { /*freopen("table.in","r",stdin); freopen("table.out","w",stdout);*/ int xx,yy,m,d; LL x; scanf("%d%d",&q,&n); init(); while (q--) { scanf("%d%d%lld%d",&xx,&yy,&x,&m); x%=p; d=gcd(xx,yy); f[d]=(LL)x*inv[xx]%p*inv[yy]%p; update(d); solve(m); } }