bzoj4652 [Noi2016]循环之美

Description

牛牛是一个热爱算法设计的高中生。在他设计的算法中，常常会使用带小数的数进行计算。牛牛认为，如果在k进制下，一个数的小数部分是纯循环的，那么它就是美的。现在，牛牛想知道：对于已知的十进制数n和m，在k进制下，有多少个数值上互不相等的纯循环小数，可以用分数 xy 表示,其中 1≤x≤n,1≤y≤m，且 x,y是整数。一个数是纯循环的，当且仅当其可以写成以下形式：a.c1˙c2c3…cp-1cp˙其中，a 是一个整数，p≥1；对于 1≤i≤p，ci是 kk 进制下的一位数字。例如，在十进制下，0.45454545……=0.4˙5˙是纯循环的，它可以用 5/11、10/22 等分数表示；在十进制下，0.1666666……=0.16˙则不是纯循环的，它可以用 1/6 等分数表示。需要特别注意的是，我们认为一个整数是纯循环的，因为它的小数部分可以表示成 0 的循环或是 k?1 的循环；而一个小数部分非 0 的有限小数不是纯循环的。

Input

只有一行，包含三个十进制数N,M,K意义如题所述,保证 1≤n≤10^9,1≤m≤10^9,2≤k≤2000

Output

一行一个整数，表示满足条件的美的数的个数。

Sample Input

2 6 10

Sample Output

4
explanation
满足条件的数分别是：
1/1=1.0000……
1/3=0.3333……
2/1=2.0000……
2/3=0.6666……
1/1 和 2/2 虽然都是纯循环小数，但因为它们相等，因此只计数一次；同样，1/3 和 2/6 也只计数一次。

正解：莫比乌斯函数+杜教筛。

首先，小学奥数学过一个东西。一个分数是纯循坏小数，当且仅当这个分数的分母没有2或5作为因子，而10分解以后正好是2和5。那么我们类比一下：k进制分数是纯循环小数是不是就是分母没有k的因数。也就是说，k和分母要互质。

所以我们要求的就是：$Ans=sum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)==1,gcd(j,k)==1]$

然后一脸懵逼，果断24分暴力。（套路用不上了。。）

24分暴力：

 1 //It is made by wfj_2048~
 2 #include <algorithm>
 3 #include <iostream>
 4 #include <complex>
 5 #include <cstring>
 6 #include <cstdlib>
 7 #include <cstdio>
 8 #include <vector>
 9 #include <cmath>
