[NOI2016]循环之美

Description

牛牛是一个热爱算法设计的高中生。在他设计的算法中，常常会使用带小数的数进行计算。牛牛认为，如果在 k 

进制下，一个数的小数部分是纯循环的，那么它就是美的。现在，牛牛想知道：对于已知的十进制数 n 和 m，在 

kk 进制下，有多少个数值上互不相等的纯循环小数，可以用分数 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 也只计数一次。

判断一个分数$\frac{x}{y}$是否是纯循环小数,就是一直除，直到出现了相同余数x，即
$$xk^l \equiv x (\bmod y)$$　　　　　　　　　　　　
$$l\neq 0$$
要求值不重复，所以$x \perp y$
可以推得$k^l \equiv 1(\bmod y)$，条件也就是$y \perp k$
答案就是\begin{aligned}&\\&\sum_{y=1}^{m}[y\perp k]\sum_{x=1}^{n}[x\perp y]\end{aligned}
莫比乌斯反演得
\begin{aligned}&\sum_{y=1}^{m}[y\perp k] \sum_{x=1}^{n}\sum_{d|x,d|y}\mu(d)\\=&\sum_{y=1}^{m}[y\perp k]\sum_{d|y}^{n}\mu(d)\lfloor \frac{n}{d}\rfloor \end{aligned}
交换顺序
\begin{aligned}&\sum_{d=1}^{n}\mu(d)\lfloor \frac{n}{d} \rfloor \sum_{y=1}^{m}[d\ |\ y][y\perp k]\\=&\sum_{d=1}^{n}\mu(d)\lfloor \frac{n}{d}\rfloor \sum_{i=1}^{\lfloor \frac{m}{d} \rfloor}[id\perp k]\\=& \sum_{d=1}^{n}[d\perp k]\mu(d)\lfloor \frac{n}{d}\rfloor \sum_{i=1}^{\lfloor \frac{m}{d} \rfloor}[i\perp k]\end{aligned}
令f[n]=$\sum_{d=1}^{n}[i\perp k]$
由同余性质：
gcd(a,b)=gcd(a%b,b)
$f(n)=\left \lfloor \frac{n}{k} \right \rfloor f(k)+f(n \bmod k)$
这样84分就有了O(n)计算就行
设$g(n,k)=\sum_{i=1}^n[i\perp k]\mu(i)$
考虑将k分解成pcq （p是素数）
于是可以将只与q互质的答案减去不与p互质的答案，不与素数p互质就是p的倍数
\begin{aligned}  g(n,k)&=\sum_{i=1}^n[i\perp q]\mu(i)-\sum_{y=1}^{\lfloor\frac{n}{p}\rfloor}[py\perp q]\mu(py) \\&=g(n,q)- \sum_{y=1}^{\lfloor\frac{n}{p}\rfloor}[y\perp q]\mu(py)\end{aligned}
显然当$p\perp y$的时候$\mu(py)=\mu(p)\mu(y)$，否则$\mu(py)=0$
还有$\mu(p)=-1$
\begin{aligned} g(n,k)&=g(n,q)- \sum_{y=1}^{\lfloor\frac{n}{p}\rfloor}[y\perp p][y\perp q]\mu(p)\mu(y)\\&=g(n,q)-\mu(p)\sum_{y=1}^{\lfloor\frac{n}{p}\rfloor}[y\perp k]\mu(y)\\&=g(n,q)+g(\lfloor\frac{n}{p}\rfloor,k) \end{aligned}
递归求解容易发现边界情况就是$n=0$或者$k=1$。$n=0$的时候直接返回$0$就可以了，$k=1$的时候就是莫比乌斯函数的前缀和，用杜教筛求出
递归时记忆化，直接用map
由于每次递归要么$n$会变成$\lfloor \frac{n}{p} \rfloor$，有$O(\sqrt{n})$种取值；要么$p$少了一个质因数，有$\omega(k)$种取值，所以总共有$O(\omega(k)\sqrt{n})$种取值，记忆化搜索即可。其中$\omega(k)$表示$k$的不同的质因子数目。于是最后的总复杂度为$O(\omega(k)\sqrt{n}+n^{\frac{2}{3}})$

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<map>
using namespace std;
typedef long long lol;
lol mu[8000001],prime[4000001];
lol tot,N=8000000,cnt;
lol d[10001],f[2001],ans;
bool vis[8000001];
map<lol,lol>M1;
map<lol,map<lol,lol> >M2;
int gcd(int a,int b)
{
  if (!b) return a;
  return gcd(b,a%b);
}
void pre(lol k)
{lol i,j;
  mu[1]=1;
  for (i=2;i<=N;i++)
    {
      if (vis[i]==0)
    {
      prime[++tot]=i;
      mu[i]=-1;
    }
      for (j=1;j<=tot;j++)
    {
      if (i*prime[j]>N) break;
      vis[i*prime[j]]=1;
      if (i%prime[j]==0)
        {
          mu[i*prime[j]]=0;
          break;
        }
      else mu[i*prime[j]]=-mu[i];
    }
    }
  for (i=1;i<=N;i++)
    mu[i]+=mu[i-1];
  lol kk=k;
  for (i=2;i<=k;i++)
    {
      if (kk%i==0) d[++cnt]=i;
      while (kk%i==0) kk/=i;
    }
}
lol get_mu(lol n)
{lol i,pos;
  if (n<=N) return mu[n];
  if (M1[n]) return M1[n];
  lol s=1;
  for (i=2;i<=n;i=pos+1)
    {
      pos=n/(n/i);
      s-=(pos-i+1)*get_mu(n/i);
    }
  return M1[n]=s;
}
lol get_s(lol n,lol k)
{
  if (n<=1) return n;
  if (k==0) return get_mu(n);
  if (M2[n][k]) return M2[n][k];
  return M2[n][k]=get_s(n,k-1)+get_s(n/d[k],k);
}
lol get_f(lol n,lol k)
{
  return (n/k)*f[k]+f[n%k];
}
int main()
{lol n,m,k;
  lol i,pos;
  cin>>n>>m>>k;
  for (i=1;i<=k;i++)
    {
      if (gcd(i,k)==1) f[i]=1;
      f[i]+=f[i-1];
    }
  N=min(N,max(n,m));
  pre(k);
  for (i=1;i<=min(n,m);i=pos+1)
    {
      pos=min(n/(n/i),m/(m/i));
      ans+=(n/i)*(get_s(pos,cnt)-get_s(i-1,cnt))*get_f(m/i,k);
    }
  cout<<ans;
}