LGOJ3327 【SDOI2015】约数个数和

又是一道卡常好题

坑掉我的 (define space int space long space long) 感觉出题人并没有获得什么快乐……

Description

link

题意概述:

(d(x))(x) 的约数个数,求:

[sum_{i=1}^n sum_{j=1}^n d(ij) ]

多组询问,(1leq n,m,T leq 5 imes10^4)

Solution

首先给一个引理:

[d(ij)=sum_{x|i}sum_{y|j} [gcd(x,y)=1] ]

证明:参考(@Siyuan)的博客,好像可以用映射法?

[Begin ]

我们把式子转化一下:

[sum_{i=1}^n sum_{j=1}^n sum_{x|i}sum_{y|j} [gcd(x,y)=1] ]

(四个(sum)好恐怖?)

然后改变求和的次序,首先枚举因数 (x)(y)

(这一步转化可以考虑用(1)(a)(b)的倍数个数来理解)

[sum_{x=1}^n sum_{y=1}^n lfloor frac{n}{x} floor lfloor frac{m}{y} floor[gcd(x,y)=1] ]

换一下,我们把这个 (x)(y) 换成 (i)(j)

[sum_{i=1}^n sum_{j=1}^n lfloor frac{n}{i} floor lfloor frac{m}{j} floor[gcd(i,j)=1] ]

下一步开始反演,根据套路式,设:

[f(x)=sum_{i=1}^n sum_{j=1}^n lfloor frac{n}{i} floor lfloor frac{m}{j} floor[gcd(i,j)=1] ]

[g(x)=sum_{x|d} f(d) ]

所以得到关于 (g(x)) 的表达式:

[g(x)=sum_{i=1}^n sum_{j=1}^n lfloor frac{n}{i} floor lfloor frac{m}{j} floor[x|gcd(i,j)] ]

(x) 提出来,同时去掉 (gcd)

[g(x)=sum_{i=1}^{frac{n}{x}}sum_{j=1}^{frac{m}{x}} lfloor frac{n}{ix} floor lfloor frac{m}{jx} floor ]

再根据(f(x))的定义,得到答案为(f(1))

又有:

[f(x)=sum_{n|d} mu(frac{d}{n}) g(d) ]

所以有

[Ans=sum_{i=1}^n mu(i)g(i) ]

我们预处理(s(x)=sum^{x}_ {i=1} lfloor frac{x}{i} floor)然后就可以快速求解本题

[Finished ]

Code

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define reg register
namespace yspm{
	inline int read()
	{
		int res=0,f=1; char k;
		while(!isdigit(k=getchar())) if(k=='-') f=-1;
		while(isdigit(k)) res=res*10+k-'0',k=getchar();
		return res*f;
	}
	const int N=5e4+10;
	int tot,mu[N],pri[N];ll s[N];
	bool fl[N];
	inline void prework()
	{
		mu[1]=1;
		for(reg int i=2;i<N;++i)
		{
			if(!fl[i]) mu[i]=-1,pri[++tot]=i;
			for(reg int j=1;j<=tot&&i*pri[j]<N;++j)
			{
				fl[i*pri[j]]=1;
				if(i%pri[j]==0){mu[i*pri[j]]=0; break;}
				else mu[i*pri[j]]=-mu[i];
			}
		} for(reg int i=1;i<N;++i) mu[i]+=mu[i-1];
		for(reg int x=1;x<N;++x)
		{
			ll res=0; 
			for(reg int l=1,r;l<=x;l=r+1){r=x/(x/l); res+=1ll*(r-l+1)*(x/l);}
			s[x]=res;
		}
		return ;
	}
	inline void work()
	{
		int n=read(),m=read(); if(n>m) swap(n,m);
		ll ans=0; for(reg int l=1,r;l<=n;l=r+1){r=min(n/(n/l),m/(m/l)); ans+=1ll*(mu[r]-mu[l-1])*s[n/l]*s[m/l];}		
		return printf("%lld
",ans),void();
	}
	signed main()
	{
		prework(); int T=read(); while(T--) work();
		return 0;
	}
}
signed main(){return yspm::main();} 
原文地址:https://www.cnblogs.com/yspm/p/12273318.html