洛谷 P3700

又是一道需要一些观察的数论 hot tea……

注意到题目中 (b·f(a,a+b)=(a+b)·f(a,b)) 这个柿子长得有点像求解 (gcd) 的辗转相除法，因此考虑从这方面入手解决这道题，不难发现对于两个数 ((a,b))，通过辗转相除法咱们总可以得到 (f(gcd(a,b),gcd(a,b))) 处。那么我们考虑归纳求出 (f(a,b)) 是个什么东西，记 (d=gcd(a,b))，由于我们不管怎么辗转相除都只能除到 (f(d,d)) 处，因此我们保持 (f(d,d)) 不变好了，那么 (f(d,2d)=2f(d,d),f(d,3d)=3f(d,d),cdots,f(d,kd)=kf(d,d))，而 (f((k+1)d,kd)=(k+1)f(d,kd)=k(k+1)f(d,d))，(f((2k+1)d,kd)=dfrac{2k+1}{k+1}·f((k+1)d,kd)=(2k+1)kf(d,d))，如此归纳下去不难得出：

[f(ad,bd)=abf(d,d)(aperp b) ]

也就是说我们只用维护所有 (f(d,d)) 的变化即可。

注意到当我们改变某个 (f(a,b)=x) 时候，由于 (f(a,b)=dfrac{ab}{gcd(a,b)^2}f(gcd(a,b),gcd(a,b)))，因此改变 (f(a,b)) 会影响 (f(gcd(a,b),gcd(a,b)))，也进而影响其他满足 (gcd(x,y)=gcd(a,b)) 的 (f(x,y)) 的值。而对于 (d egcd(a,b))，改变 (f(a,b)) 显然是不会影响 (f(d,d)) 的值的，因此改变 (f(a,b)) 只会影响一个 (f(d,d)) 的值。

接下来考虑求解答案。

[egin{aligned} ans&=sumlimits_{i=1}^ksumlimits_{j=1}^kdfrac{ij}{gcd(i,j)^2}f(gcd(i,j),gcd(i,j))\ &=sumlimits_{d=1}^kf(d,d)dfrac{1}{d^2}sumlimits_{i=1}^{n/d}sumlimits_{j=1}^{n/d}id·jd[gcd(i,j)=1]\ &=sumlimits_{d=1}^kf(d,d)sumlimits_{i=1}^{n/d}sumlimits_{j=1}^{n/d}ij[gcd(i,j)=1]\ end{aligned} ]

设

[g(x)=sumlimits_{i=1}^xsumlimits_{j=1}^xij[gcd(i,j)=1] ]

那么

[ans=sumlimits_{d=1}^kf(d,d)g(lfloordfrac{k}{d} floor) ]

我们先不考虑 (g(x)) 怎么求，假设我们已经求得了所有 (g(x))，那么显然可以对 (lfloordfrac{k}{d} floor) 整除分块，那么我们需要求出 (f(d,d)) 的前缀和，又因为要修改，因此我们需要一个支持修改、查询前缀和的数据结构。有人肯定会想树状数组，不过对于此题而言，由于修改次数为 (m)，查询次数为 (msqrt{n})，因此直接 BIT 复杂度为 (msqrt{n}log n) 无法通过，因此考虑根号平衡的思想，众所周知分块可以实现 (mathcal O(sqrt{n})) 修改 (mathcal O(1)) 查询前缀和，使用这样的数据结构维护一下即可 (mathcal O(msqrt{n})) 查询。

最后考虑 (g(x)) 怎么求，直接整除分块是 (nsqrt{n}) 的，无法通过。因此考虑不反演。首先证明一个 Lemma：

[sumlimits_{i=1}^xi[gcd(i,x)=1]=dfrac{x(varphi(x)+epsilon(x))}{2} ]

证明大概就 (forall ile x)，若 (iperp x)，那么 ((x-i)perp x)，那么我们就对于所有 (ile x,iperp x)，将 (i) 与 (x-i) 配对形成一个 (x) 即可，这样共可以配成 (dfrac{varphi(x)}{2}) 对，然后对 (x=1) 和 (x=2) 特判一下即可。

又这个引理可知

[egin{aligned} g(x)&=g(x-1)+2xsumlimits_{i=1}^xi[gcd(i,x)=1]-epsilon(x)\ &=g(x-1)+x^2(varphi(x)+epsilon(x))-epsilon(x)\ &=g(x-1)+x^2varphi(x) end{aligned} ]

因此

[g(x)=sumlimits_{i=1}^xi^2varphi(i) ]

线筛算算即可。

const int MAXN=4e6;
const int BLK=2000;
const int MOD=1e9+7;
void pls(int &x,int v){((x+=v)>=MOD)&&(x-=MOD);}
int n,qu,phi[MAXN+5],pr[MAXN/8+5],prcnt=0,s[MAXN+5];
bitset<MAXN+5> vis;
void sieve(int n){
	phi[1]=1;
	for(int i=2;i<=n;i++){
		if(!vis[i]) phi[i]=i-1,pr[++prcnt]=i;
		for(int j=1;j<=prcnt&&pr[j]*i<=n;j++){
			vis[pr[j]*i]=1;
			if(i%pr[j]==0){phi[i*pr[j]]=phi[i]*pr[j];break;}
			else phi[i*pr[j]]=phi[i]*phi[pr[j]];
		}
	}
	for(int i=1;i<=n;i++) s[i]=(s[i-1]+1ll*i*i%MOD*phi[i])%MOD;
}
int L[BLK+5],R[BLK+5],bel[MAXN+5],blk_cnt,blk_sz;
int sum_blk[BLK+5],sum_tot[MAXN+5];
void add(int x,int v){
	for(int i=bel[x];i<=blk_cnt;i++) pls(sum_blk[i],v);
	for(int i=x;i<=R[bel[x]];i++) pls(sum_tot[i],v);
}
int ask(int x){if(!x) return 0;return (sum_blk[bel[x]-1]+sum_tot[x])%MOD;}
int f[MAXN+5];
int main(){
	scanf("%d%d",&qu,&n);sieve(n);
//	for(int i=1;i<=n;i++) printf("%d
",s[i]);
	blk_sz=(int)sqrt(n);blk_cnt=(n-1)/blk_sz+1;
	int sum=0;
	for(int i=1;i<=blk_cnt;i++){
		L[i]=(i-1)*blk_sz+1;R[i]=min(i*blk_sz,n);int sum0=0;
		for(int j=L[i];j<=R[i];j++){
			bel[j]=i;sum0=(sum0+1ll*j*j)%MOD;
			sum_tot[j]=sum0;f[j]=1ll*j*j%MOD;
//			printf("%d %d
",j,sum0);
		} pls(sum,sum0);sum_blk[i]=sum;
//		printf("%d %d
",i,sum);
	}
	while(qu--){
		int a,b,k;ll v;scanf("%d%d%lld%d",&a,&b,&v,&k);
		int d=__gcd(a,b);add(d,(MOD-f[d])%MOD);
		f[d]=(v/(1ll*a*b/d/d))%MOD;add(d,f[d]);
		int res=0;
		for(int l=1,r;l<=k;l=r+1){
			r=k/(k/l);
//			printf("%d %d %d
",l,r,(ask(r)-ask(l-1)+MOD)%MOD);
			res=(res+1ll*(ask(r)-ask(l-1)+MOD)*s[k/l])%MOD;
		} printf("%d
",res);
	}
	return 0;
}