NOIP模拟赛20191112 T3 c【推式子(莫比乌斯反演)】

题目描述：

求：
$$\sum _{d|n,d|s}\sum _{g|\frac{n}{d}}\mu (g)\ast \lfloor \frac{m}{dg}\rfloor$$
其中$n,m,s$满足$\{\begin{array}{l}n=n_1\cdot n_2\cdot n_3\\ m=m_1\cdot m_2\cdot m_3\\ s=s_1\cdot s_2\cdot s_3\end{array}$
多组数据，数据组数$\le 1000$，$n_1,n_2,n_3,m_1,m_2,m_3,s_1,s_2,s_3\le 10^6$

题目分析：

看到这个东西并且联想到NOIP的标题以及yali标准的模拟赛样式就明白这一定是把莫比乌斯反演的式子推到某一步进行了奇奇怪怪的变换然后拿来恶心人防AK。

于是我们开始按照题解推式子：
$$\sum _{d|n,d|s}\sum _{g|\frac{n}{d}}\mu (g)\ast \lfloor \frac{m}{dg}\rfloor$$
首先把$\mu$函数的系数想办法搞掉：
$$\sum _{d|n,d|s}\sum _{g|\frac{n}{d}}\sum _{i=1}^{\lfloor \frac{m}{dg}\rfloor }\mu (g)$$
联想到${\sum }_{g|k}\mu (g)=[k=1]$，把第二个∑放到里层，原来的后两式相当于枚举$\lfloor \frac{m}{d}\rfloor$中g的倍数，转换一下可变为：
$$\sum _{d|n,d|s}\sum _{i=1}^{\lfloor \frac{m}{d}\rfloor }\sum _{g|\frac{n}{d},g|i}\mu (g)$$
整除两个数就相当于整除它们的gcd，于是：
$$\sum _{d|n,d|s}\sum _{i=1}^{\lfloor \frac{m}{d}\rfloor }[gcd(i,\frac{n}{d})=1]$$
先不看$d|s$的条件，观察式子，发现其实$d$就是在枚举$[1,m]$中的数与$n$的$gcd$，如果没有$d|s$的条件，则式子就 $=m$
现在加上$d|s$的条件，那么只有与$n$的$gcd$为$s$的约数的数会被统计到，即：
$$\sum _{i=1}^m[gcd(i,n)|s]$$
到这里可以$O(m\log m)$枚举计算，结合部分分拿到70分。

Code：

#include<bits/stdc++.h>
#define LL long long
#define x first
#define y second
using namespace std;
int T,num,bit[1<<15];
LL n,m,s,ans,f[1<<15];
vector<pair<int,int> >a,b;
void Divide(int x,vector<pair<int,int> >&q){
	for(int i=q.size()-1;i>=0;i--) while(x%q[i].x==0) q[i].y++,x/=q[i].x;
	for(int i=2;i*i<=x;i++) if(x%i==0){
		q.push_back(make_pair(i,0));
		while(x%i==0) q.back().y++,x/=i;
	}
	if(x>1) q.push_back(make_pair(x,1));
}
int main()
{
	//freopen("c.in","r",stdin);
	//freopen("c.out","w",stdout);
	scanf("%d",&T);
	while(T--){
		n=m=s=1,ans=num=0,a.clear(),b.clear();
		for(int i=0,x;i<3;i++) scanf("%d",&x),Divide(x,a);
		for(int i=0,x;i<3;i++) scanf("%d",&x),m*=x;
		for(int i=0,x;i<3;i++) scanf("%d",&x),Divide(x,b);
		sort(a.begin(),a.end()),sort(b.begin(),b.end());
		for(int i=0,j=0;i<a.size();i++){
			while(j<b.size()&&b[j].x<a[i].x) j++;
			if(j==b.size()||b[j].x>a[i].x||b[j].y<a[i].y){
				f[1<<num]=a[i].x;int ed=j==b.size()||b[j].x>a[i].x?0:b[j].y;
				for(int k=1;k<=ed;k++) f[1<<num]*=a[i].x;
				num++;
			}
		}
		ans=m,f[0]=1;
		for(int s=1;s<1<<num;s++){
			bit[s]=bit[s>>1]+(s&1);
			f[s]=f[s^(s&-s)]*f[s&-s];
			ans+=(bit[s]&1?-1:1)*m/f[s];
		}
		printf("%lld\n",ans);
	}
}