「CQOI2017」小Q的表格 题解

八万年没写过莫反了……随便写点。

先简单的列出性质：

1. \(f(a,b)=f(b,a)\)；
2. \(b f(a,a+b) = (a+b)f(a,b)\)。

先考虑这个 \(f\) 两个恒等式的形式。将 \(f\) 视作一个状态，那么有转移 \((a,a+b) \to (a,b)\)，且 \((a,b) \to (b,a)\)。

不难联想到辗转相减求最大公约数的过程，最终 \((a,b)\) 一定可以转移到 \((\gcd(a,b),\gcd(a,b))\) 上。

那么，表格上的所有数都可以用对应的一个对角线上的数乘上一个系数表示。修改一个数的话，我们只用管它的对角线元素就行了。

注意到 \(b f(a,a+b) = (a+b)f(a,b)\)，进行简单处理后可以得到 \(\dfrac{f(a,a+b)}{a(a+b)} = \dfrac{f(a,b)}{ab}\)。于是思考辗转相减的过程，易得 \(\dfrac{f(a,b)}{ab} = \dfrac{f(\gcd(a,b),\gcd(a,b))}{\gcd(a,b)^2}\)。

下面的 \(i/j\) 均表示 \(\left\lfloor \dfrac{i}{j} \right\rfloor\)。

然后求和的工作话，我们先记 \(g(i) = f(i,i)\)，直接推式子：

\[\begin{aligned} \sum_{i=1}^k \sum_{j=1}^k f(i,j) &= \sum_{d=1}^n g(d) \sum_{i=1}^n \dfrac{i}{d} \sum_{j=1}^n \dfrac{j}{d} [\gcd(i,j)=d] \\ &= \sum_{d=1}^n g(d) \sum_{i=1}^{n/d} i \sum_{j=1}^{n/d} j[\gcd(i,j)=1] \end{aligned} \]

注意到后面的 \(\displaystyle \sum_{j=1}^{\lfloor \frac{n}{d} \rfloor} j[\gcd(i,j)=1]\) 看起来就非常可以预处理。研究一下这个东西。

首先这个上界不优美。注意到在原式中，\(i=j \neq 1\) 的时候不会造成贡献，并且 \(i \neq j\) 的贡献会算两次。我们强制 \(i>j\)，然后乘 \(2\) 即可。（下面用 \(\frac{i}{D}\) 的原因是其一定是一个整数）

\[\begin{aligned} 2\sum_{j=1}^i j[\gcd(i,j)=1] &= 2\sum_{j=1}^i j \sum_{D|\gcd(i,j)} \mu(D) \\ &= 2\sum_{D|i} \mu(D) \sum_{j=1}^{\frac{i}{D}} jD \\ &= 2\sum_{D|i} \mu(D) D \sum_{j=1}^{\frac{i}{D}} j \\ &= 2\sum_{D|i} \mu(D) D \dfrac{\left(\frac{i}{D}+1\right)\times \frac{i}{D}}{2} \\ &= \sum_{D|i} \mu(D) D \left(\frac{i}{D}+1\right)\times \frac{i}{D} \\ &= i\sum_{D|i} \mu(D) \left(\frac{i}{D}+1\right) \end{aligned} \]

推到这里基本可以结束了。注意到和式是 \(\mu * (id+1)\) 的形式，那么这个函数就相当于 \(\varphi + \epsilon\)。但是注意到，在 \(i=1\) 的时候，我们将 \(i=j=1\) 的情况算了两次……处理方式是直接把 \(\epsilon\) 直接阉割了，这样就不会重了。

那么继续推之前的式子：

\[\begin{aligned} \sum_{i=1}^k \sum_{j=1}^k f(i,j) &= \sum_{d=1}^n g(d) \sum_{i=1}^n \dfrac{i}{d} \sum_{j=1}^n \dfrac{j}{d} [\gcd(i,j)=d] \\ &= \sum_{d=1}^n g(d) \sum_{i=1}^{n/d} i^2 \varphi(i) \end{aligned} \]

