1. 程式人生 > >2018徐州ICPC網路賽D Easy Math 莫比烏斯反演/杜教篩

2018徐州ICPC網路賽D Easy Math 莫比烏斯反演/杜教篩

遞推題解

#include<stdio.h>
#include<string.h>
#include<map>
#include<algorithm>
using namespace std;
const int MAX=1e7+5;
bool vis[MAX];
int prim[MAX],mu[MAX],fac[MAX],tot,pcnt;
map<long long,long long>dp;
void moblus()
{
    mu[1]=1;tot=0;
    for(int i=2;i<MAX;i++)
    {
        if(!vis[i]) {prim[tot++]=i;mu[i]=-1;}
        for(int j=0;j<tot&&i*prim[j]<MAX;j++)
        {
            vis[i*prim[j]]=1;
            if(i%prim[j]) mu[i*prim[j]]=-mu[i];
            else {mu[i*prim[j]]=0;break;}
        }
    }
    for(int i=2;i<MAX;++i) mu[i]+=mu[i-1];
}
long long M(long long x)
{
    if(x<MAX) return mu[x];
    if(dp[x]) return dp[x];
    long long sum=1;
    for(long long l=2,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        sum-=M(x/l)*(r-l+1);
    }
    return dp[x]=sum;
}
long long query(long long m,long long n)
{
    if(m==0) return 0;
    if(n==1) return M(m);
    for(int i=0;i<pcnt;++i)
        if(n%fac[i]==0)
            return query(m/fac[i],n)-query(m,n/fac[i]);
    return query(m/n,n)-query(m,1);
}
int main()
{
    long long n,m;
    moblus();
    while(~scanf("%lld%lld",&m,&n))
    {
        long long x=n;bool f=0;pcnt=0;
        for(long long i=2;i*i<=x;++i)
        {
            if(x%(i*i)==0) {f=1;break;}
            if(x%i==0) {fac[pcnt++]=i;x/=i;}
        }
        if(f) {printf("0\n");continue;}
        if(x>1) fac[pcnt++]=x;
        dp.clear();
        printf("%lld\n",query(m,n));
    }
    return 0;
}