[洛谷P3700] CQOI2017 小Q的表格

问题描述

小Q是个程序员。

作为一个年轻的程序员，小Q总是被老C欺负，老C经常把一些麻烦的任务交给小Q来处理。每当小Q不知道如何解决时，就只好向你求助。

为了完成任务，小Q需要列一个表格，表格有无穷多行，无穷多列，行和列都从1开始标号。为了完成任务，表格里面每个格子都填了一个整数，为了方便描述，小Q把第a行第b列的整数记为f(a,b)。为了完成任务，这个表格要满足一些条件：

(1)对任意的正整数a,b，都要满足f(a,b)=f(b,a)；

(2)对任意的正整数a,b，都要满足b×f(a,a+b)=(a+b)×f(a,b)。

为了完成任务，一开始表格里面的数很有规律，第a行第b列的数恰好等于a×b，显然一开始是满足上述两个条件的。为了完成任务，小Q需要不断的修改表格里面的数，每当修改了一个格子的数之后，为了让表格继续满足上述两个条件，小Q还需要把这次修改能够波及到的全部格子里都改为恰当的数。由于某种神奇的力量驱使，已经确保了每一轮修改之后所有格子里的数仍然都是整数。为了完成任务，小Q还需要随时获取前k行前k列这个有限区域内所有数的和是多少，答案可能比较大，只需要算出答案mod1,000,000,007之后的结果。

输入格式

输入文件第1行包含两个整数m,n，表示共有m次操作，所有操作和查询涉及到的行编号和列编号都不超过n。

接下来m行，每行4个整数a,b,x,k，表示把第a行b列的数改成x，然后把它能够波及到的所有格子全部修改，保证修改之后所有格子的数仍然都是整数，修改完成后计算前k行前k列里所有数的和。

输出格式

输出共m行，每次操作之后输出1行，表示答案mod1,000,000,007之后的结果。

样例输入

3 3
1 1 1 2
2 2 4 3
1 2 4 2

样例输出

9
36
14

数据范围

对于 100%的测试点，1 ≤ m ≤ 10^4, 1 ≤ a,b,k ≤ n ≤ 4×10^6, 0 ≤ x < 10^18。

解析

将第二个条件转化一下可得：\(\frac{f(a,b)}{ab}=\frac{f(a,a+b)}{a(a+b)}\)。这和更相减损法类似，我们可以得到 \(\frac{f(a,b)}{ab}=\frac{f(d,d)}{d^2}\ \ (d=\gcd(a,b))\)。推一推式子：

\[\begin{aligned} \sum_{i=1}^k\sum_{j=1}^kf(i,j)&=\sum_{d=1}^kf(d,d)\sum_{i=1}^{\lfloor \frac{k}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{k}{d}\rfloor}ij[\gcd(i,j)=1]\\ \end{aligned} \]

交换 \(i,j\) 对结果没有影响，我们可以只计算一半，注意减掉重复计算的 \(i=j\) 的情况。我们对右边搞一下：

\[\begin{aligned} \sum_{i=1}^n\sum_{j=1}^nij[\gcd(i,j)=1]&=2\sum_{i=1}^n\sum_{j=1}^iij\sum_{d|i,d|j}\mu(d)-\sum_{i=1}^n i^2[\gcd(i,i)=1]\\ &=2\sum_{i=1}^ni\sum_{d|i} \mu(d)\sum_{j=1}^{i/d}jd-1\\ &=2\sum_{i=1}^ni\sum_{d|i}\mu(d)\frac{\frac{i}{d}\times(\frac{i}{d}+1)}{2}\times d-1\\ &=\sum_{i=1}^ni^2\sum_{d|i}\mu(d)(\frac{i}{d}+1)-1\\ &=\sum_{i=1}^ni^2 (\sum_{d|i}\mu(d)\times \frac{i}{d}+\sum_{d|i}\mu(d))-1\\ &=\sum_{i=1}^ni^2+\phi(i)+\sum_{i=1}^ni^2[i=1]-1\\ &=\sum_{i=1}^n i^2+\phi(i) \end{aligned} \]

那么原式等于

\[\sum_{d=1}^kf(d,d)\sum_{i=1}^{\lfloor \frac{k}{d}\rfloor} i^2+\phi(i) \]

后面可以数论分块搞。而 \(f\) 的修改可以通过分块维护前缀和来维护。具体的，每一块的最后一个点维护整个序列的前缀和，其余点只维护块内前缀和。

代码

#include <iostream>
#include <cstdio>
#include <cmath>
#define N 4000002
using namespace std;
const int mod=1000000007;
int n,m,i,j,prime[N],cnt,b[N],gap;
long long phi[N],sum[N],s[N];
bool vis[N];
long long read()
{
	char c=getchar();
	long long w=0;
	while(c<'0'||c>'9') c=getchar();
	while(c<='9'&&c>='0'){
		w=w*10+c-'0';
		c=getchar();
	}
	return w;
}
int gcd(int a,int b)
{
	if(!b) return a;
	return gcd(b,a%b);
}
int poww(int a,int b)
{
	long long ans=1,base=a;
	while(b){
		if(b&1) ans=ans*base%mod;
		base=base*base%mod;
		b>>=1;
	}
	return ans;
}
void change(int x,int v)
{
	for(int i=x;b[i]==b[x];i++) sum[i]=(sum[i]+v)%mod;
	for(int i=(b[x]+1)*gap;i<=n;i+=gap) sum[i]=(sum[i]+v)%mod;
	if(n%gap!=0&&b[n]!=b[x]) sum[n]=(sum[n]+v)%mod;
}
int ask(int x)
{
	if(x%gap==0||x==n) return sum[x];
	return (sum[x]+sum[(b[x]-1)*gap])%mod;
}
signed main()
{
	m=read();n=read();
	phi[1]=1;
	for(i=2;i<=n;i++){
		if(!vis[i]){
			prime[++cnt]=i;
			phi[i]=i-1;
		}
		for(j=1;j<=cnt;j++){
			if(i*prime[j]>n) break;
			vis[i*prime[j]]=1;
			if(i%prime[j]) phi[i*prime[j]]=phi[i]*(prime[j]-1)%mod;
			else{
				phi[i*prime[j]]=prime[j]*phi[i]%mod;
				break;
			}
		}
	}
	for(i=1;i<=n;i++) phi[i]=phi[i]*i%mod*i%mod;
	for(i=1;i<=n;i++) phi[i]=(phi[i-1]+phi[i])%mod;
	gap=sqrt(1.0*n);
	for(i=1;i<=n;i++) b[i]=(i-1)/gap+1;
	long long tmp=0;
	for(i=1;i<=n;i++){
		tmp=(tmp+1LL*i*i%mod)%mod;s[i]=1LL*i*i%mod;
		if(i%gap==0||i==n) sum[i]=tmp;
		else if((i-1)%gap!=0) sum[i]=(sum[i-1]+1LL*i*i%mod)%mod;
		else sum[i]=1LL*i*i%mod;
	}
	for(i=1;i<=m;i++){
		int a=read(),b=read();
		long long x=read(),ans=0;
		int k=read();
		int d=gcd(a,b);
		long long val=1LL*x%mod*d%mod*d%mod*poww(1LL*a*b%mod,mod-2)%mod;
		change(d,(val-s[d]+mod)%mod);s[d]=val;
		for(int l=1,r;l<=k;l=r+1){
			r=k/(k/l);
			ans=(ans+1LL*(ask(r)-ask(l-1)+mod)%mod*phi[k/l]%mod)%mod;
		}
		printf("%lld\n",ans);
	}
	return 0;
}