「SDOI2015」约数个数和 / 「bzoj4176」Lucas的数论 - 莫比乌斯反演+杜教筛 | Bill Yang's Blog

「SDOI2015」约数个数和 / 「bzoj4176」Lucas的数论 - 莫比乌斯反演+杜教筛

题目大意(「SDOI2015」约数个数和)

    设$d(x)$为$x$的约数个数,给定$n,m$,求:

    $n,m\le50000$,多组数据$t\le50000$

题目大意(「bzoj4176」Lucas的数论)

    设$d(x)$为$x$的约数个数,给定$n$,求:

    $n\le10^9$


题目分析

两道题中后者是前者的加强版,不妨直接考虑前者的式子(后者将$m$改成$n$):

结论

证明:
设$xy=\prod\limits_{i=1}^kp_i^{a_i}$,且$x=\prod\limits_{i=1}^kp_i^{b_i}$,$y=\prod\limits_{i=1}^kp_i^{a_i-b_i}$。
设$i=\prod\limits_{i=1}^kp_i^{c_i}$,$j=\prod\limits_{i=1}^kp_i^{d_i}$。
若$(i,j)=1$,那么$c_i=0$或$d_i=0$,那么分别$d_i$有$a_i-b_i+1$种取值,$c_i$有$b_i+1$种取值,共有$a_i+1$种取值。
因此满足条件的$i,j$对数有$\prod\limits_{i=1}^k(x_i+1)$,与约数个数定理的形式相同,证毕。

令$f(x)=\sum\limits_{i=1}^x\lfloor\frac xi\rfloor$,则:

显然,若能够求出$f(x)$,剩余部分可以使用分块跳跃$O(\sqrt n)$解决(线性筛$\mu$)。

那么对于第一题,只需要暴力$O(n\sqrt n)$预处理$f(x)$,就可以在$O(n\sqrt n+T\sqrt n)$的时间内解决。

这显然不能满足第二题的需求。
考虑,有效的$f(x)$只有$O(\sqrt n)$个,因此若仅求出这$O(\sqrt n)$个$f(x)$,即可降低时间,在$O(n+T\sqrt n)$的时间内解决。

这样第一题就可以很快解决了。
但仍然不能满足第二题$n\le10^9$的限制。
主要的问题有两个,$\mu$不能线性筛,$f(x)$无法快速求出。

不妨考虑杜教筛。
预处理$4000000$以内的$\mu(x)$,剩下的调用莫比乌斯函数杜教筛解决。
再考虑刚才对$f(x)$的优化,其时间复杂度并不是$O(\sqrt n\times\sqrt n)=O(n)$,因为分块的范围在变化。
通过积分可知总时间复杂度为$O(n^{\frac34})$,若不会积分也可以暴力用程序验证时间复杂度。

1
2
3
4
5
6
7
8
9
10
11
12
int main() {
LL sum=0,n=1000000000;
LL next=1;
for(int i=1; i<=n; i=next+1) {
next=n/(n/i);
if(i>4000000)sum+=pow(i,2.0/3.0);
if(next>4000000)sum+=pow(next,2.0/3.0);
sum+=sqrt(n/i);
}
printf("%lld\n",sum);
return 0;
}

于是本题解决了。


代码

「SDOI2015」约数个数和 $O(n\sqrt n+T\sqrt n)$

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
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
#include<algorithm>
#include<iostream>
#include<iomanip>
#include<cstring>
#include<cstdlib>
#include<vector>
#include<cstdio>
#include<cmath>
#include<queue>
using namespace std;

typedef long long LL;

inline const LL Get_Int() {
LL num=0,bj=1;
char x=getchar();
while(x<'0'||x>'9') {
if(x=='-')bj=-1;
x=getchar();
}
while(x>='0'&&x<='9') {
num=num*10+x-'0';
x=getchar();
}
return num*bj;
}

const int maxn=50005;

LL t,n,m,cnt=0,vst[maxn],Prime[maxn],u[maxn],f[maxn];

