【CQOI2017】小Q的表格（数论，分块）

题意：

有一个无限大的整数表格 $f$ 满足以下两条法则：

• $f(a,b)=f(b,a)$。
• $b\times f(a,a+b)=(a+b)\times f(a,b)$。

初始时 $f(a,b)=a\times b$。有 $m$ 次修改，每次修改会改变某个位置，并将所有此次修改会影响到的所有位置按照法则重置一遍。每次修改完还会询问你前 $k$ 行前 $k$ 列里所有数的和对 $10^9+7$ 取模的结果（$k$ 不固定）。

$m\le 10^4$，修改位置行数、列数和 $k$ 均不超过 $4\times 10^6$。

题解：

找性质的方向和感觉还是不太对，推了半天看题解才反应过来那两条法则实际上在讲什么……

我们先将第二条法则移项得到 $\frac{f(a,a+b)}{a+b}=\frac{f(a,b)}{b}\to \frac{f(a,b)}{b}=\frac{f(a,b-a)}{b-a}$，同时再结合第一条法则得到 $\frac{f(a,b)}{a}=\frac{f(a-b,a)}{a-b}$。二者结合一下，得到 $\frac{f(a,b)}{ab}=\frac{f(a,b-a)}{a(b-a)}=\frac{f(a-b,a)}{(a-b)a}$，类似一个辗转相除的过程，最后得到：
$$\begin{array}{rlr}\frac{f(a,b)}{ab}& =\frac{f(d,d)}{d^2}& \\ f(d\cdot a,d\cdot b)& =ab\cdot f(d,d)& [\gcd (a,b)=1]\end{array}$$
也就是说，本质上是所有的 $f(d,d)$ 确定了所有的 $f(a,b)$，而一次修改相当于修改了一次 $f(d,d)$。

答案即为：
$$\begin{array}{rl}& \sum _{d=1}^{k}f(d,d)\sum _{a=1}^{\lfloor k/d\rfloor }\sum _{b=1}^{\lfloor k/d\rfloor }ab[(a,b)=1]\end{array}$$
考虑先求出 $g(n)={\sum }_{a=1}^n{\sum }_{b=1}^{n}ab[(a,b)=1]$，这是与 $f$ 无关的，可以提前预处理出来：
$$\begin{array}{rl}g(n)& =g(n-1)-[n=1]+2\sum _{i=1}^{n}ni[(n,i)=1]\\ & =g(n-1)-[n=1]+2n\sum _{i=1}^{n}i\sum _{j|n,j|i}\mu (j)\\ & =g(n-1)-[n=1]+2n\sum _{j|n}\mu (j)\sum _{i=1}^{\lfloor n/j\rfloor }ij\\ & =g(n-1)-[n=1]+2n\sum _{j|n}\mu (j)\cdot j\cdot S(\lfloor n/j\rfloor )\end{array}$$
注意到使 $\mu (j)\ne 0$ 的 $j$ 的取值只有 $2^{\omega (n)}$ 种（其中 $\omega (n)$ 表示 $n$ 的不同质因子个数），而由于 $2^{\omega (n)}\le d(n)$，所以 ${\sum }_{i=1}^n{2}^{\omega (i)}\le {\sum }_{i=1}^{n}d(n)=O(n\log n)$。进一步地，发现当 $n=10^6$ 时 ${\sum }_{i=1}^n{2}^{\omega (i)}$ 约为 $8\times 10^6$，所以递推时只要暴力枚举每个 $j$ 即可。

---

看完题解后发现我们还能进一步简化 $g(n)$：
$$\begin{array}{rl}g(n)& =g(n-1)-[n=1]+2n\sum _{j|n}\mu (j)\cdot j\cdot S(\lfloor n/j\rfloor )\\ & =g(n-1)-[n=1]+2n\sum _{j|n}\mu (j)\cdot j\cdot S(n/j)\\ & =g(n-1)-[n=1]+n\sum _{j|n}\mu (j)\cdot j\cdot \frac{n}{j}\cdot \left(\frac{n}{j}+1\right)\\ & =g(n-1)-[n=1]+n^2\sum _{j|n}\mu (j)\left(\frac{n}{j}+1\right)\\ & =g(n-1)+n^2\sum _{j|n}\mu (j)\frac{n}{j}\\ & =g(n-1)+n^2\phi (n)\end{array}$$
于是现在的递推时严格 $O(k)$ 的了。

