bzoj4652 [Noi2016]循环之美

http://www.elijahqi.win/archives/3364
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 也只计数一次。
HINT

Source
感觉很毒瘤的一题
首先考虑到底如何才会形成纯循环小数 可以知道 商循环的时候不一定是纯循环小数 只有当余数循环的时候才是纯循环小数 那么所以如何快速的求出第几位之后的余数那么就是
设该小数的分数形式为$\frac{x}{y}$那么$x*k^i\%y$就是我们想要的第i位上的商是多少了
由题面提示所知题中求的纯循环小数一定是第二次出现相同的数字的时候是和第一次的整数相同的数
那么设纯循环小数的循环节长度为L 那么可以知道 $x*k^l \equiv x (mod\ y)$由于gcd(x,y)==1
所以$k^l\equiv 1(mod\ y)$ ->gcd(k,y)=1
那么于是原题转化为$\sum\limits_{x=1}^{n}\sum\limits_{y=1}^{m} [gcd(x,y)==1]\times[gcd(y,k)==1]$
$\sum\limits_{y=1}^{m}[gcd(y,k)==1]\sum\limits_{x=1}^{n}[gcd(x,y)==1]$
考虑后面的一步怎么求 用$\varepsilon$替换一下变成
$\sum\limits_{y=1}^{m}[gcd(y,k)==1]\sum\limits_{x=1}^{n}\sum\limits_{d|x,d|y} \mu(d)$
考虑把d移动到前面去 变成
$\sum\limits_{d=1}^{min(n,m)}\mu(d)[gcd(d,k)==1]\sum\limits_{d|x}^{n}\sum\limits_{d|y}^{m}[gcd(y,k)==1]$
$\sum\limits_{d=1}^{min(n,m)}\mu(d)[gcd(d,k)==1]\sum\limits_{x=1}^{\frac{n}{d}}\sum\limits_{y=1}^{\frac{m}{d}}[gcd(y,k)==1]$
$\sum\limits_{d=1}^{min(n,m)}\mu(d)[gcd(d,k)==1]\lfloor\frac{n}{d}\rfloor\sum\limits_{y=1}^{\frac{m}{d}}[gcd(y,k)==1]$
考虑现在我们已经将式子分成了两部分 后半部分
$\sum\limits_{y=1}^{\frac{m}{d}}[gcd(y,k)==1]$考虑针对因为其中一个固定是k那么显然这个式子呈现k个一循环的样子即当$y<=k且gcd(y,k)==1$时那么即使y+=k最后的答案gcd(y,k)也还是1 于是我们可以考虑o(k) 预处理所有的gcd的和然后询问的时候整块的就直接乘法 零散的就分开即可 设后一段为函数$f[x]$
那么$f[x]=\lfloor\frac{x}{k}\rfloor\times f[k]+f[x\%k]$
那么每次后面的这段式子都可以o(1)回答了 那么前面的怎么办
考虑也是一样先化简
$\sum\limits_{i=1}^{n}\mu(i)[gcd(i,k)==1]$
$\sum\limits_{i=1}^{n}\mu(i)\sum\limits_{d|i,d|k}\mu(d)$
$\sum\limits_{d|k}\mu(d) \sum\limits_{d|i}^{i<=n}\mu(i)$
$\sum\limits_{d|k}\mu(d) \sum\limits_{i=1}^{\frac{n}{d}}\mu(i\times d)$
考虑gcd(i,d)==1的情况 如果不==1那么对答案也不会有贡献
所以$\sum\limits_{d|k}\mu(d)\sum\limits_{i=1}^{\frac{n}{d}}\mu(i)\times \mu(d)$
$\sum\limits_{d|k}\mu(d)^2\sum\limits_{i=1}^{\frac{n}{d}}\mu(i)[gcd(i,d)==1]$
好了这个式子现在可以考虑递归做下去了 如果k=1那么显然可以杜教筛搞定
如果k!=1那么递归做下去即可 然后用hash把同一个x,k记忆化即可

#include<queue>
#include<cstdio>
#include<cctype>
#include<algorithm>
#define ll long long
using namespace std;
inline char gc(){
    static char now[1<<16],*S,*T;
    if (T==S){T=(S=now)+fread(now,1,1<<16,stdin);if (T==S) return EOF;}
    return *S++;
}
inline int read(){
    int x=0,f=1;char ch=gc();
    while(!isdigit(ch)) {if (ch=='-') f=-1;ch=gc();}
    while(isdigit(ch)) x=x*10+ch-'0',ch=gc();
    return x*f;
}
const int N=1e6+10;
const int mod=1e6+3;
int n,m,k,sum[N],prime[N],tot,f[N],mu[N],h[N],num;
ll ans;bool not_prime[N];
inline int gcd(int x,int y){return !y?x:gcd(y,x%y);}
inline ll calc(int x){return (ll)x/k*f[k]+f[x%k];}
inline void init(){
    for (int i=1;i<=k;++i) f[i]=f[i-1]+(gcd(i,k)==1);mu[1]=1;
    for (int i=2;i<=1e6;++i){
        if (!not_prime[i]) prime[++tot]=i,mu[i]=-1;
        for (int j=1;prime[j]*i<=1e6;++j){
            not_prime[prime[j]*i]=1;
            if (i%prime[j]==0){mu[prime[j]*i]=0;break;}
            else mu[prime[j]*i]=-mu[i];
        }
    }for (int i=1;i<=1e6;++i) sum[i]=sum[i-1]+mu[i];
}
struct node{
    int x,k,next,v;
}data[6*N];
inline void insert1(int x,ll k,int v){
    int kk=k;k=k*1e9+x;static int tmp;tmp=k%mod;
    for (int i=h[tmp];i;i=data[i].next){
        if (data[i].x==x&&data[i].k==kk) return;
    }data[++num].next=h[tmp];h[tmp]=num;data[num].x=x;data[num].k=kk;data[num].v=v;
}
inline int query(int x,ll k){
    int kk=k;k=k*1e9+x;static int tmp;tmp=k%mod;
    for (int i=h[tmp];i;i=data[i].next){
        if(x==data[i].x&&kk==data[i].k) return data[i].v;
    }return -1;
}
inline int gao(int x,int k){
    if(k==1&&x<=1e6) return sum[x];if(!x) return 0;
    int tmp=query(x,k);if (tmp!=-1) return tmp;tmp=0;
    if (k==1){tmp=1;int last;
        for (int i=2;i<=x;i=last+1) last=x/(x/i),tmp-=(last-i+1)*gao(x/i,1);
    }else{
        for (int i=1;i*i<=k;++i){
            if (k%i==0){
                if (mu[i]) tmp+=mu[i]*mu[i]*gao(x/i,i);
                if (i*i!=k&&mu[k/i]) tmp+=mu[k/i]*mu[k/i]*gao(x/(k/i),k/i);
            }
        }
    }insert1(x,k,tmp);return tmp;
}
int main(){
    freopen("bzoj4652.in","r",stdin);
    n=read();m=read();k=read();init();
    int nn=min(n,m),last;int lst=0,now=0;
    for (int i=1;i<=nn;i=last+1){
        last=min(n/(n/i),m/(m/i));now=gao(last,k);
        ans+=(now-lst)*(n/i)*calc(m/i);lst=now;
    }printf("%lld\n",ans);//printf("%d\n",num);
    return 0;
}