[NOI2016]循环之美

题目

这个题让我们求的就是这个柿子

\[\sum_{i=1}^n\sum_{j=1}^m[i\perp j][j\perp k] \]

非常简单啊，\(i\perp j\)保证了\(\frac{i}{j}\)是一个最简分数，\(j\perp k\)保证\(\frac{1}{j}\)是一个无限不循环小数

于是我们要求的就是上面那个柿子啦

首先先看一个简单的问题，求\(1\)到\(n\)中和\(k\)互质的数的个数

我们设

\[F(d)=\sum_{i=1}^n[d|(i,k)]=[d|k]\times \left \lfloor \frac{n}{d} \right \rfloor \]

反演可得，我们要求的东西就是

\[\sum_{d|k}\mu(d)\left \lfloor \frac{n}{d} \right \rfloor \]

我们可以在\(\sigma(k)\)的时间内求出这个东西

现在我们可以来求我们的柿子啦

交换一下求和符号

\[\sum_{j=1}^m[j\perp k]\sum_{i=1}^n[i\perp j] \]

于是后面可以写成

\[\sum_{j=1}^m[j\perp k]\sum_{d|j}\mu(d)\left \lfloor \frac{n}{d} \right \rfloor \]

到前面来枚举\(d\)

\[\sum_{d=1}^n\left \lfloor \frac{n}{d} \right \rfloor\mu(d)\sum_{d|j,j<=m}[j\perp k] \]

\[=\sum_{d=1}^n\left \lfloor \frac{n}{d} \right \rfloor\mu(d)\sum_{t=1}^{\left \lfloor \frac{m}{d} \right \rfloor}[td\perp k] \]

我们发现\(td\)和\(k\)互质完全可以表示成\(t\perp k\)且\(d\perp k\)，于是我们又能把柿子写成这个样子

\[\sum_{d=1}^n\left \lfloor \frac{n}{d} \right \rfloor S(\left \lfloor \frac{m}{d} \right \rfloor)\mu(d)[d\perp k] \]

其中\(S(n)=\sum_{i=1}^n[i\perp k]\)

我们发现我们现在只需要求出\(\mu(d)[d\perp k]\)的前缀和就可以整除分块啦

数据范围不小啊，考虑杜教筛

。。。之后就考虑不出来了，之后成爷随便一化就秒掉了

我们设\(f(i)=\mu(i)[i\perp k]\)，\(g(i)=[i\perp k]\)

我们考虑一下\(f\times g\)

\[(f\times g)(n)=\sum_{d|n}[d\perp k][\frac{n}{d}\perp k]\mu(d) \]

我们发现只有当\(n\perp k\)的时候才能保证\(d\perp k\)和\(\frac{n}{d}\perp k\)同时满足，于是我们可以直接把\([d\perp k][\frac{n}{d}\perp k]\)变成\([n\perp k]\)提到外面来，变成

\[(f\times g)(n)=[n\perp k]\sum_{d|n}\mu(d)=[n=1] \]

至于\(g(i)\)的前缀和，发现那不就是\(S(i)\)吗

于是我们现在就可以杜教筛啦，复杂度是\(O(\sigma(k)n^{\frac{2}{3}})\)

不过有一点需要注意，就是当同时对\(n,m\)整除分块的时候不能使用类似于min_25筛的那一个\(trick\)，只能使用\(map\)或\(hash\)来记忆化

代码

#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<tr1/unordered_map>
#define re register
#define LL long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
using namespace std::tr1;
const int maxn=2000005;
int f[maxn],p[maxn>>1],mu[maxn],pre[maxn];
int d[2005],n,m,k,H,M,Sqr,tot;
unordered_map<int,int> vis,ma,F;
inline int calc(int n) {
	if(F[n]) return F[n];
    int ans=0;
    for(re int i=1;i<=tot;i++) 
        ans+=mu[d[i]]*(n/d[i]);
    return F[n]=ans;
}
inline int gcd(int a,int b) {return !b?a:gcd(b,a%b);}
inline int solve(int n) {
    if(n<=M) return pre[n];
    if(vis[n]) return ma[n];
    int ans=1;
    for(re int l=2,r;l<=n;l=r+1) {
        r=n/(n/l);
        ans-=solve(n/l)*(calc(r)-calc(l-1));
    }
    vis[n]=1,ma[n]=ans;
    return ans;
}
int main() {
    scanf("%d%d%d",&n,&m,&k);M=min(max(n,m),2000000);
    f[1]=1,mu[1]=1;
    for(re int i=2;i<=M;i++) {
        if(!f[i]) p[++p[0]]=i,mu[i]=-1;
        for(re int j=1;j<=p[0]&&p[j]*i<=M;j++) {
            f[p[j]*i]=1;
            if(i%p[j]==0) break;
            mu[p[j]*i]=-1*mu[i];
        }
    }
    for(re int i=1;i<=k;i++) {
        if(k%i||!mu[i]) continue;
        d[++tot]=i; 
    }
    for(re int i=1;i<=M;i++) 
        if(gcd(i,k)!=1) pre[i]=pre[i-1];
            else pre[i]=pre[i-1]+mu[i];
    LL ans=0;
    H=min(n,m);
    for(re int l=1,r;l<=H;l=r+1) {
        r=min(n/(n/l),m/(m/l));
        ans+=1ll*(n/l)*calc(m/l)*(solve(r)-solve(l-1));
    }
    printf("%lld\n",ans);
    return 0;
}