ACM-ICPC 2018 徐州赛区网络预赛 D.easy math

---

分析

其实要去求的就是1-m之间与n互质的数的莫比乌斯函数之和。
这样我们可以枚举n的因数d，然后再容斥地加上（或减去）1-m之间与n的gcd为d的倍数的数的莫比乌斯函数之和。
$ ans(m,n)=\sum _{i=1}^m \mu(in) =\mu(n)\sum {i=1}^m \mu(i)[gcd(i,n)==1] =\mu(n)\sum\mu(d)\sum {i=1}^m \mu(i)[d|gcd(i,n)] =\mu(n)\sum\mu(d)\sum {i=1}^{ \lfloor \frac{m}{d} \rfloor } \mu(id) =\mu(n)\sum\mu(d)ans(\lfloor \frac{m}{d} \rfloor,d) $
参考：https://www.zybuluo.com/yang12138/note/1277248

代码

#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N=1e6+200;
bool isprime[N];
int prime[N],pz,mob[N];
int mob_sum[N];
unordered_map<int,int> Mob_sum;
void init() {
	memset(isprime,true,sizeof isprime);
	isprime[0]=isprime[1]=false;
	mob[1]=mob_sum[1]=1;
	for(int i = 2; i < N; ++i) {
		if(isprime[i]) prime[pz++]=i,mob[i]=-1;
		mob_sum[i]=mob_sum[i-1]+mob[i];
		for(int j = 0; j < pz&&(ll)i*prime[j]<N; ++j) {
			isprime[i*prime[j]]=false;
			mob[i*prime[j]]=-mob[i];
			if(i%prime[j]==0) {
				mob[i*prime[j]]=0;
				break;
			}
		}
	}
}
ll calc(ll n) {
	if(n<N) return mob_sum[n];
	if(Mob_sum.count(n)) return Mob_sum[n];
	ll ans=1;
	for(int i = 2; i <= n; ++i) {
		ll t=n/i,nxt=n/t;
		ans-=(nxt-i+1)*calc(t);
		i=nxt;
	}
	return Mob_sum[n]=ans;
}
ll m,n,a[20],tot;
ll mu(ll x) {
	int ans=1;
	for(int i = 0; i < tot; ++i) {
		if(x%a[i]==0) ans=-ans;
	}
	return ans;
}
map<pair<ll,ll>,ll> mp;
ll solve(ll m, ll n) {
	if(m==0) return 0;
	if(mp.count({m,n})) {
		return mp[{m,n}];
	}
	if(m==1) return mp[{m,n}]=mu(n);
	if(n==1) return mp[{m,n}]=calc(m);
	ll ans=0;
	for(int i = 1; (ll)i*i <= n; ++i) {
		if(n%i) continue;
		ans+=mu(i)*solve(m/i,i);
		ll d=n/i;
		if(d!=i) {
			ans+=mu(d)*solve(m/d,d);
		}
	}
	return mp[{m,n}]=ans*mu(n);
}
int main() {
	init();
	cin>>m>>n;
	ll tmp=n;
	for(int i = 0; (ll)prime[i]*prime[i]<=tmp; ++i) {
		int p=prime[i];
		if(tmp%p==0) {
			tmp/=p;
			a[tot++]=p;
			if(tmp%p==0) {
				puts("0");
				return 0;
			}
		}
	}
	if(tmp>1) a[tot++]=tmp;
	cout<<solve(m,n)<<endl;
	return 0;
}