luogu P1587 [NOI2016]循环之美

传送门

首先要知道什么样的数才是"纯循环数".打表可以发现,这样的数当且仅当分母(y)和(k)互质,这是因为,首先考虑除法过程,每次先给当前余数(*k),然后对分母做带余除法,如果把每次的余数写成一个序列,那么"纯循环数"就要使的第一个出现过两次的余数正好为序列中第一个余数.设第一个余数为(x),循环节长度为(a),那么会有(xk^aequiv x pmod y),即(k^aequiv 1 pmod y),注意到(ige 1)时得到的(k^i mod y)都是(gcd(k,y))的倍数,所以必须要有(gcd(k,y)=1)

然后要使得数值不同,所以其实要求的是这个东西

(sum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)=1][gcd(j,k)=1])

先把只有(j)和(k)的提前

(sum_{j=1}^{m}[gcd(j,k)=1]sum_{i=1}^{n}[gcd(i,j)=1])

然后把([gcd(j,k)=1])转化一下

(sum_{j=1}^{m}sum_{d|j,d|k}mu(d)sum_{i=1}^{n}[gcd(i,j)=1])

把(d)提前

(sum_{d=1}^{m}mu(d)sum_{j=1}^{lfloorfrac{m}{d} floor}sum_{i=1}^{n}[gcd(i,jd)=1])

(sum_{d=1}^{m}mu(d)sum_{j=1}^{lfloorfrac{m}{d} floor}sum_{i=1}^{n}[gcd(i,j)=1][gcd(i,d)=1])

后面那个是不是有点眼熟?我们如果记(s(i,j,k)=sum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)=1][gcd(j,k)=1]),那么就能得到

(s(i,j,k)=sum_{d|k}mu(d)s(lfloorfrac{m}{d} floor,n,d))

这个递归处理就好了,边界就是(m=0)或者(n=0)时为(0),(k=1)时为(sum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)=1]),数论分块+杜教筛即可.复杂度大概是(O(log nlog m sqrt{n}+n^{frac{2}{3}})),这个还是比较慢的,加了记忆化都要1800多ms 也可能是我写丑了

// luogu-judger-enable-o2
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<vector>
#include<cmath>
#include<ctime>
#include<queue>
#include<map>
#include<set>
#define LL long long
#define db double

using namespace std;
const int N=5e6+10,M=2000+10;
int rd()
{
    int x=0,w=1;char ch=0;
    while(ch<'0'||ch>'9'){if(ch=='-') w=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+(ch^48);ch=getchar();}
    return x*w;
}
int n,m,k,prm[N],tt;
bool pp[N];
LL mu[N];
map<int,LL> f;
int lm;
LL siv(int nn)
{
    if(nn<=N-10) return mu[nn];
    if(f.count(nn)) return f[nn];
    LL &an=f[nn];
    an=1;
    for(int i=2,j;i<=nn;i=j+1)
    {
        j=nn/(nn/i);
        an-=1ll*siv(nn/i)*(j-i+1);
    }
    return an;
}
LL ff(int nn,int mm)
{
    lm=min(nn,mm);
    LL an=0;
    for(int i=1,j;i<=lm;i=j+1)
    {
        j=min(nn/(nn/i),mm/(mm/i));
        an+=1ll*(siv(j)-siv(i-1))*(nn/i)*(mm/i);
    }
    return an;
}
vector<int> dd[M];
struct node
{
    int n,m,k;
    bool operator < (const node &bb) const {return n!=bb.n?n<bb.n:(m!=bb.m?m<bb.m:k<bb.k);}
};
map<node,LL> g;
LL sov(int n,int m,int k)
{
    if(n<=0||m<=0) return 0;
    node nw=(node){n,m,k};
    if(g.count(nw)) return g[nw];
    if(k==1) return g[nw]=ff(n,m);
    LL an=0;
    vector<int>::iterator it;
    for(it=dd[k].begin();it!=dd[k].end();++it)
    {
        int i=*it;
        an+=1ll*sov(m/i,n,i)*(mu[i]-mu[i-1]);
    }
    return g[nw]=an;
}

int main()
{
    mu[1]=1;
    for(int i=2;i<=N-10;++i)
    {
        if(!pp[i]) prm[++tt]=i,mu[i]=-1;
        for(int j=1;j<=tt&&i*prm[j]<=N-10;++j)
        {
            pp[i*prm[j]]=1;
            if(i%prm[j]==0) break;
            mu[i*prm[j]]=-mu[i];
        }
    }
    for(int i=2;i<=N-10;++i) mu[i]+=mu[i-1];
    n=rd(),m=rd(),k=rd();
    for(int i=1;i<=k;++i)
        if(k%i==0) dd[k].push_back(i);
    int nn=dd[k].size();
    for(int i=0;i<nn-1;++i)
    {
        int x=dd[k][i];
        vector<int>::iterator it;
        for(it=dd[k].begin();it!=dd[k].end();++it)
        {
            int y=*it;
            if(x<y) break;
            if(x%y==0) dd[x].push_back(y);
        }
    }
    printf("%lld
",sov(n,m,k));
    return 0;
}