后面的和式是个前缀和形式。然后预处理 \(i^2 \varphi(i)\) 的前缀和非常简单，外层直接整除分块做就可以了。

但是还有个修改。这是个方程形式，采用之前的那个关于 \(f(a,b),g(\gcd(a,b))\) 的等式，采用树状数组求出新的 \(g(\gcd(a,b))\)，直接改就行了。求 \(g\) 的区间和也可以用树状数组。

但是有点慢……还凑合，就这样吧。时间复杂度 \(O(n \sqrt n \log n)\)。不懂这个题的神性在哪里。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
char buf[1<<18],*p1=buf,*p2=buf;
#define getchar() (p1==p2 && (p2=(p1=buf)+fread(buf,1,1<<18,stdin),p1==p2)?EOF:*p1++)
LL read()
{
	LL x=0;
	char c=getchar();
	while(c<'0' || c>'9')	c=getchar();
	while(c>='0' && c<='9')	x=(x<<1)+(x<<3)+(c^'0'),c=getchar();
	return x;
}
void write(LL x)
{
	if(x>9)	write(x/10);
	putchar(x%10+'0');
}
const LL MOD=1e9+7;
inline LL Add(LL x,LL y){return x+y>=MOD?x+y-MOD:x+y;}
inline LL Sub(LL x,LL y){return x<y?x-y+MOD:x-y;}
inline LL Mul(LL x,LL y){return 1ll*x*y%MOD;}
LL QuickPow(LL x,LL p)
{
	LL ans=1,base=x;
	while(p)
	{
		if(p&1)	ans=Mul(ans,base);
		base=Mul(base,base);
		p>>=1;
	}
	return ans;
}
LL N,m;
inline LL lowbit(LL x){return x&(-x);}
struct BinaryIndexedTree{
	LL tr[4000005];
	void modify(LL x,LL val){for(LL i=x;i<=N;i+=lowbit(i))	tr[i]=Add(tr[i],val);}
	LL query(LL x){LL ans=0;for(LL i=x;i;i^=lowbit(i))	ans=Add(ans,tr[i]);return ans;}
	LL query(LL l,LL r){return Sub(query(r),query(l-1));}
}bit;
LL cnt,prime[4000005];
bool vis[4000005];
LL phi[4000005],f[4000005];
void shai(LL up)
{
	vis[phi[1]=1]=true;
	for(LL i=2;i<=up;++i)
	{
		if(!vis[i])	phi[i]=i-1,prime[++cnt]=i;
		for(LL j=1;j<=cnt && i*prime[j]<=up;++j)
		{
			vis[i*prime[j]]=true;
			if(i%prime[j]==0)
			{
				phi[i*prime[j]]=Mul(phi[i],prime[j]);
				break;
			}
			phi[i*prime[j]]=Mul(phi[i],prime[j]-1);
		}
	}
	for(LL i=1;i<=up;++i)	f[i]=Mul(Mul(i,i),phi[i]),f[i]=Add(f[i],f[i-1]);
}
LL gcd(LL a,LL b){return !b?a:gcd(b,a%b);}
LL g[4000005];
int main(){
	shai(4000000);
	m=read(),N=read();
	for(LL i=1;i<=N;++i)	bit.modify(i,g[i]=Mul(i,i));
	while(m-->0)
	{
		LL a=read(),b=read(),x=read(),n=read();
		LL d=gcd(a,b);
		bit.modify(d,MOD-g[d]);
		long long tmp=x;
		tmp/=a/d;
		tmp/=b/d;
		tmp%=MOD;
		g[d]=tmp;
		bit.modify(d,g[d]);
		LL ans=0;
		for(LL l=1,r;l<=n;l=r+1)	r=n/(n/l),ans=Add(ans,Mul(f[n/l],bit.query(l,r)));
		write(ans),puts("");
	}
	return 0;
}