【BZOJ4815】[CQOI2017] 小Q的表格（莫比乌斯反演+分块）
大致题意： 给定一个\(n\times n\)的表格，初始状态下\((a,b)\)格子的值为\(a\times b\)。每次修改一个位置上的值，并要求你通过修改其他位置上的值使得表格满足\(f(a,b)=f(b,a)\)且\(b\times f(a,a+b)=(a+b)\times f(a,b)\)。每次修改后询问前\(k\)行\(k\)列内数字的和。
前言
我还是太弱了，式子推到一半就推不下去，最后还是只能去跪膜闪指导的博客。
这题还真是挺神奇的啊。
修改
首先，经过简单的推导可以发现，修改了\(f(a,b)\)之后，设\(b\equiv m(mod\ a),a\equiv n(mod\ b)\)，则首先受到影响的就是\(f(a,m),f(m,a),f(b,n),f(n,b)\)。
然后进一步搞下去，就会发现这是个类似于辗转相减法的东西，设\(gcd(a,b)=gcd(x,y)\)，则最终受影响的应该就是所有的\(f(x,y)\)。
令\(gcd(a,b)=d\)，而我们有\(x\times y\times f(d,d)=d^2\times f(x,y)\)（显然可以从原式推得），因此：
也就是说，对于一次修改，我们可以快速计算出\(f(d,d)\)修改后的值，而其他所有的值都可以从\(f(d,d)\)推导得出。
于是，原本很复杂的题面，立刻简单了许多。
询问
又是一波大力推式子：
显然，我们可以枚举\(gcd\)，然后就能把从上面得到的\(f(x,y)=\frac{f(d,d)\times xy}{d^2}\)这个式子代入，得到：
按照莫比乌斯反演的一般套路，有\(\sum_{p|i,p|j}\mu(p)=[gcd(i,j)=1]\)，然后先枚举\(p\)，得到：
从后式中提出\(p\)，并设\(S(n)=\sum_{i=1}^ni\)，于是就有：
这看起来似乎不能再进一步化简了，但复杂度依旧堪忧。
事实上，还有一种很神奇的差分套路。
显然有\(\lfloor\frac n{i}\rfloor-\lfloor\frac {n-1}i\rfloor=[i|n]\)。
因此，如果我们设\(F(x)=\sum_{p=1}^x\mu(p)p^2S(\lfloor\frac xp\rfloor)^2\)，并用\(F(x)\)减去\(F(x-1)\)（差分），就能得到：
单独拿出其中的\(S(\frac xp)^2-S(\frac xp-1)^2\)来分析，对于\((\sum_{i=1}^ai)^2-(\sum_{i=1}^{a-1}i)^2\)，运用平方差公式+高斯求和得到：
因此，原式就等同于：
然后我们就发现了，式子中的\(\sum_{p|n}\mu(p)\times Id(\frac xp)\)显然是一个狄利克雷卷积啊！
而我们知道\(\mu*Id=\phi\)，因此就得到了很简单的一个式子：
那么只要一遍线性筛筛出\(\phi\)，就能方便地计算出\(F(x)\)了。
最终，我们只要通过分块实现对\(f(d,d)\)的\(O(\sqrt n)\)修改、\(O(1)\)求值，并用除法分块计算答案即可。
代码
#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 4000000
#define X 1000000007
#define Qinv(x) Qpow(x,X-2)
#define Inc(x,y) ((x+=(y))>=X&&(x-=X))
#define XSum(x,y) ((x)+(y)>=X?(x)+(y)-X:(x)+(y))
#define XSub(x,y) ((x)<(y)?(x)-(y)+X:(x)-(y))
using namespace std;
int n;
class LinearSieve//线性筛
{
private:
int Pt,P[N+5],phi[N+5];
public:
int f[N+5];
I void Sieve()
{
RI i,j,t;for(phi[1]=1,i=2;i<=n;++i)//筛phi
for(!P[i]&&(phi[P[++Pt]=i]=i-1),j=1;j<=Pt&&(t=i*P[j])<=n;++j)
if(P[t]=1,i%P[j]) phi[t]=phi[i]*(P[j]-1);else {phi[t]=phi[i]*P[j];break;}
for(i=1;i<=n;++i) f[i]=(1LL*i*i%X*phi[i]+f[i-1])%X;//累加求出f
}
}S;
class BlockSolver//分块
{
private:
int sz,a[N+5],bl[N+5],V[N+5],BV[N+5];
I void BF(CI l,CI r,CI v) {for(RI i=l;i<=r;++i) Inc(V[i],v);}//散块暴力
public:
I void Init()
{
sz=sqrt(n);for(RI i=1;i<=n;++i)//预处理，赋初始值
bl[i]=(i-1)/sz+1,a[i]=1LL*i*i%X,V[i]=XSum(V[i-1],a[i]);
}
I void U(CI p,CI nv)//修改
{
RI v=XSub(nv,a[p]);a[p]=nv;//计算变化值
if(bl[p]==bl[n]) return BF(p,n,v);BF(p,bl[p]*sz,v),BF((bl[n]-1)*sz+1,n,v);
for(RI i=bl[p]+1;i^bl[n];++i) Inc(BV[i],v);//整块标记
}
I int Q(CI x) {return XSum(V[x],BV[bl[x]]);}//询问
}B;
I int gcd(CI x,CI y) {return y?gcd(y,x%y):x;}
I int Qpow(RI x,RI y) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
int main()
{
RI Qt,a,b,k,t;long long x;scanf("%d%d",&Qt,&n),S.Sieve(),B.Init();W(Qt--)
{
scanf("%d%d%lld%d",&a,&b,&x,&k),t=gcd(a,b),B.U(t,(x%X)*t%X*t%X*Qinv(1LL*a*b%X)%X);//计算f(d,d)
for(t=0,a=1;a<=k;a=b+1) b=k/(k/a),t=(1LL*(B.Q(b)-B.Q(a-1)+X)*S.f[k/a]+t)%X;printf("%d\n",t);//除法分块求答案
}return 0;
}
