【BZOJ4652】循环之美（NOI2016）-杜教筛

测试地址：循环之美
做法：本题需要用到杜教筛。
首先，思考一下纯循环小数这个条件有什么意义。注意到，如果一个数w$w$的循环节长度为x$x$，那么我们用kxw$k^{x}w$减去w$w$可以得到一个整数u$u$，那么w$w$就可以表示成ukx−1$\frac{u}{k^x-1}$，注意到分数能表达成这种形式是满足题目条件的充要条件，那么条件就变成一个分数能不能表示成这种形式，也就是存不存在分母w|kx−1$w|k^x-1$，x≥1$x\ge 1$的情况。
注意到w|kx−1$w|k^x-1$，也就意味着kxmodw=1$k^x\mathrm{mod}w=1$，注意到这个东西可能成立当且仅当gcd(w,k)=1$gcd(w,k)=1$，即w,k$w,k$互质，那么我们就把题目条件转化成了一个可以列式的条件了。
综上所述，这题要求的式子就是：
ans=∑ni=1∑mj=1[gcd(i,j)=1][gcd(j,k)=1]$ans=\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=1][gcd(j,k)=1]$
接下来有几种推式子的方法，由于我太菜了，只会比较简单的不能通过全部测试点的方法（在BZOJ的总时限下勉强能过），推法如下：
常规把[gcd(j,k)=1]$[gcd(j,k)=1]$处理掉，得到：
ans=∑d|kμ(d)∑ni=1∑d|j[gcd(i,j)=1]$ans=\sum _{d|k}\mu (d)\sum _{i=1}^n\sum _{d|j}[gcd(i,j)=1]$
把用jd$\frac{j}{d}$替换j$j$，得：
ans=∑d|kμ(d)∑ni=1∑⌊md⌋j=1[gcd(i,j)=1][gcd(i,d)=1]$ans=\sum _{d|k}\mu (d)\sum _{i=1}^n\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }[gcd(i,j)=1][gcd(i,d)=1]$
我们发现后面式子的形式和我们一开始列出的式子非常相似，那么我们令ans=g(n,m,k)$ans=g(n,m,k)$，那么有：
g(n,m,k)=∑d|kμ(d)g(⌊md⌋,n,d)$g(n,m,k)=\sum _{d|k}\mu (d)g(\lfloor \frac{m}{d}\rfloor ,n,d)$
边界为g(0,m,k)=g(n,0,k)=0$g(0,m,k)=g(n,0,k)=0$，而g(n,m,1)=∑ni=1∑mj=1[gcd(i,j)=1]$g(n,m,1)=\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=1]$这个式子我们在其他题目里推了好多遍了，直接杜教筛处理即可，这样我们就可以递归计算ans$ans$了。这个方法可以拿到很高的分数：92分，但不幸的是，由于这个方法不是很优或者本人太菜，这里给的代码并不能过最大的两个点，有待学习。
以下是本人代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll hashsiz=7000003;
ll n,m,k,limit;
ll prime[1100010],mu[1100010],summu[1100010];
ll hashlist[7000010]={0},hashval[7000010],inhashlist=0;
bool vis[1100010]={0};

void calc_mu()
{
    mu[1]=1;
    prime[0]=0;
    for(ll i=2;i<=limit;i++)
    {
        if (!vis[i])
        {
            prime[++prime[0]]=i;
            mu[i]=-1;
        }
        for(ll j=1;j<=prime[0]&&i*prime[j]<=limit;j++)
        {
            vis[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
    }
    summu[0]=0;
    for(int i=1;i<=limit;i++)
        summu[i]=summu[i-1]+mu[i];
}

ll hashfind(ll x)
{
    ll pos=x%hashsiz;
    while(hashlist[pos]&&hashlist[pos]!=x)
    {
        pos++;
        if (pos>=hashsiz) pos-=hashsiz;
    }
    if (hashlist[pos]==x) return pos;
    else return -1;
}

void hashinsert(ll x,ll val)
{
    if (inhashlist>=hashsiz) return;
    inhashlist++;
    ll pos=x%hashsiz;
    while(hashlist[pos]&&hashlist[pos]!=x)
    {
        pos++;
        if (pos>=hashsiz) pos-=hashsiz;
    }
    hashlist[pos]=x;
    hashval[pos]=val;
}

ll sum_mu(ll n)
{
    if (n<=limit) return summu[n];
    int pos=hashfind(n);
    if (pos!=-1) return hashval[pos];
    ll ans=1;
    for(ll i=n;i>=2;i=n/(n/i+1))
    {
        ll l=max(n/(n/i+1)+1,2ll),r=i;
        ans-=(r-l+1)*sum_mu(n/i);
    }
    hashinsert(n,ans);
    return ans;
}

ll sum(ll n,ll m)
{
    ll ans=0;
    for(ll i=min(n,m);i>=1;i=max(n/(n/i+1),m/(m/i+1)))
    {
        ll l=max(n/(n/i+1),m/(m/i+1))+1,r=i;
        ans+=(sum_mu(r)-sum_mu(l-1))*(n/i)*(m/i);
    }
    return ans;
}

ll calc(ll n,ll m,ll k)
{
    ll ans=0;
    if (!n||!m) return 0;
    if (k==1) return sum(n,m);
    for(ll i=1;i*i<=k;i++)
        if (k%i==0)
        {
            ans+=mu[i]*calc(m/i,n,i);
            if (i*i!=k) ans+=mu[k/i]*calc(m/(k/i),n,k/i);
        }
    return ans;
}

int main()
{
    scanf("%lld%lld%lld",&n,&m,&k);
    limit=1;
    while(limit*limit*limit<=max(n,m)) limit++;
    limit*=limit;
    limit=max(limit,k);
    calc_mu();
    printf("%lld",calc(n,m,k));

    return 0;
}