题解 洛谷 P1587 【[NOI2016]循环之美】

从纯循环与不能重复两个条件可以发现，题目本质上要我们求：

$$\sum_{x=1}^n \sum_{y=1}^m [x \perp y][y \perp k]$$

即是：

$$\sum_{x=1}^n \sum_{y=1}^m [y \perp k][gcd(x,y)==1]$$

将 $[gcd(x,y)==1]$ 反演得：

$$\sum_{x=1}^n \sum_{y=1}^m [y \perp k] \sum_{d|x,d|y} \mu(d)$$

把 $d$ 提到前面：

$$\sum_{d=1}^{min(n,m)} \mu(d)[d \perp k]\lfloor \frac{n}{d} \rfloor \sum_{y=1}^{\lfloor \frac{m}{d} \rfloor} [y \perp k]$$

设：

$$f(n,k) = \sum_{i=1}^n [i \perp k]$$

$$g(n,k) = \sum_{i=1}^n \mu(i)[i \perp k]$$

易知:

$$f(n,1) = n$$

$$g(n,1) = \sum_{i=1}^n \mu(i)$$

前者可以直接算，后者可以杜教筛。

考虑到：$k$ 中有无平方因子与互质的条件无关，可假设 $k$ 没有平方因子。

设 $k$ 某个素因子为 $p$，考虑 $f(n,k)$ 较于 $f(n,k/p)$ 的变化。

因为多了一个素因子，故需要减去含有因子 $p$ 的一部分，得：

$$f(n,k) = \sum_{i=1}^n [i \perp \frac{k}{p}] - \sum_{i=1}^n [i \perp \frac{k}{p}][p|i]$$

$$ = f(n,k/p) - \sum_{i=1}^{\lfloor \frac{n}{p} \rfloor} [ip \perp \frac{k}{p}]$$

$$ = f(n,k/p) - \sum_{i=1}^{\lfloor \frac{n}{p} \rfloor} [i \perp \frac{k}{p}]$$

$$ = f(n,k/p) - f(\lfloor n/p \rfloor,k/p)$$

同理可知 $g(n,k)$ 与 $g(n,k/p)$ 的关系：

$$g(n,k) = \sum_{i=1}^n \mu(i)[i \perp \frac{k}{p}] - \sum_{i=1}^n \mu(i)[i \perp \frac{k}{p}][p|i]$$

$$ = g(n,k/p) - \sum_{i=1}^{\lfloor \frac{n}{p} \rfloor} \mu(ip)[ip \perp \frac{k}{p}]$$

$$ = g(n,k/p) - \sum_{i=1}^{\lfloor \frac{n}{p} \rfloor} \mu(i)\mu(p)[i \perp p][i \perp \frac{k}{p}]$$

$$ = g(n,k/p) - \mu(p)\sum_{i=1}^{\lfloor \frac{n}{p} \rfloor} \mu(i)[i \perp k]$$

$$ = g(n,k/p) + g(\lfloor n/p \rfloor,k)$$

式中用到结论：

$$\mu(ab) = \mu(a)\mu(b)[a \perp b]$$

所以，对于每个 $f(n,k)$、$g(n,k)$ 可以递推得出，然后搞个数论分块，就可以了。

代码如下：

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;

const int maxnn=1000000010;
const int maxn=1000010;

map<pair<int,int>,int> f,g;
int n,m,k;
ll ans;
int tot,prime[maxn/10],mu[maxn],hash[maxn];
int smu[2*maxnn/maxn];

void Euler()
{
	mu[1]=1;
	for(int i=2;i<=maxn;i++)
	{
		if(!hash[i])
		{
			mu[i]=-1;
			prime[++tot]=i;
		} 
		for(int j=1;j<=tot&&i*prime[j]<=maxn;j++)
		{
			int c=i*prime[j];
			hash[c]=1;
			if(i%prime[j]) mu[c]=-mu[i];
			else break;
		}
	}
	for(int i=1;i<=maxn;i++) mu[i]+=mu[i-1];
}

int calc_mu(int d)
{
	if(d<=maxn) return mu[d];
	else
	{
		if(smu[n/d+m/d]) return smu[n/d+m/d];//防止 n，m 断点不同冲突
		else 
		{
			int z=1;
			for(int i=2;i<=d;i++)
			{
				int t=d/i,j=d/t;
				z-=(j-i+1)*calc_mu(t);
				i=j;
			}
			smu[n/d+m/d]=z;
			return z;
		}
	}
}

int calc_f(int q,int r)
{
	if(q==0) return 0;
	if(r==1) return q;
	else
	{
		if(f[make_pair(q,r)]) return f[make_pair(q,r)];
		else
		{
			int &z=f[make_pair(q,r)];
			for(int i=1;prime[i]<=k;i++)
			{
				int c=prime[i];
				if(r%c==0)
				{
					z+=calc_f(q,r/c);
					if((r/c)%c!=0) z-=calc_f(q/c,r/c);//是否为唯一素因子
					return z;
				}
			}
			z=calc_f(q,1)-calc_f(q/r,1);
			return z;
		}
	}
}

int calc_g(int q,int r)
{
	if(q==0) return 0;
	if(r==1) return calc_mu(q);
	else
	{
		if(g[make_pair(q,r)]) return g[make_pair(q,r)];
		else
		{
			int &z=g[make_pair(q,r)];
			for(int i=1;prime[i]<=k;i++)
			{
				int c=prime[i];
				if(r%c==0)
				{
					z=calc_g(q,r/c);
					if((r/c)%c!=0) z+=calc_g(q/c,r);//是否为唯一素因子
					return z;
				}
			}
			z=calc_g(q,1)+calc_g(q/r,r);
			return z;
		}
	}	
}

int main()
{
	scanf("%d %d %d",&n,&m,&k);
	Euler();
	for(int i=1;i<=min(n,m);)
	{
		int e1=n/(n/i),e2=m/(m/i);
		e1=min(e1,e2);
		ans+=(ll)(n/i)*(calc_g(e1,k)-calc_g(i-1,k))*calc_f(m/i,k);
		i=e1+1;
	}
	printf("%lld",ans);
	return 0;
}