BZOJ 4815: [Cqoi2017]小Q的表格

推到一半被一个函数卡住真是太屑了的说（菜是原罪）

首先对于(b imes f(a,a+b)=(a+b) imes f(a,b))我们把它化成一个整齐一点的式子：(frac{f(a,a+b)}{a(a+b)}=frac{f(a,b)}{ab})

这样是不是看出点什么了，然后考虑(f)里的变化，是不是相当于一个数减去另一个数

如果我们减到不能减，然后再换过来继续减，直到两数相等，容易发现此时我们相当于在进行辗转相减的过程，实际上容易发现会递归到(gcd(a,b))，因此我们有：

[f(a,b)=frac{ab imes f(d,d)}{d^2},d=gcd(a,b) ]

那么就意为着我们只需要考虑(f(d,d))对棋盘的影响即可，考虑对于前(n)行前(n)列的棋盘求和：

[ans=sum_{i=1}^nsum_{j=1}^n frac{ij imes f(gcd(i,j),gcd(i,j))}{gcd(i,j)^2}\=sum_{d=1}^n frac{f(d,d)}{d^2} sum_{i=1}^nsum_{j=1}^n ij[gcd(i,j)=d]\=sum_{d=1}^n f(d,d) sum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{n}{d} floor} ij[gcd(i,j)=1]\=sum_{d=1}^n f(d,d) sum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{n}{d} floor} ijsum_{k|i,k|j}mu(k)\=sum_{d=1}^n f(d,d) sum_{k=1}^{lfloorfrac{n}{d} floor}mu(k) imes k^2 imes (sum_{i=1}^{lfloorfrac{n}{kd} floor} i) imes (sum_{j=1}^{lfloorfrac{n}{kd} floor} j)\=sum_{d=1}^n f(d,d) sum_{k=1}^{lfloorfrac{n}{d} floor}mu(k) imes k^2 imes S(lfloorfrac{n}{kd} floor)^2\ ]

其中(S(n)=frac{n(n+1)}{2})，前面的部分可以求出前缀和然后除法分块，那么问题是后面的这部分怎么做（到这里就不会了）

看了题解发现方法确实是妙，我们令(G(n)=sum_{i=1}^n mu(i) imes i^2 imes S(lfloorfrac{n}{i} floor)^2)，发现主要是里面(S(lfloorfrac{n}{i} floor))有个下取整很难搞

考虑运用关于下取整的技巧——差分：

[lfloorfrac{n}{i} floor-lfloorfrac{n-1}{i} floor=[i|n] ]

然后就有：

[G(n)-G(n-1)=sum_{i|n} mu(i) imes i^2 imes [S(lfloorfrac{n}{i} floor)^2-S(lfloorfrac{n}{i} floor-1)^2]\=sum_{i|n}mu(i) imes i^2 imes(frac{n}{i})^3\=n^2sum_{i|n}mu(i) imesfrac{n}{i}\=n^2 imes varphi(n) ]

第二步的转化自已化开一下就出来了，最后一步的原因是(varphiast I=idLeftrightarrow varphiast Iast mu=idast muLeftrightarrow varphiast epsilon=idast muLeftrightarrowvarphi=idast mu)，就进行了这一步化(mu)为(varphi)

于是(G(n))就可以(O(n))预处理了

考虑我们还需要修改，显然每次修改都只需要修改一个对应的(f(d,d))即可，而我们同时也要维护(f(d,d))的前缀和

如果直接用树状数组维护单次询问的复杂度达到(O(sqrt nlog n))，这样可能会被卡掉

但是我们套路地发现这个修改的复杂度可以大一些，因为前面不需要乘(sqrt n)的询问次数，因此我们容易想到用分块(O(sqrt n))修改(O(1))询问，总复杂度就是(O(msqrt n))的了

PS：这题不知道为什么BZOJ上一直RE，Luogu和LOJ上都没问题，下载了数据也是对的233

#include<cstdio>
#include<cmath>
#define RI register int
#define CI const int&
using namespace std;
const int N=4e6,mod=1e9+7;
int n,m,a,b,k,prime[N+5],cnt,phi[N+5],g[N+5]; bool vis[N+5]; long long x;
inline int gcd(CI x,CI y)
{
	return y?gcd(y,x%y):x;
}
inline int sum(CI x,CI y)
{
	int t=x+y; return t>=mod?t-mod:t;
}
inline int sub(CI x,CI y)
{
	int t=x-y; return t<0?t+mod:t;
}
#define P(x) prime[x]
inline void init(CI n)
{
	RI i,j; for (vis[1]=phi[1]=1,i=2;i<=n;++i)
	{
		if (!vis[i]) prime[++cnt]=i,phi[i]=i-1;
		int tp; for (j=1;j<=cnt&&(tp=i*P(j))<=n;++j)
		{
			vis[tp]=1; if (i%P(j)) phi[tp]=phi[i]*(P(j)-1);
			else { phi[tp]=phi[i]*P(j); break; }
		}
	}
	for (i=1;i<=n;++i) g[i]=sum(g[i-1],1LL*i*i%mod*phi[i]%mod);
}
#undef P
class Value_Blocks
{
	private:
		static const int S=2005;
		int a[N+5],pfx[S][S],tag[S],bel[N+5],L[S],R[S],size;
	public:
		inline void init(CI n)
		{
			RI i,j; for (size=(int)sqrt(n),i=1;i<=n;++i) bel[i]=(i-1)/size+1;
			for (i=1;i<=bel[n];++i) L[i]=(i-1)*size+1,R[i]=i*size; R[n]=n;
			for (i=1;i<=bel[n];++i) for (j=L[i];j<=R[i];++j)
			tag[i]=sum(tag[i],a[j]=1LL*j*j%mod),pfx[i][j-L[i]+1]=sum(pfx[i][j-L[i]],a[j]);
			for (i=1;i<=bel[n];++i) tag[i]=sum(tag[i],tag[i-1]);
		}
		inline void modify(CI x,CI y)
		{
			RI i; int mv=sub(y,a[x]); a[x]=y;
			for (i=x;i<=R[bel[x]];++i) pfx[bel[x]][i-L[bel[x]]+1]=sum(pfx[bel[x]][i-L[bel[x]]+1],mv);
			for (i=bel[x];i<=bel[n];++i) tag[i]=sum(tag[i],mv);
		}
		inline int query(CI x)
		{
			//int ret=0; for (RI i=1;i<=x;++i) ret=sum(ret,a[i]); return ret;
			return x?sum(tag[bel[x]-1],pfx[bel[x]][x-L[bel[x]]+1]):0;
		}
}B;
int main()
{
	for (scanf("%d%d",&m,&n),init(n),B.init(n);m;--m)
	{
		scanf("%d%d%lld%d",&a,&b,&x,&k); int d=gcd(a,b),ans=0;
		B.modify(d,(x/(1LL*(a/d)*(b/d)))%mod); for (RI l=1,r;l<=k;l=r+1)
		r=k/(k/l),ans=sum(ans,1LL*sub(B.query(r),B.query(l-1))*g[k/l]%mod);
		printf("%d
",ans);
	}
	return 0;
}