P3327 [SDOI2015]约数个数和 题解

   日期:2020-05-02     浏览:118    评论:0    
核心提示:博客园同步原题链接简要题意:求∑i=1n∑j=1md(ij)\\sum_{i=1}^n \\sum_{j数据结构与算法

博客园同步

原题链接

简要题意:

∑ i = 1 n ∑ j = 1 m d ( i j ) \sum_{i=1}^n \sum_{j=1}^m d(ij) i=1nj=1md(ij)

其中, d ( x ) d(x) d(x) 表示 x x x 的因数个数。

算法一

爆搜。

时间复杂度: O ( T n m n m ) O(Tnm \sqrt{nm}) O(Tnmnm ).( T T T 的飞起)

期望得分: 0 p t 0pt 0pt.

算法二

考虑每个数作为其它数因数所产生的贡献。

时间复杂度: O ( T × min ⁡ ( n , m ) ) O(T \times \min(n,m)) O(T×min(n,m)). (海星吧)

期望得分: 0 p t 0pt 0pt.

算法三

一言不合就推式子

∑ i = 1 n ∑ j = 1 m d ( i j ) \sum_{i=1}^n \sum_{j=1}^m d(ij) i=1nj=1md(ij)

= ∑ i = 1 n ∑ j = 1 m ∑ x ∣ i ∑ y ∣ j [ gcd ⁡ ( x , y ) = = 1 ] = \sum_{i=1}^n \sum_{j=1}^m \sum_{x|i} \sum_{y|j} [\gcd(x,y)==1] =i=1nj=1mxiyj[gcd(x,y)==1]

= ∑ x = 1 n ∑ y = 1 m [ gcd ⁡ ( x , y ) = = 1 ] ∑ i = 1 ⌊ n x ⌋ ∑ j = 1 ⌊ m y ⌋ 1 = \sum_{x=1}^n \sum_{y=1}^m [\gcd(x,y)==1] \sum_{i=1}^{\lfloor \frac{n}{x} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{y} \rfloor} 1 =x=1ny=1m[gcd(x,y)==1]i=1xnj=1ym1

= ∑ x = 1 n ∑ y = 1 m [ gcd ⁡ ( x , y ) = = 1 ] ⌊ n x ⌋ × ⌊ m y ⌋ = \sum_{x=1}^n \sum_{y=1}^m [\gcd(x,y)==1] {\lfloor \frac{n}{x} \rfloor} \times {\lfloor \frac{m}{y} \rfloor} =x=1ny=1m[gcd(x,y)==1]xn×ym

= ∑ i = 1 n ∑ j = 1 m [ gcd ⁡ ( i , j ) = = 1 ] ⌊ n i ⌋ × ⌊ m j ⌋ = \sum_{i=1}^n \sum_{j=1}^m [\gcd(i,j)==1] {\lfloor \frac{n}{i} \rfloor} \times {\lfloor \frac{m}{j} \rfloor} =i=1nj=1m[gcd(i,j)==1]in×jm

显然到了这里式子推完,我们需要 莫比乌斯反演

f x = ∑ i = 1 n ∑ j = 1 m [ gcd ⁡ ( i , j ) = = x ] ⌊ n i ⌋ × ⌊ m j ⌋ f_x = \sum_{i=1}^n \sum_{j=1}^m [\gcd(i,j)==x] {\lfloor \frac{n}{i} \rfloor} \times {\lfloor \frac{m}{j} \rfloor} fx=i=1nj=1m[gcd(i,j)==x]in×jm

g x = ∑ x ∣ i f i g_x = \sum_{x|i} f_i gx=xifi

g x = ∑ i = 1 n ∑ j = 1 m [ x ∣ gcd ⁡ ( i , j ) ] ⌊ n i ⌋ × ⌊ m j ⌋ g_x = \sum_{i=1}^n \sum_{j=1}^m [x | \gcd(i,j)] {\lfloor \frac{n}{i} \rfloor} \times {\lfloor \frac{m}{j} \rfloor} gx=i=1nj=1m[xgcd(i,j)]in×jm

= ∑ i = 1 n ∑ j = 1 m ⌊ n i x ⌋ × ⌊ m j x ⌋ = \sum_{i=1}^n \sum_{j=1}^m {\lfloor \frac{n}{ix} \rfloor} \times {\lfloor \frac{m}{jx} \rfloor} =i=1nj=1mixn×jxm

答案为 f 1 f_1 f1,由 莫比乌斯反演 得:

f x = ∑ x ∣ d μ d x g d f_x = \sum_{x|d} \mu_\frac{d}{x} g_d fx=xdμxdgd

f 1 = ∑ x ∣ 1 μ d 1 g d f_1 = \sum_{x|1} \mu_\frac{d}{1} g_d f1=x1μ1dgd

= ∑ x = 1 n μ d g d = ∑ i = 1 n μ i g i = \sum_{x=1}^n \mu_d g_d = \sum_{i=1}^n \mu_i g_i =x=1nμdgd=i=1nμigi

再瞅一眼这个 g g g 的定义:

g x = ∑ i = 1 n ∑ j = 1 m ⌊ n i x ⌋ × ⌊ m j x ⌋ g_x = \sum_{i=1}^n \sum_{j=1}^m {\lfloor \frac{n}{ix} \rfloor} \times {\lfloor \frac{m}{jx} \rfloor} gx=i=1nj=1mixn×jxm

首先我们需要预处理

s x = ∑ i = 1 x ⌊ x i ⌋ s_x = \sum_{i=1}^x \lfloor \frac{x}{i} \rfloor sx=i=1xix

可得:

g x = s ⌊ n x ⌋ × s ⌊ m x ⌋ g_x = s_{\lfloor \frac{n}{x} \rfloor} \times s_{\lfloor \frac{m}{x} \rfloor} gx=sxn×sxm

显然 s s s 可以 整除分块,所以 g g g 也可以 整除分块,预处理 μ \mu μ 即可。

时间复杂度: O ( n n + T n ) O(n \sqrt{n} + T \sqrt{n}) O(nn +Tn )

实际得分: 100 p t s 100pts 100pts.

#pragma GCC optimize(2)
#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
const int N=1e5+1;

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

int mu[N],prime[N],tot=0; bool h[N];
ll s[N],ans=0; int T,n,m;

inline void Euler() {
	mu[1]=1;
	for(int i=2;i<N;i++) {
		if(!h[i]) prime[++tot]=i,mu[i]=-1;
		for(int j=1;j<=tot && i*prime[j]<N;j++) {
			h[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<N;i++) mu[i]+=mu[i-1]; //线性筛并做前缀和
	for(int i=1;i<N;i++) {
		ll k=0;
		for(int j=1;j<=i;) {
			int t=i/(i/j); k+=1ll*(t-j+1)*(i/j);
			j=t+1;
		} s[i]=k;
	} //整除分块预处理 s
}

int main() {
	Euler(); T=read(); while(T--) {
		n=read(),m=read();
		ans=0; for(int i=1;i<=min(n,m);) {
			int t=min(n/(n/i),m/(m/i));
			ans+=1ll*(mu[t]-mu[i-1])*s[n/i]*s[m/i];
			i=t+1; //整除分块计算答案
		} printf("%lld\n",ans);
	}
	return 0;
}

 
打赏
 本文转载自:网络 
所有权利归属于原作者,如文章来源标示错误或侵犯了您的权利请联系微信13520258486
更多>最近资讯中心
更多>最新资讯中心
0相关评论

推荐图文
推荐资讯中心
点击排行
最新信息
新手指南
采购商服务
供应商服务
交易安全
关注我们
手机网站:
新浪微博:
微信关注:

13520258486

周一至周五 9:00-18:00
(其他时间联系在线客服)

24小时在线客服