[BZOJ2154]Crash的数字表格

题面戳我
题意:求

\[\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)\quad mod \quad 20101009 \]

(我不知道那个是不是质数不过有关系吗)
(SYC手玩告诉我这个数是质数,真是太厉害了)

sol

首先

\[ans=\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)=\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{ij}{\gcd(i,j)} \]

这一步没问题吧
然后我们提出\(\gcd\)(假设\(n \le m\)

\[=\sum_{d=1}^{n}\frac{f(n,m,d)}{d} \]

其中$$f(n,m,d)=\sum_{i=1}{n}\sum_{j=1}ij*[\gcd(i,j)==d]$$
(理解一下。。。)
然后可以知道

\[f(n,m,d)=d^2f(\lfloor\frac nd\rfloor,\lfloor\frac md\rfloor,1) \]

其实就是把\(d\)提出来了而已
所以

\[ans=\sum_{d=1}^{n}d*f(\lfloor\frac nd\rfloor,\lfloor\frac md\rfloor,1) \]

接下来我们考虑\(f\)怎么求,想方设法往莫比乌斯反演上靠
我们令

\[F(n,m,d)=\sum_{i=1}^{n}\sum_{j=1}^{m}ij*[d|\gcd(i,j)] \]

这样子显然就对应了一个关系是

\[F(x)=\sum_{x|d}^{n}f(d) \]

(你就当做前面的\(n\)\(m\)看不见嘛。。。)
所以反演一下就有

\[f(x)=\sum_{x|d}^{n}\mu(d)F(\frac dx) \]

这里求的是\(f(1)\)所以

\[f(1)=\sum_{d=1}^{n}\mu(d)F(d) \]

现在我们又要考虑\(F\)怎么求(。。。)
\(F\)里面的\(d\)提出来(就像前面提\(f\)一样)

\[F(n,m,d)=d^2\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}ij \]

因为\([1|\gcd(i,j)]\)恒成立
那后面这个东西不就可以直接\(O(1)\)算啦?(别告诉我这个都不会算)
还是讲一下吧

\[\sum_{i=1}^{x}\sum_{j=1}^{y}ij=\frac {x(x+1)}{2} \frac {y(y+1)}{2} \]

(其实就是两个等差数列求和再乘起来)
所以总结一下,形成了下面两个式子

\[ans=\sum_{d=1}^{n}d*f(\lfloor\frac nd\rfloor,\lfloor\frac md\rfloor,1) \]

\[f(n,m,1)=\sum_{d=1}^{n}\mu(d)d^2\frac {\lfloor\frac nd\rfloor(\lfloor\frac nd\rfloor+1)}{2} \frac {\lfloor\frac md\rfloor(\lfloor\frac md\rfloor+1)}{2} \]

