题目
题目描述
小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之后的结果。
输入输出样例
输入样例#1：
3 3
1 1 1 2
2 2 4 3
1 2 4 2
输出样例#1：
9
36
14
输入样例#2：
4 125
1 2 4 8
1 3 9 27
1 4 16 64
1 5 25 125
输出样例#2：
2073
316642
12157159
213336861
说明
【输入输出样例 1 说明】
一开始，表格的前3行前3列如图中左边所示。前2次操作后表格没有变化，第3次操作之后的表格如右边所示。
1 2 3 2 4 6
2 4 6 4 4 12
3 6 9 6 12 9
对于 100%的测试点，1 ≤ m ≤ 10^4, 1 ≤ a,b,k ≤ n ≤ 4×10^6, 0 ≤ x < 10^18。
题解
先看第二个式子：$b\times f(a,a+b)=(a+b)\times f(a,b)$
可以化成：$\frac{f(a,a+b)}{a+b}=\frac{f(a,b)}{b}$
即$\frac{f(a,a+b)}{a(a+b)}=\frac{f(a,b)}{ab}$
设$c=a+b$，有$\frac{f(a,c)}{ac}=\frac{f(a,c-a)}{a(c-a)}=\frac{f(a,c-2a)}{a(c-2a)}=\frac{f(a,c\%a)}{a(c\%a)}=\frac{f(c\%a,a)}{(c\%a)a}$
最后有点像求$\gcd$，所以可以直接化成$\frac{f(\gcd(a,c),\gcd(a,c))}{\gcd(a,c)^2}$
所以$\frac{f(a,b)}{ab}=\frac{f(\gcd(a,b),\gcd(a,b))}{\gcd(a,b)^2}$
发现位置$(a,b)$的值可以通过位置$(\gcd(a,b),\gcd(a,b))$求出
则$f(\gcd(a,b),\gcd(a,b))=\frac{f(a,b)\times \gcd(a,b)^2}{ab}$，就可以维护$f(x,x)$的值了，这里的维护先按下不表
题目要求的是$\sum_{i=1}^{k}\sum_{j=1}^{k}f(i,j)$
即$\sum_{i=1}^{k}\sum_{j=1}^{k}\frac{f(\gcd(i,j),\gcd(i,j))\times ij}{\gcd(i,j)^2}$
$=\sum_{d=1}^{k}f(d,d)\sum_{i=1}^{k}\sum_{j=1}^{k}[\gcd(i,j)==d]\frac{ij}{d^2}$
$=\sum_{d=1}^{k}f(d,d)\sum_{i=1}^{\left\lfloor\frac{k}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{k}{d}\right\rfloor}[\gcd(i,j)==1]ij$
后面的部分是个莫比乌斯反演的形式！根据套路可以化成：
$\sum_{d=1}^{k}f_{d,d}\sum_{x=1}^{\left\lfloor\frac{k}{d}\right\rfloor}\mu(x)x^2\sum_{i=1}^{\left\lfloor\frac{k}{dx}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{k}{dx}\right\rfloor}ij$
然后怎么办……
可以整除分块套整除分块，再用一个支持单点修改+区间查询的什么玩意维护一下$f(x,x)$，能得$50$分的好成绩（并不）
当然你可以枚举$T=dx$，然后式子就成了$\sum_{T=1}^{k}\sum_{x\vert T}f(x,x)\mu(\frac{T}{x})(\frac{T}{x})^2$（大概），反正这玩意我不知道怎么快速维护……
再来看看这个式子
$\sum_{d=1}^{k}f(d,d)\sum_{i=1}^{\left\lfloor\frac{k}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{k}{d}\right\rfloor}[\gcd(i,j)==1]ij$
发现$i,j$范围是一样的，则可以化成$\sum_{d=1}^{k}f(d,d)\times 2\sum_{i=1}^{\left\lfloor\frac{k}{d}\right\rfloor}i\sum_{j=1}^{i}[\gcd(i,j)==1]j$
后面部分是求小于$i$与$i$互质的数的和，有公式$\sum_{i=1}^{n}[\gcd(i,n)==1]i=\frac{\phi(n)}{2}n$
因为如果$x$与$n$互质的话，$n-x$与$n$也互质，这样可以两两配对，共$\frac{\phi(n)}{2}$对
这样式子就成了$\sum_{d=1}^{k}f(d,d)\sum_{i=1}^{\left\lfloor\frac{k}{d}\right\rfloor}\phi(i)i^2$
后面的部分可以预处理出来！这样就只需要一个整除分块了
现在考虑一个支持单点修改和维护前缀和数据结构来维护$f(x,x)$
- 我会树状数组！$O(\log n)$的查询，$O(\log n)$的修改，这样一次查询的复杂度就是$O(\sqrt n\times \log n)$，有点卡时间但是树状数组常数小（你懂我意思吧）
- 我会分块！$O(1)$查询，$O(\sqrt n)$修改，这样一次查询复杂度就是$O(\sqrt n)$的，时间很充裕
实测分块比树状数组快$1s$左右，这里提供的是分块代码
代码
# include<iostream>
# include<cstring>
# include<cstdio>
# include<cmath>
# include<algorithm>
using namespace std;
const int MAX=4e6+5,mod=1e9+7,M=2e3+5;
int m,n,tot,k;
int phi[MAX],pm[MAX],val[MAX],sum[MAX],s[MAX],pos[MAX],bl[M],br[M],add[M];
bool use[MAX];
int gcd(int a,int b) {return b?gcd(b,a%b):a;}
void ss()
{
phi[1]=sum[1]=1,val[1]=1,k=sqrt(n),pos[1]=bl[1]=br[1]=1;
for(int i=2,x;i<=n;++i)
{
val[i]=val[i-1]+1ll*i*i%mod,pos[i]=x=(i-1)/k+1,br[x]=i;
if(val[i]>=mod) val[i]-=mod;
if(!bl[x]) bl[x]=i;
if(!use[i]) pm[++tot]=i,phi[i]=i-1;
for(int j=1;j<=tot&&1ll*pm[j]*i<=n;++j)
{
use[pm[j]*i]=1;
if(i%pm[j]==0) {phi[pm[j]*i]=phi[i]*pm[j];break;}
phi[pm[j]*i]=phi[i]*(pm[j]-1);
}
sum[i]=(sum[i-1]+1ll*phi[i]*i%mod*i%mod)%mod+mod;
if(sum[i]>=mod) sum[i]-=mod;
}
}
int read()
{
int x(0);
char ch=getchar();
for(;!isdigit(ch);ch=getchar());
for(;isdigit(ch);x=(1ll*x*10%mod+ch-48)%mod,ch=getchar());
return x;
}
int Pow(int a,int b)
{
int ans=1;
for(;b;b>>=1,a=1ll*a*a%mod)
if(b&1) ans=1ll*ans*a%mod;
return ans;
}
void change(int l,int r,int d)
{
int Pl=pos[l],Pr=pos[r];
if(Pl==Pr)
{
if(l==bl[Pl]&&r==br[Pr])
{
add[Pl]+=d;
if(add[Pl]>=mod) add[Pl]-=mod;
return;
}
for(int i=l;i<=r;++i)
{
val[i]+=d;
if(val[i]>=mod) val[i]-=d;
}
return;
}
int L=Pl+(l!=bl[Pl]),R=Pr-(r!=br[Pr]);
for(int i=L;i<=R;++i)
{
add[i]+=d;
if(add[i]>=mod) add[i]-=mod;
}
if(L!=Pl) for(int i=l;i<=br[Pl];++i)
{
val[i]+=d;
if(val[i]>=mod) val[i]-=mod;
}
if(R!=Pr) for(int i=bl[Pr];i<=r;++i)
{
val[i]+=d;
if(val[i]>=mod) val[i]-=mod;
}
}
int Cal(int x,int a,int b,int d) {return 1ll*x*d%mod*d%mod*Pow(a,mod-2)%mod*Pow(b,mod-2)%mod;}
int Cal_(int x,int d) {return (((x-val[d]+mod)%mod-add[pos[d]]+mod)%mod+(val[d-1]+add[pos[d-1]])%mod)%mod;}
int main()
{
m=read(),n=read(),ss();
for(int i=1,a,b,x,k,d,ans;i<=m;++i)
{
a=read(),b=read(),x=read(),k=read(),d=gcd(a,b),x=Cal(x,a,b,d),change(d,n,Cal_(x,d)),ans=0;
for(int l=1,r,N;l<=k;l=r+1)
{
N=k/l,r=k/N,ans+=1ll*((val[r]+add[pos[r]])%mod-(val[l-1]+add[pos[l-1]])%mod+mod)%mod*sum[N]%mod;
if(ans>=mod) ans-=mod;
}
printf("%d\n",ans);
}
return 0;
}
