「bzoj2671」Calc - 数论+莫比乌斯反演 | Bill Yang's Blog

「bzoj2671」Calc - 数论+莫比乌斯反演

题目大意

    给出$N$,统计满足下面条件的数对$(a,b)$的个数:

  1. $1\le a\lt b\le N$
  2. $a+b$ 整除 $a\cdot b$

题目分析

做这道题时跟着CQzhangyu学,然后发现他给出了假的式子和假的反演,代码却是对的。

前置技能

  • 使用$(a,b)$表示$\gcd(a,b)$
  • 若$(a,b)=1$,则$(a+b,ab)=1$
    证明(来自Hineven):
    反证法。
    设$(a+b,ab)=k\gt1$,因为$(a,b)$互质,故$k\mid a$或$k\mid b$。
    设$k\mid a$,则:

正式分析

设$d=\gcd(a,b)$
则:

问题就变成了问有多少对$a,b$满足$(a,b)=1$且$b\times d\le n$满足$a+b\mid d$,又因为$d\gt b$所以b最多只能是$m=\sqrt n$。

最后的$\lfloor \frac{n}{(a+b)b}\rfloor$是因为,$b\times d\le n$,故$d\le\lfloor\frac nb\rfloor$,又因为$a+b\mid d$,则合法的数对有$\left\lfloor\frac{\lfloor\frac nb\rfloor}{a+b}\right\rfloor=\left\lfloor\frac{n}{(a+b)b}\right\rfloor$个。
反演一波:

不难发现,最后这一部分:

可以枚举$a+b$整除分块。

时间复杂度由简单积分可知为$O(n^{\frac34}-n^{\frac12})$。


代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
inline int Get_Int() {
int num=0,bj=1;
char x=getchar();
while(!isdigit(x)) {if(x=='-')bj=-1;x=getchar();}
while(isdigit(x)) {num=num*10+x-'0';x=getchar();}
return num*bj;
}
const int maxn=1000005;
int vst[maxn],p[maxn],u[maxn],cnt=0;
void Mobius_Table(int n) {
u[1]=1;
for(int i=2; i<=n; i++) {
if(!vst[i]) {p[++cnt]=i;u[i]=-1;}
for(int j=1; j<=cnt&&i*p[j]<=n; j++) {
vst[i*p[j]]=1;
if(i%p[j]==0) {u[i*p[j]]=0;break;}
u[i*p[j]]=-u[i];
}
}
}
LL Cal(int n,int m) {
LL ans=0;
for(int b=1; b<=m; b++) {
int up=n/b;
for(int i=b+1,last; i<(b<<1)&&i<=up; i=last+1) {
last=min((b<<1)-1,up/(up/i));
ans+=1ll*(last-i+1)*(up/i);
}
}
return ans;
}
int n,m;
int main() {
n=Get_Int();
m=sqrt(n);
Mobius_Table(m);
LL ans=0;
for(int d=1; d<=m; d++)ans+=u[d]*Cal(n/d/d,m/d);
printf("%lld\n",ans);
return 0;
}
本篇文章有用吗?有用还不快点击!
0%