直接搞的化复杂度俨然\(O(n^2)\)但是发现两式都可以数论分块所以复杂度就是\(O(\sqrt n *\sqrt n)=O(n)\)
(莫比乌斯的线性筛做到输入的\(n\)就行了别每次都做到\(10^7\)

code

#include<cstdio>
#include<algorithm>
using namespace std;
#define ll long long
const int mod = 20101009;
const int N = 10000000;

int n,m,ans,mu[N+5],pri[N+5],tot,zhi[N+5],s[N+5];
void Mobius()
{
	mu[1]=zhi[1]=1;
	for (int i=2;i<=n;i++)
	{
		if (!zhi[i]) pri[++tot]=i,mu[i]=-1;
		for (int j=1;j<=tot&&i*pri[j]<=N;j++)
		{
			zhi[i*pri[j]]=1;
			if (i%pri[j]) mu[i*pri[j]]=-mu[i];
			else {mu[i*pri[j]]=0;break;}
		}
	}
	for (int i=1;i<=n;i++)
		s[i]=(s[i-1]+1ll*mu[i]*i*i%mod+mod)%mod;
}
int Sum(int x,int y){return (1ll*x*(x+1)/2%mod)*(1ll*y*(y+1)/2%mod)%mod;}
int calc(int a,int b)
{
	int i=1,res=0;
	while (i<=a)
	{
		int j=min(a/(a/i),b/(b/i));
		res=(res+1ll*(s[j]-s[i-1]+mod)%mod*Sum(a/i,b/i)%mod)%mod;
		i=j+1;
	}
	return res;
}
int main()
{
	scanf("%d %d",&n,&m);
	if (n>m) swap(n,m);
	Mobius();
	int i=1;
	while (i<=n)
	{
		int j=min(n/(n/i),m/(m/i));
		ans=(ans+1ll*(j+i)*(j-i+1)/2%mod*calc(n/i,m/i)%mod)%mod;
		i=j+1;
	}
	printf("%d\n",ans);
	return 0;
}

然而

你以为这样就完了么?
题面戳我(权限题别问我是怎么搞到的)
就是和上面相同,然后一万组询问。
呵呵,\(O(n)\)的死掉了吧。
所以,我们进一步考虑这道题怎么优化到单组询问\(O(\sqrt n)\)

\[ans=\sum_{d=1}^{n}d*\sum_{i=1}^{n/d}\mu(i)i^2\frac {\lfloor\frac n{id}\rfloor(\lfloor\frac n{id}\rfloor+1)}{2} \frac {\lfloor\frac m{id}\rfloor(\lfloor\frac m{id}\rfloor+1)}{2} \]

答案式是长这样的没错吧(别比了,就是从上面蒯下来的,只是合并了而已)
这道题的套路一样,我们令\(T=id\),然后考虑每一项$$\frac {\lfloor\frac nT\rfloor(\lfloor\frac nT\rfloor+1)}{2} \frac {\lfloor\frac mT\rfloor(\lfloor\frac mT\rfloor+1)}{2}$$的贡献。
所以式子就化成

\[ans=\sum_{T=1}^{n}\frac {\lfloor\frac nT\rfloor(\lfloor\frac nT\rfloor+1)}{2} \frac {\lfloor\frac mT\rfloor(\lfloor\frac mT\rfloor+1)}{2}\sum_{d|T}\mu(\frac Td)\frac {T^2}{d} \]

(这一步请自行手玩。。。)
然后后面那一坨是一个积性函数(可参照狄利克雷卷积的形式),所以线性筛即可。

code

#include<cstdio>
#include<algorithm>
using namespace std;
const int mod = 100000009;
const int N = 10000000;
int gi()
{
	int x=0,w=1;char ch=getchar();
	while ((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
	if (ch=='-') w=0,ch=getchar();
	while (ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
	return w?x:-x;
}
int pri[N+5],tot,zhi[N+5],h[N];
void Mobius()
{
	h[1]=1;
	for (int i=2;i<=N;i++)
	{
		if (!zhi[i]) pri[++tot]=i,h[i]=(i-1ll*i*i%mod+mod)%mod;
		for (int j=1;j<=tot&&i*pri[j]<=N;j++)
		{
			zhi[i*pri[j]]=1;
			if (i%pri[j]) h[i*pri[j]]=1ll*h[i]*h[pri[j]]%mod;
			else {h[i*pri[j]]=1ll*h[i]*pri[j]%mod;break;}
		}
	}
	for (int i=1;i<=N;i++) h[i]=(h[i]+h[i-1])%mod;
}
int Sum(int x,int y){return (1ll*x*(x+1)/2%mod)*(1ll*y*(y+1)/2%mod)%mod;}
int main()
{
	Mobius();
	int T=gi();
	while (T--)
	{
		int n=gi(),m=gi();
		if (n>m) swap(n,m);
		int i=1,ans=0;
		while (i<=n)
		{
			int j=min(n/(n/i),m/(m/i));
			ans=(ans+1ll*Sum(n/i,m/i)*(h[j]-h[i-1]+mod)%mod)%mod;
			i=j+1;
		}
		printf("%d\n",ans);
	}
	return 0;
}
posted @ 2018-01-09 09:27  租酥雨  阅读(335)  评论(0编辑  收藏  举报