UOJ221 循环之美 【莫比乌斯反演】【分块】【杜教筛】

题目链接：http://uoj.ac/problem/221

题意：求有多少对$(x,y)$，满足$1\le x\le n$,$1\le y\le m$，$\frac{x}{y}$在$k$进制下是纯循环小数，且$\frac{x}{y}$数值上两两不同

题解：
有一个很重要的结论：$\frac{x}{y}$在$k$进制下为纯循环小数当且仅当$gcd(y,k)==1$，证明见 https://www.cnblogs.com/lcf-2000/p/6250330.html
所以$$ans={\Sigma }_{x=1}^n{\Sigma }_{y=1}^m[gcd(x,y)==1][gcd(y,k)==1]$$
将$[gcd(x,y)==1]$部分反演可得
$$ans={\Sigma }_{x=1}^n{\Sigma }_{y=1}^m[gcd(y,k)==1]{\Sigma }_{d|x,d|y}\mu (d)$$
$$={\Sigma }_{d=1}^{min(n,m)}\mu (d)[gcd(d,k)==1]\lfloor \frac{n}{d}\rfloor {\Sigma }_{y=1}^{\lfloor \frac{m}{d}\rfloor }[gcd(y,k)==1]$$
令$f(n,k)={\Sigma }_{i=1}^n[gcd(n,k)==1]$，那么显然
$f(n,k)=\lfloor \frac{n}{k}\rfloor f(k,k)+f(n\mathrm{m}\mathrm{o}\mathrm{d}k,k)$
$f(n,k)(n\le k)$可以O(k)预处理，对于每次询问就是$O(nlogn)$（gcd带一个log），可以获得84pts
84pts代码（写起来很简单）：

// by Balloons
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define mpr make_pair
#define debug() puts("okkkkkkkk")
#define rep(i,a,b) for(int (i)=(a);(i)<=(b);(i)++)

using namespace std;

typedef long long LL;

const int inf = 1 << 30;
const int maxn=2e7+5;

int n,m,k;
int notpm[maxn],pm[maxn],pcnt=0,f[maxn];
int mu[maxn];

void xxs(){
	notpm[1]=1;mu[1]=1;
	for(int i=2;i<=maxn-5;i++){
		if(!notpm[i]){
			pm[++pcnt]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=pcnt && i*pm[j]<=maxn-5;j++){
			notpm[i*pm[j]]=1;
			if(i%pm[j]==0){
				mu[i*pm[j]]=0;
				break;
			}
			mu[i*pm[j]]=-mu[i];
		}
	}
}

int gcd(int a,int b){return !b?a:gcd(b,a%b);}

int main(){
	scanf("%d%d%d",&n,&m,&k);
	xxs();
	for(int i=1;i<=k;i++){
		f[i]=f[i-1]+(gcd(i,k)==1);
	}
	LL ans=0;
	for(int i=1;i<=n;i++){
		if(gcd(i,k)!=1)continue;
		ans+=1ll * mu[i] * n/i * ((m/i)/k * f[k] + f[m/i%k]);
	}
	printf("%lld\n",ans);

	return 0;
}

