BZOJ 2154 Crash的数字表格

BZOJ 2154 Crash的数字表格

化简式子即可：

\[\begin{aligned} \sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)&=\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{i\times j}{gcd(i,j)} \\ &=\sum_{d=1}^{\min(n,m)}d\times\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}i\times j\times [gcd(i,j)=1] \end{aligned} \]

设 \(f(x,y)=\sum_{i=1}^x\sum_{j=1}^{m}{i\times j\times [gcd(i,j)=1]}\)

\[\begin{aligned} f(x,y)&=\sum_{i=1}^x\sum_{j=1}^yi\times j\sum_{d|gcd(i,j)}\mu(d) \\ &=\sum_d^{\min(x,y)}\mu(d)\times d^2\times\sum_{i=1}^{\lfloor\frac{x}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{y}{d}\rfloor}i\times j \\ &=\sum_{d}^{\min(n,m)}\mu(d)\times d^2\times(\lfloor\frac{x}{d}\rfloor\times(\lfloor\frac{x}{d}\rfloor+1)/2+\lfloor\frac{y}{d}\rfloor\times(\lfloor\frac{y}{d}\rfloor+1)/2) \end{aligned} \]

那么答案就是：

\[\sum_{d=1}^{\min(n,m)}\times d\times f(\lfloor\frac{x}{d}\rfloor,\lfloor\frac{y}{d}\rfloor) \]

注意到外层内层都有 \(\lfloor\frac{x}{d}\rfloor,\lfloor\frac{y}{d}\rfloor\) 这两个式子，我们可以用数论分块进行化简，同时用两个指针维护所在块的左指针，由于两个指针互不干扰，每个指针最多跳 \(\sqrt{n}\) 次，所以复杂度仍然是 \(\sqrt{n}\)。那么外层的复杂度是\(\sqrt(n)\)。内层 \(\mu(d)\times d^2\) 可以用前缀和维护然后 \(O(1)\) 查询，对于后面那一串仍然是数论分块维护，那么内层复杂度也是 \(\sqrt(n)\)。总时间复杂度为 \(O(n)\)。

数论分块套数论分块绝了

代码如下：

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int MAXN = 1e7+5;
const int MX = 1e7;
const int MOD = 20101009;
int prime[MAXN],pc,p[MAXN],n,m;
ll mu[MAXN],inv2;
ll qpw(ll x,ll b)
{
	ll now=1,tmp=x,ans=1;
	while(now<=b)
	{
		if(now&b) ans=ans*tmp%MOD;
		tmp=tmp*tmp%MOD;
		now<<=1;
	}
	return ans;
}
ll f(int x,int y)
{
	ll ans=0;
	int l1=1,l2=1,now=1,up=min(x,y);
	while(now<=up)
	{
		int r;
		r=min(x/(x/l1),y/(y/l2));
		r=min(r,up);
		ans=(ans+(mu[r]-mu[now-1]+MOD)%MOD*(x/l1+1)%MOD*(x/l1)%MOD*inv2%MOD*(y/l2+1)%MOD*(y/l2)%MOD*inv2%MOD)%MOD;
		now=r+1;
		if(x/(x/l1)<=y/(y/l2)) l1=x/(x/l1)+1;
		else l2=y/(y/l2)+1;
	}
	return ans;
}
int main()
{
	inv2=qpw(2,MOD-2);
	scanf("%d %d",&n,&m);
	int up=max(n,m);
	mu[1]=1;
	for(int i=2;i<=MX;++i)
	{
		if(!p[i]) prime[++pc]=i,mu[i]=-1;
		for(int j=1;j<=pc&&1ll*i*prime[j]<=MX;++j)
		{
			p[i*prime[j]]=1;
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				break;
			}
			mu[i*prime[j]]=-mu[i];
		}
	}
	for(int i=1;i<=MX;++i) mu[i]=mu[i]*i%MOD*i%MOD;
	for(int i=1;i<=MX;++i) mu[i]=(mu[i]+mu[i-1])%MOD;
	up=min(n,m);
	int l1=1,l2=1,r1=1,r2=1,now=1;
	ll ans=0;
	while(now<=up)
	{
		int r;
		r=min(n/(n/l1),m/(m/l2));
		r=min(up,r);
		ans=(ans+1ll*(r-now+1)*(now+r)/2*(f(n/l1,m/l2)))%MOD;
		now=r+1;
		if(n/(n/l1)<=m/(m/l2)) l1=n/(n/l1)+1;
		else l2=m/(m/l2)+1;
	}
	printf("%lld\n",ans);
	return 0;
}