# bzoj 4652: [Noi2016]循环之美

（公式怎么打啊。。。）

 1 #include<bits/stdc++.h>
2 #define LL long long
3 using namespace std;
4
5 const int maxn=5000010;
6 const int mod=100007;
7 const int maxk=2010;
8
9 int n,m,k,prime_k[maxk],cnt;
10 int prime[maxn],f[maxk];
12 int h1[11][mod],n1[maxn],t1[maxn],t2;
13 bool vis[maxn],huzhi[maxk];
14 LL ans,zto[maxn],zt1[maxn],mo[maxn];
15
16 void mobius()
17 {
18     mo[1]=1;
19     for (int i=2; i<maxn; i++)
20     {
21         if (!vis[i]) prime[++prime[0]]=i,mo[i]=-1;
22         for (int j=1; j<=prime[0] && i*prime[j]<maxn; j++)
23         {
24             vis[i*prime[j]]=1;
25             if (i%prime[j]==0)
26             {
27                 mo[i*prime[j]]=0;
28                 break;
29             }
30             mo[i*prime[j]]=-mo[i];
31         }
32     }
33     for (int i=2; i<maxn; i++) mo[i]+=mo[i-1];
34     for (int i=1; i<=k; i++) f[i]=f[i-1]+(__gcd(i,k)==1);
35 }
36
37 LL F(int x){return (x/k)*f[k]+f[x%k];}
38
39 LL solve_mo(int x)
40 {
41     if (x<maxn) return mo[x];
42     for (int i=head[x%mod];i;i=next[i]) if (to[i]==x) return zto[i];
44     for (int i=2,nt=0; nt<x; i=nt+1) nt=x/(x/i),zto[now]-=(nt-i+1)*solve_mo(x/i);
45     return zto[now];
46 }
47 LL G(int x, int y)
48 {
49     if (!x) return solve_mo(y);
50     if (y<=1) return y;
51     for (int i=h1[x][y%mod];i;i=n1[i]) if (t1[i]==y) return zt1[i];
52     int now=++t2; t1[t2]=y; n1[t2]=h1[x][y%mod]; h1[x][y%mod]=t2;
53     return zt1[now]=G(x-1,y)+G(x,y/prime_k[x]);
54 }
55 int main()
56 {
57     scanf("%d%d%d",&n,&m,&k); mobius();
58     for (int i=1; prime[i]<=k; i++) if (k%prime[i]==0) prime_k[++cnt]=prime[i];
59     for (int d=1,l=min(n,m),pos; d<=l; d=pos+1)
60     {
61         pos=min(n/(n/d),m/(m/d));// cout<<"!";
62         ans+=(G(cnt,pos)-G(cnt,d-1))*(LL)(n/d)*F(m/d);
63     }
64     printf("%lld\n",ans);
65     return 0;
66 }

posted @ 2017-04-12 21:26  ws_ccd  阅读(219)  评论(0编辑  收藏  举报