---

于是我们求出了所有的 $g(n)$，现在考虑怎么求答案。

如果我们能做到快速求 $f(d,d)$ 的前缀和的话，单次询问可以整除分块做到 $O(\sqrt{k})$，足以通过。

现在的问题是处理 $f$，我们需要 $O(m)$ 次单点修改 $f$，$O(m\sqrt{k})$ 次询问 $f$ 的前缀和。

为了均衡时间复杂度，我们考虑分块，这样可以做到 $O(\sqrt{k})$ 单点修改，$O(1)$ 前缀查询。总时间复杂度 $O(k+m\sqrt{k})$。

#include<bits/stdc++.h>

#define SN 2010
#define N 4000010
#define ll long long

using namespace std;

namespace modular
{
	const int mod=1000000007;
	inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
	inline int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
	inline int mul(int x,int y){return 1ll*x*y%mod;}
	inline void Add(int &x,int y){x=x+y>=mod?x+y-mod:x+y;}
	inline void Dec(int &x,int y){x=x-y<0?x-y+mod:x-y;}
	inline void Mul(int &x,int y){x=1ll*x*y%mod;}
}using namespace modular;

inline int poww(int a,int b)
{
	int ans=1;
	while(b)
	{
		if(b&1) ans=mul(ans,a);
		a=mul(a,a);
		b>>=1;
	}
	return ans;
}

inline ll read()
{
	ll x=0;
	int f=1;
	char ch=getchar();
	while(ch<'0'||ch>'9')
	{
		if(ch=='-') f=-1;
		ch=getchar();
	}
	while(ch>='0'&&ch<='9')
	{
		x=(x<<1)+(x<<3)+(ch^'0');
		ch=getchar();
	}
	return x*f;
}

int m,n,f[N];
int len,B,bl[N],sumf[N],tag[SN];
int cnt,prime[N],phi[N],g[N];
bool notprime[N];

void init()
{
	phi[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(!notprime[i])
		{
			prime[++cnt]=i;
			phi[i]=i-1;
		}
		for(int j=1,v;j<=cnt&&(v=i*prime[j])<=n;j++)
		{
			notprime[v]=1;
			if(!(i%prime[j]))
			{
				phi[v]=phi[i]*prime[j];
				break;
			}
			phi[v]=phi[i]*phi[prime[j]];
		}
	}
	for(int i=1;i<=n;i++)
		g[i]=add(g[i-1],mul(mul(i,i),phi[i]));
	len=sqrt(n);
	for(int i=1;i<=n;i++)
	{
		bl[i]=(i-1)/len+1;
		sumf[i]=add(sumf[i-1],f[i]=mul(i,i));
	}
	B=bl[n];
}

int gcd(int a,int b)
{
	return b?gcd(b,a%b):a;
}

void update(int x,int y)
{
	int b=bl[x];
	for(int i=x,r=min(n,b*len);i<=r;i++) Add(sumf[i],y);
	for(int i=b+1;i<=B;i++) Add(tag[i],y);
}

int query(int l,int r)
{
	return dec(add(sumf[r],tag[bl[r]]),add(sumf[l-1],tag[bl[l-1]]));
}

int main()
{
	m=read(),n=read();
	init();
	while(m--)
	{
		int a=read(),b=read(),k;
		ll y=read();k=read();
		int d=gcd(a,b);
		int fd=(y/(a/d)/(b/d))%mod;
		update(d,dec(fd,f[d]));
		f[d]=fd;
		int ans=0;
		for(int l=1,r;l<=k;l=r+1)
		{
			r=k/(k/l);
			Add(ans,mul(query(l,r),g[k/l]));
		}
		printf("%d\n",ans);
	}
	return 0;
}