bzoj 4652: [Noi2016]循环之美

额，，网上一堆题解，，随便一找就找到笨蒟蒻扒的了。

这个比较神奇的是纯循环小数就是[(y,k)=1]，题解有证明这个的，貌似就是k进制下的类似循环节，不会不会。。

然后这道题就变成了求这个东西：∑(x<=n)∑(y<=m) [(x,y)=1][(y,k)=1]

显然要反演一下，出来这个东西∑（d<=n）[(d,k)==1] μ(d) (n/d) ∑ (j<=m/d) [(j,k)==1]

然后k固定，所以搞一个∑（d<=n）[(d,k)==1] μ(d)和 ∑ (j<=m/d) [(j,k)==1]，这两部分有神奇的递推式，然后各种各种，，不会不会。。。

（公式怎么打啊。。。）

这位神犇写的很详细 %%%% http://www.cnblogs.com/lcf-2000/p/6250330.html

 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];
11 int head[mod],next[maxn],to[maxn],tt;
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];
43     int now=++tt; to[tt]=x; next[tt]=head[x%mod]; head[x%mod]=tt; zto[now]=1;
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
",ans);
65     return 0;
66 }

小技巧：这里的记忆化搜索用链表代替了map（会不会map也是差不多的实现原理啊2333）