void Mobius_Table(int n) {
u[1]=1;
for(int i=2; i<=n; i++) {
if(!vst[i]) {
Prime[++cnt]=i;
u[i]=-1;
}
for(int j=1; j<=cnt&&Prime[j]*i<=n; j++) {
vst[Prime[j]*i]=1;
if(i%Prime[j]==0) {
u[i*Prime[j]]=0;
break;
}
u[i*Prime[j]]=-u[i];
}
}
}

LL Cal(LL n,LL m) {
LL next=1,ans=0;
for(int i=1; i<=min(n,m); i=next+1) {
next=min(n/(n/i),m/(m/i));
ans+=(u[next]-u[i-1])*f[n/i]*f[m/i];
}
return ans;
}

int main() {
Mobius_Table(50000);
for(int i=1; i<=50000; i++)u[i]+=u[i-1];
for(int i=1; i<=50000; i++) {
LL next=1;
for(int j=1; j<=i; j=next+1) {
next=i/(i/j);
f[i]+=(i/j)*(next-j+1);
}
}
t=Get_Int();
while(t--) {
n=Get_Int();
m=Get_Int();
printf("%lld\n",Cal(n,m));
}
return 0;
}

bzoj4176

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
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
#include<algorithm>
#include<iostream>
#include<iomanip>
#include<cstring>
#include<cstdlib>
#include<climits>
#include<vector>
#include<cstdio>
#include<cmath>
#include<queue>
using namespace std;

typedef long long LL;
#define pii pair<LL,LL>
#define mp make_pair

inline const LL Get_Int() {
LL num=0,bj=1;
char x=getchar();
while(x<'0'||x>'9') {
if(x=='-')bj=-1;
x=getchar();
}
while(x>='0'&&x<='9') {
num=num*10+x-'0';
x=getchar();
}
return num*bj;
}

const int limit=4000000;
const LL mod=1000000007,hashmod=999599;

LL n;
bool vst[limit+5];
int cnt=0,Prime[limit+5],Mobius[limit+5];
vector<pii>Hash[hashmod+5];

void Mobius_Table(int n) {
Mobius[1]=1;
for(int i=2; i<=n; i++) {
if(!vst[i]) {
Prime[++cnt]=i;
Mobius[i]=-1;
}
for(int j=1; j<=cnt&&i*Prime[j]<=n; j++) {
vst[i*Prime[j]]=1;
if(i%Prime[j]==0) {
Mobius[i*Prime[j]]=0;
break;
}
Mobius[i*Prime[j]]=-Mobius[i];
}
}
for(int i=2; i<=n; i++)Mobius[i]+=Mobius[i-1];
}

void Add(int x,LL p,LL v) {
Hash[x].push_back(mp(p,v));
}

LL Cal_u(LL n) {
if(n<=limit)return Mobius[n];
for(pii x:Hash[n%hashmod])if(x.first==n)return x.second;
LL Next,ans=0;
for(LL i=2; i<=n; i=Next+1) {
Next=n/(n/i);
ans+=(Next-i+1)*Cal_u(n/i);
}
ans=1-ans;
Add(n%hashmod,n,ans);
return ans;
}

LL Cal_f(LL x) {
LL next=1,ans=0;
for(LL i=1; i<=x; i=next+1) {
next=x/(x/i);
ans=(ans+(x/i)*(next-i+1)%mod)%mod;
}
return ans;
}

LL Cal(LL n) {
LL next=1,ans=0;
for(LL i=1; i<=n; i=next+1) {
next=n/(n/i);
LL f=Cal_f(n/i);
ans=(ans+(Cal_u(next)-Cal_u(i-1)+mod)%mod*f%mod*f%mod)%mod;
}
return ans;
}

int main() {
Mobius_Table(limit);
n=Get_Int();
printf("%lld\n",Cal(n));
return 0;
}

姥爷们赏瓶冰阔落吧~
0%