继续考虑优化式子：
令$g(n,k)={\Sigma }_{d=1}^n\mu (d)[gcd(d,k)==1]$
考虑 $k$的一个质因子$q$，则显然$k=p^c\ast q((q,k)==1)$
因为$[1,n]$中与$k$互质的数才能加入$g$中，所以考虑如何求出与$k$互质的数，显然等价于求$p,q$与$k$互质的数
容斥原理易得答案即为$(x,q)==1$的个数减去$(x,p)̸=1$的个数$1\le x\le n$
因为与$q$互质的数可以写成$p^{c}y((y,q)==1,(p,q)==1)$的形式，而若$c\ge 2$，那么$\mu (p^{c}y)$一定为0，不加入答案，所以$c==1$，即所有与$q$互质的数一定可以写成$py$的形式，根据上面的容斥原理我们可以列出一个式子：
$$g(n,k)={\Sigma }_{i=1}^n[gcd(i,q)==1]\mu (i)-{\Sigma }_{y=1}^{\lfloor \frac{n}{p}\rfloor }[(gcd(py,q)==1]\mu (py)$$
$$=g(n,q)-{\Sigma }_{y=1}^{\lfloor \frac{n}{p}\rfloor }[gcd(y,q)==1]\mu (py)$$
根据积性函数性质$\mu (py)=\mu (p)\mu (y)[gcd(p,y)==1]$
所以原式$$g(n,k)=g(n,q)-{\Sigma }_{y=1}^{\lfloor \frac{n}{p}\rfloor }[gcd(y,p)==1][gcd(y,q)==1]\mu (p)\mu (y)$$
$$=g(n,q)-\mu (p){\Sigma }_{y=1}^{\lfloor \frac{n}{p}\rfloor }[y\perp k]\mu (y)$$
$$=g(n,q)+g(\lfloor \frac{n}{p}\rfloor ,k)$$
$f,g$均可以递归求解
为了方便，$g(n,k)$在求解的过程中$k$代表的是剩余的质因数个数
$g$递归的过程中是用了类似前向星的方式建立分层关系的

// by Balloons
#include <cstdio>
#include <cstring>
#include <map>
#include <iostream>
#include <algorithm>
#define mpr make_pair
#define debug() puts("okkkkkkkk")
#define rep(i,a,b) for(int (i)=(a);(i)<=(b);(i)++)

using namespace std;

typedef long long LL;
#define int LL

const int inf = 1 << 30,mod=100007;
const int maxn=5e6+5;	// 1e7+5??

int notpm[maxn],pm[maxn],pcnt=0,f[maxn];
int n,m,k;
int mu[maxn],pz[maxn],yz;
int hd[11][mod+1],nxt[mod+1],cur[maxn],clo=0,curres[maxn];

int gcd(int a,int b){return !b?a:gcd(b,a%b);}

void xxs(){
	notpm[1]=mu[1]=1;
	for(int i=2;i<=maxn-5;i++){
		if(!notpm[i]){
			pm[++pcnt]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=pcnt && pm[j]*i<=maxn-5;j++){
			notpm[i*pm[j]]=1;
			if(i%pm[j]==0){
				mu[i*pm[j]]=0;
				break;
			}
			mu[i*pm[j]]=-mu[i];
		}
	}
	for(int i=1;i<=maxn-5;i++)mu[i]+=mu[i-1];
	for(int i=1;i<=maxn-5;i++)f[i]=f[i-1]+(gcd(i,k)==1);
}

int gf(int x){	// f[x]
	return x/k*f[k]+f[x%k];
}

map<int,int>exi;
int djs(int x){	// 杜教筛求莫比乌斯函数前缀和 
	if(x <= maxn-5)return mu[x];
	if(exi[x])return exi[x];
	
	int ans=1;
	for(int i=2,j=0;i<=x;i=j+1){
		j=x/(x/i);
		ans-=(j-i+1)*djs(x/i);
	}
	exi[x]=ans;
	return exi[x];
}

map<int,int>hs;
int g(int y,int x){	// g(y,x)=sum_{i=1}^y [(i,x')==1]mu(i)  x 表示x'剩下的因子数而不是原数x'
	if(!x){
		int tmp=djs(y);	 // 边界情况 mu的前缀和，杜教筛优化 
//		printf("hello %d\n",tmp);
		return tmp;
	}
	if(y<=1)return y;
	
	// 对于每个 x y 分开讨论，一开始写的map一直不对 
	for(int i=hd[x][y%mod];i;i=nxt[i])if(cur[i]==y)return curres[i];	// 类似前向星遍历 
	int now=++clo;cur[clo]=y;nxt[clo]=hd[x][y%mod];hd[x][y%mod]=clo;	// 存边 
	int tmp=g(y,x-1)+g(y/pz[x],x);
//	printf("%d\n",tmp);
	return curres[now]=tmp;	// 记忆化 
}

signed main(){
	scanf("%lld%lld%lld",&n,&m,&k);
	xxs();
	for(int i=1;pm[i]<=k;i++)if(k%pm[i]==0)pz[++yz]=pm[i];	// 求出 k 的所有质因数 
	int ans=0;
	for(int i=1,j=0;i<=n && i<=m;i=j+1){
		j=min(n/(n/i),m/(m/i));
		ans+=(g(j,yz)-g(i-1,yz))*(n/i)*gf(m/i);	// 根据推的式子分块 
//		printf("%d %d %d %d\n",ans,g(yz,j),g(yz,i-1),gf(m/i));
	}
	printf("%lld\n",ans);

	return 0;
}