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

Posted 2021-01-05 sciorz

---

分析

其实要去求的就是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_{d|n}mu(d)sum {i=1}^m mu(i)[d|gcd(i,n)] =mu(n)sum{d|n}mu(d)sum {i=1}^{ lfloor frac{m}{d} floor } mu(id) =mu(n)sum{d|n}mu(d)ans(lfloor frac{m}{d} floor,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;
}