「NOI2018模拟2」「JZOJ5134」「Codeforces 235E」除数 divisor - 莫比乌斯反演+数论/乱搞 | Bill Yang's Blog

「NOI2018模拟2」「JZOJ5134」「Codeforces 235E」除数 divisor - 莫比乌斯反演+数论/乱搞

题目大意

    求:


题目分析

在做此题前建议先完成这两道题

标解看这两篇好了:12

接下来讲一讲我的乱搞方法。

首先这个结论相信大家都知道:

不知道也没关系,我在这儿证过。

显然这个结论可以推广到三维:

故原题的式子变为:

把后面那坨$\sum_{z=1}^C[(y,z)=1,(x,z)=1]\lfloor\frac Cz\rfloor$拿出来莫比乌斯反演。

令$f(x)=\sum_{i=1}^x\lfloor\frac xi\rfloor$,这个$f(x)$就是$1\sim x$的约数个数和,可以线性筛出来。
那么原式就变为了:

然而这个东西是没什么性质的,暴力求解是$O(n)$的,显然不能接受。
这个时候需要玄学,考虑化简前的式子:

若某两个数$x,y$满足条件,将$x,y$质因数分解得到形似$p_1^{a_1}\cdot p_2^{a_2}\cdots$的式子,此时将所有的$a_i$都变为$1$,得到的$x’,y’$仍然满足条件$[(y’,z)=1,(x’,z)=1]$,而后面的计算式$\lfloor\frac Cz\rfloor$与$x,y$无关,因此对于化简后得到$x’,y’$的$x,y$,答案是一样的,故可以记忆化一下,转化为求解:

然后就可以暴力计算了,暴力枚举约数计算(这里的暴力不是枚举约数,而是每次去掉最小约数的DFS)。
经过打表发现,DFS只进入了一千万次,只执行了六千万次,也就是说每次的递归只有$6$层,$O(n^2)$的常数级别。
故可以通过了,时间瓶颈在于求$\gcd$,可以记忆化一下。
时间复杂度近似$O(n^2)$。


代码

代码中$mds[],d[]$用于线性筛约数个数,$mind[]$是最小质因数,$simp[x]$即为$x’$。

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
#include<bits/stdc++.h>
using namespace std;
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=5005,mod=1<<30;
int cnt=0,p[maxn],d[maxn],u[maxn],mds[maxn],mind[maxn],simp[maxn],f[maxn];
bool vst[maxn];
void Sieve(int n) {
d[1]=u[1]=mind[1]=simp[1]=1;
for(int i=2; i<=n; i++) {
if(!vst[i]) {
p[++cnt]=i;
mds[i]=1,d[i]=2;
mind[i]=simp[i]=i;
u[i]=-1;
}
for(int j=1; j<=cnt&&i*p[j]<=n; j++) {
vst[i*p[j]]=1;
mind[i*p[j]]=p[j];
if(i%p[j]==0) {
mds[i*p[j]]=mds[i]+1;
d[i*p[j]]=d[i]/(mds[i]+1)*(mds[i]+2);
simp[i*p[j]]=simp[i];
u[i*p[j]]=0;
break;
}
mds[i*p[j]]=1;
d[i*p[j]]=d[i]*d[p[j]];
simp[i*p[j]]=simp[i]*p[j];
u[i*p[j]]=-u[i];
}
}
for(int i=1; i<=n; i++)f[i]=f[i-1]+d[i];
}
int Gcd[maxn][maxn];
int gcd(int a,int b) {
if(Gcd[a][b])return Gcd[a][b];
return Gcd[a][b]=!b?a:gcd(b,a%b);
}
int A,B,C,M[maxn*maxn];
int Dfs(int x,int y,int d) {
if(d>C)return 0;
if(x==1&&y==1)return u[d]*f[C/d];
else if(x>1)return Dfs(x/mind[x],y,d)+Dfs(x/mind[x],y,d*mind[x]);
else return Dfs(x,y/mind[y],d)+Dfs(x,y/mind[y],d*mind[y]);
}
int main() {
Sieve(5000);
A=Get_Int();
B=Get_Int();
C=Get_Int();
int ans=0;
for(int x=1; x<=A; x++)
for(int y=1; y<=B; y++) {
if(gcd(simp[x],simp[y])!=1)continue;
int T=simp[x]*simp[y];
if(!M[T])M[T]=Dfs(x,y,1);
ans+=(A/x)*(B/y)*M[T];
}
printf("%d\n",ans&(mod-1));
return 0;
}

本篇文章有用吗?有用还不快点击!
0%