10 #include <queue>
11 #include <stack>
12 #include <map>
13 #include <set>
14 #define inf (1<<30)
15 #define il inline
16 #define RG register
17 #define ll long long
18 
19 using namespace std;
20 
21 int n,m,k,ans;
22 
23 il int gi(){
24     RG int x=0,q=1; RG char ch=getchar(); while ((ch<'0' || ch>'9') && ch!='-') ch=getchar();
25     if (ch=='-') q=-1,ch=getchar(); while (ch>='0' && ch<='9') x=x*10+ch-48,ch=getchar(); return q*x;
26 }
27 
28 il int gcd(RG int a,RG int b){ return b ? gcd(b,a%b) : a; }
29 
30 il void work(){
31     n=gi(),m=gi(),k=gi();
32     for (RG int i=1;i<=m;++i){
33     if (gcd(i,k)!=1) continue;
34     for (RG int j=1;j<=n;++j) if (gcd(i,j)==1) ans++;
35     }
36     printf("%d
",ans); return;
37 }
38 
39 int main(){
40     work();
41     return 0;
42 }

其实上式是可以化简的，我们考虑如何化简这个式子。

$Ans=sum_{i=1}^{m}[gcd(i,k)==1]sum_{j=1}^{n}[gcd(i,j)==1]$

$Ans=sum_{i=1}^{m}[gcd(i,k)==1]sum_{j=1}^{n}sum_{d|i,d|j}mu(d)$

$Ans=sum_{i=1}^{m}[gcd(i,k)==1]sum_{d|i}^{n}mu(d)left lfloor frac{n}{d} ight floor$

其实到这里，就是莫比乌斯函数的24分暴力了。。但这不是跟$O(n^{2})$暴力一样的分吗。。

我们把$mu$提前枚举，那么$Ans=sum_{d=1}^{min(n,m)} mu(d)left lfloor frac{n}{d} ight floor sum_{i=1}^{m} [d|i,gcd(i,k)==1]$

$Ans=sum_{d=1}^{min(n,m)} mu(d)left lfloor frac{n}{d} ight floor sum_{i=1}^{frac{m}{d}} [gcd(i*d,k)==1]$

$Ans=sum_{d=1}^{min(n,m)} [gcd(d,k)==1]mu(d)left lfloor frac{n}{d} ight floor sum_{i=1}^{frac{m}{d}} [gcd(i,k)==1]$

我们设$f(n)=sum_{i=1}^{n}[gcd(i,k)==1]$。根据欧几里得定理，$f(n)=left lfloor frac{n}{k} ight floor f(k)+f(n mod k)$

于是我们预处理出$f$在$[1,k]$的取值，就能把复杂度降到$O(nlogn)$，实际上可以过84分，相当于$O(n)$的复杂度。

84分暴力：

 1 //It is made by wfj_2048~
 2 #include <algorithm>
 3 #include <iostream>
 4 #include <complex>
 5 #include <cstring>
 6 #include <cstdlib>
 7 #include <cstdio>
 8 #include <vector>
 9 #include <cmath>
10 #include <queue>
11 #include <stack>
12 #include <map>
13 #include <set>
14 #define inf (1<<30)
15 #define N (20000010)
16 #define il inline
17 #define RG register
18 #define ll long long
19 
20 using namespace std;
21 
22 int f[2010],vis[N],mu[N],prime[N],n,m,k,nn,cnt;
23 ll ans;
24 
25 il int gi(){
26     RG int x=0,q=1; RG char ch=getchar(); while ((ch<'0' || ch>'9') && ch!='-') ch=getchar();
27     if (ch=='-') q=-1,ch=getchar(); while (ch>='0' && ch<='9') x=x*10+ch-48,ch=getchar(); return q*x;
28 }
29 
30 il int gcd(RG int a,RG int b){ return b ? gcd(b,a%b) : a; }
31 
32 il void pre(){
33     vis[1]=mu[1]=1;
34     for (RG int i=2;i<=nn;++i){
35     if (!vis[i]) mu[i]=-1,prime[++cnt]=i;
36     for (RG int j=1,k=i*prime[j];j<=cnt && k<=nn;++j,k=i*prime[j]){
37         vis[k]=1; if (i%prime[j]) mu[k]=-mu[i]; else break;
38     }
39     }
40     for (RG int i=1;i<=k;++i) f[i]=f[i-1]+(gcd(i,k)==1); return;
41 }
42 
43 il void work(){
44     n=gi(),m=gi(),k=gi(),nn=min(n,m); pre();
45     for (RG int d=1;d<=nn;++d){
46     if (gcd(d,k)!=1) continue;
47     ans+=(ll)mu[d]*(ll)(n/d)*(ll)((ll)(m/d/k)*(ll)f[k]+(ll)f[m/d%k]);
48     }
49     printf("%lld
",ans); return;
50 }
51 
52 int main(){
53     work();
54     return 0;
55 }

100分：

自己写不动了。。LCF学长的题解写得挺好的：http://www.cnblogs.com/lcf-2000/p/6250330.html

摘自学长博客：

我们考虑接下来该如何优化。由于$lfloor frac{m}{x} floor$只有$sqrt{m}$种取值，$lfloor frac{n}{x} floor$只有$sqrt{n}$种取值，于是我们显然可以分段求和。然后，我们就需要快速求出$sum_{i=1}^n[iperp k]mu(i)$的值。

不妨设$g(n,k)=sum_{i=1}^n[iperp k]mu(i)$，我们来考虑一下这个函数如何快速求。我们先考虑$k$的一个质因数$p$，那么$k$显然可以写成$p^cq$的形式。由于在$[1,n]$的范围中只有与$k$互质的才是有效值，那么若$xperp k$，我们可以得到$xperp p$并且$xperp q$。于是，我们可以考虑从$xperp q$的取值中减去$x$不与$p$互质的部分，就可以得到$xperp k$的部分。这里如果不懂的话可以自己画一个$xperp q$，$xperp p$，$xperp k$的关系图理解一下。

由于所有与$q$互质的数一定可以写成$p^xy(yperp q)$的形式。那么我们需要减去的数一定满足$x>0$。又由于当$x>1$时$mu(p^xy)=0$，所以我们只需要考虑$x=1$的情况即可。在这种情况下，我们需要考虑的数就是$py(yperp q)$的形式，所以我们可以得到如下式子：egin{aligned} g(n,k)&=sum_{i=1}^n[iperp q]mu(i)-sum_{y=1}^{lfloorfrac{n}{p} floor}[pyperp q]mu(py) \&=g(n,q)- sum_{y=1}^{lfloorfrac{n}{p} floor}[yperp q]mu(py)end{aligned}

上面的最后一步是由于$pperp q$，所以$pyperp q$只需在保证$yperp q$即可。

我们接着来考虑后一个式子。显然当$pperp y$的时候$mu(py)=mu(p)mu(y)$，否则$mu(py)=0$。于是，我们又得到了：

egin{aligned} g(n,k)&=g(n,q)- sum_{y=1}^{lfloorfrac{n}{p} floor}[yperp p][yperp q]mu(p)mu(y)\&=g(n,q)-mu(p)sum_{y=1}^{lfloorfrac{n}{p} floor}[yperp k]mu(y)\&=g(n,q)+g(lfloorfrac{n}{p} floor,k) end{aligned}

于是我们就可以递归求解了。容易发现边界情况就是$n=0$或者$k=1$。$n=0$的时候直接返回$0$就可以了，$k=1$的时候就是莫比乌斯函数的前缀和，用杜教筛求出来就可以了。由于每次递归要么$n$会变成$lfloor frac{n}{p} floor$，有$O(sqrt{n})$种取值；要么$p$少了一个质因数，有$omega(k)$种取值，所以总共有$O(omega(k)sqrt{n})$种取值，记忆化搜索即可。其中$omega(k)$表示$k$的不同的质因子数目。于是最后的总复杂度为$O(omega(k)sqrt{n}+n^{frac{2}{3}})$，可以通过此题。

我也不写哈希了，直接强行$map$就好。。

 1 //It is made by wfj_2048~
 2 #include <algorithm>
 3 #include <iostream>
 4 #include <complex>
 5 #include <cstring>
 6 #include <cstdlib>
 7 #include <cstdio>
 8 #include <vector>
 9 #include <cmath>
10 #include <queue>
11 #include <stack>
12 #include <map>
13 #include <set>
14 #define inf (1<<30)
15 #define N (5000010)
16 #define il inline
17 #define RG register
18 #define ll long long
19 #define calc(x) ( (x/k)*f[k]+f[x%k] )
20 
21 using namespace std;
22 
23 map <ll,map<ll,ll> >gg;
24 map <ll,ll> ff;
25 
26 ll f[2010],di[2010],vis[N],mu[N],prime[N],n,m,k,nn,NN,pos,cnt,kcnt,ans;
27 
28 il ll gi(){
29     RG ll x=0,q=1; RG char ch=getchar(); while ((ch<'0' || ch>'9') && ch!='-') ch=getchar();
30     if (ch=='-') q=-1,ch=getchar(); while (ch>='0' && ch<='9') x=x*10+ch-48,ch=getchar(); return q*x;
31 }
32 
33 il ll gcd(RG ll a,RG ll b){ return b ? gcd(b,a%b) : a; }
34 
35 il void pre(){
36     vis[1]=mu[1]=1; RG ll kk=k;
37     for (RG ll i=2;i<=nn;++i){
38     if (!vis[i]) mu[i]=-1,prime[++cnt]=i;
39     for (RG ll j=1,k=i*prime[j];j<=cnt && k<=nn;++j,k=i*prime[j]){
40         vis[k]=1; if (i%prime[j]) mu[k]=-mu[i]; else break;
41     }
42     }
43     for (RG ll i=2;i<=k;++i){
44     if (!(kk%i)) di[++kcnt]=i;
45     while (!(kk%i)) kk/=i;
46     }
47     for (RG ll i=1;i<=k;++i) f[i]=f[i-1]+(gcd(i,k)==1);
48     for (RG ll i=2;i<=nn;++i) mu[i]+=mu[i-1]; return;
49 }
50 
51 il ll du(RG ll x){
52     if (x<=nn) return mu[x];
53     if (ff[x]) return ff[x];
54     RG ll ans=1,pos;
55     for (RG ll i=2;i<=x;i=pos+1)
56     pos=x/(x/i),ans-=(pos-i+1)*du(x/i);
57     return ff[x]=ans;
58 }
59 
60 il ll g(RG ll x,RG ll y){
61     if (x<=1) return x;
62     if (!y) return du(x);
63     if (gg[x][y]) return gg[x][y];
64     return gg[x][y]=g(x,y-1)+g(x/di[y],y);
65 }
66 
67 il void work(){
68     n=gi(),m=gi(),k=gi(); NN=min(n,m);
69     nn=min(NN,5000000LL); pre();
70     for (RG ll d=1;d<=NN;d=pos+1){
71     pos=min(n/(n/d),m/(m/d));
72     ans+=(g(pos,kcnt)-g(d-1,kcnt))*(ll)(n/d)*(ll)calc(m/d);
73     }
74     printf("%lld
",ans); return;
75 }
76 
77 int main(){
78     work();
79     return 0;
80 }