【BZOJ4815】[CQOI2017]小Q的表格（莫比乌斯反演，分块）
题面
题解
神仙题啊。
首先(f(a,b)=f(b,a))告诉我们矩阵只要算一半就好了。
接下来是(b*f(a,a+b)=(a+b)*f(a,b))
这个式子怎么看呢？
[egin{aligned}b*f(a,a+b)&=(a+b)*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)}{a*b}end{aligned}
]
这个式子很明显类似于辗转相减的式子，如果我们令(c=(a+b))，那么等式就可以写成
[frac{f(a,c)}{a*c}=frac{f(a,c-a)}{a*(c-a)}
]
那么进一步，意味着我们可以写成辗转相除的式子。
[frac{f(a,b)}{a*b}=frac{f(a,a\%b)}{a*(a\%b)}
]
所以，类似于求(gcd)的终止状态，我们可以利用这个写出一个等式：
[frac{f(a,b)}{a*b}=frac{f(gcd(a,b),gcd(a,b)}{gcd^2(a,b)}
]
为了方便后面直接令(d=gcd(a,b))，即等式为
[f(a,b)=frac{a*b}{d^2}f(d,d)
]
那么，不难发现，单次的修改只会对于(d)相等的位置产生影响。
考虑(ans)是个啥玩意
[egin{aligned}ans&=sum_{i=1}^ksum_{j=1}^kf(i,j)\&=sum_{d=1}^kf(d,d)sum_{i=1}^ksum_{j=1}^kfrac{i*j}{d^2}[gcd(i,j)=d]\&=sum_{d=1}^kf(d,d)sum_{i=1}^{k/d}sum_{j=1}^{k/d}ij[gcd(i,j)=1]end{aligned}
]
后面一半的式子可以用莫比乌斯反演直接计算值。
当然了，也可以这样子推:
[egin{aligned}sum_{i=1}^nsum_{j=1}^nij[gcd(i,j)=1]&=2*sum_{i=1}^nisum_{j=1}^ij[gcd(i,j)=1]\&=2*sum_{i=1}^nfrac{ivarphi(i)}{2}\&=sum_{i=1}^ni^2varphi(i)end{aligned}
]
后面那一步化简简单证明一下，假设(x)与(n)互质，那么(n-x)也与(n)互质，因此所有与(n)互质的数的和是两两配对的。
也就是后面这一部分是可以提前预处理出来的，因为是(k/d)，也就是等价于可以数论分块，所以我们现在唯一要做的就是维护(f(i,i))的前缀和。
而修改操作显然是把当前位置会影响的位置的值全部对应的扩倍。也就是会影响到(f(gcd,gcd))这个位置的值。我们要找个方法快速计算(f(i,i))的前缀和。
考虑数据范围，操作次数很少，但是操作的范围很大。单次询问的时候我们有一个(sqrt k)的数论分块的复杂度，也就是(2*10^3)，操作次数有(10^4)，所以我们显然不能带(log)的计算前缀和。那么考虑分块，将计算前缀和的复杂度将至(O(1))，而修改复杂度变为(O(sqrt n))。
#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<vector>
using namespace std;
#define ll long long
#define MOD 1000000007
#define MAX 4000400
#define BLK 2200
inline ll read()
{
ll x=0;bool t=false;char ch=getchar();
while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
if(ch=='-')t=true,ch=getchar();
while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
return t?-x:x;
}
int fpow(int a,int b)
{
int s=1;
while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}
return s;
}
int phi[MAX],pri[MAX],tt,s[MAX];
bool zs[MAX];
void Pre(int n)
{
phi[1]=1;
for(int i=2;i<=n;++i)
{
if(!zs[i])pri[++tt]=i,phi[i]=i-1;
for(int j=1;j<=tt&&i*pri[j]<=n;++j)
{
zs[i*pri[j]]=true;
if(i%pri[j])phi[i*pri[j]]=phi[i]*(pri[j]-1);
else{phi[i*pri[j]]=phi[i]*pri[j];break;}
}
}
for(int i=1;i<=n;++i)s[i]=1ll*i*i%MOD*phi[i]%MOD;
for(int i=1;i<=n;++i)s[i]=(s[i]+s[i-1])%MOD;
}
int pre[BLK],ps[BLK],bs[BLK][BLK],val[MAX];
int b[MAX],tot,blk;
int n,m;
int Query(int x){if(!x)return 0;return (pre[b[x]-1]+bs[b[x]][x-blk*(b[x]-1)])%MOD;}
void Modify(int x,int w)
{
val[x]=w;int p=x-blk*(b[x]-1);
for(int i=x;b[i]==b[x];++i,++p)
bs[b[x]][p]=(bs[b[x]][p-1]+val[i])%MOD;
ps[b[x]]=bs[b[x]][p-1];
for(int i=1;i<=tot;++i)pre[i]=(pre[i-1]+ps[i])%MOD;
}
int Calc(int k)
{
int ret=0;
for(int i=1,j;i<=k;i=j+1)
{
j=k/(k/i);
ret=(ret+1ll*(Query(j)-Query(i-1)+MOD)*s[k/i])%MOD;
}
return ret;
}
int main()
{
m=read();n=read();Pre(n);blk=sqrt(n);
for(int i=1;i<=n;++i)b[i]=(i-1)/blk+1;
tot=b[n];
for(int i=1;i<=n;++i)val[i]=1ll*i*i%MOD;
for(int i=1;i<=tot;++i)
for(int j=1,a=(i-1)*blk+1;j<=blk&&a<=n;++j,++a)
ps[i]=(ps[i]+val[a])%MOD,bs[i][j]=(bs[i][j-1]+val[a])%MOD;
for(int i=1;i<=tot;++i)pre[i]=(pre[i-1]+ps[i])%MOD;
while(m--)
{
int a=read(),b=read(),x=read()%MOD,k=read();
int d=__gcd(a,b);
Modify(d,1ll*x*d%MOD*d%MOD*fpow(1ll*a*b%MOD,MOD-2)%MOD);
printf("%d
",Calc(k));
}
return 0;
}