【数论-莫比乌斯入门】hdu 1695 GCD

莫比乌斯反演：

                                                                                      

莫比乌斯反演变形：

                               假设：n<m

                               :gcd(i,j)=d的倍数的数对的个数；

                             ：gcd(i,j)=d的数对的个数；

参考blog：https://blog.sengxian.com/algorithms/mobius-inversion-formula

kuangbin的模板代码不错！

*****************************************************************************************************************************************************************

hdu 1695 GCD

题意：给出两个区间x:(1,n)，y:(1,m);求gcd(x,y)=d的数对个数；比如d=1、(4，5)(5，4)算一对；

思路：

根据数学知识易知gcd(x,y)=d的倍数的数对个数为：

利用莫比乌斯反演，使用变形公式:

则，

那么计算f(d)的复杂度为O(n)，优化：可以对(n/dx)*(m/dx)进行分块：

找到某一段，这一段的(n/dx)都相等，(m/dx)都相等，求取μ的前缀和，批量计算这一段的答案。

复杂度变为O(sqrt(n)+sqrt(m))；

分块举例：例如，n=32,m=40,d=2； x：1~min(n,m)/d

分块代码：

last=min(n/(n/i),m/(m/i))

代码：

#include<cstdio>#include<cstring>#include<algorithm>using namespace std;typedef long long ll;const int N=100100;bool check[N+10];int prime[N+10];int mu[N+10];void Moblus(){    memset(check,false,sizeof(check));    mu[1]=1;    int tot=0;    for(int i=2;i<=N;i++)    {        if(!check[i])        {            prime[tot++]=i;            mu[i]=-1;        }        for(int j=0;j<tot;j++)        {            if(i*prime[j]>N) break;            check[i*prime[j]]=true;            if(i%prime[j]==0)            {                mu[i*prime[j]]=0;                break;            }            else            {                mu[i*prime[j]]=-mu[i];            }        }    }}int sum[N+10];int main(){    int t;    Moblus();    sum[0]=0;    for(int i=1;i<=N;i++)        sum[i]=sum[i-1]+mu[i];//    for(int i=1;i<=5;i++)//        printf("%d\n",sum[i]);    scanf("%d",&t);    int co=0;    while(t--)    {        int a,n,c,m,k;        scanf("%d%d%d%d%d",&a,&n,&c,&m,&k);        if(n>m) swap(n,m);        if(k==0)        {            printf("Case %d: 0\n",++co);            continue;        }        else            n/=k,m/=k;        ll ans1=0,ans2=0,ans=0;        for(int i=1,la=1;i<=n;i=la+1)        {            la=n/(n/i);                   //分块；            ans1+=(ll)(sum[la]-sum[i-1])*(n/i)*(n/i);        }        ans1/=2;//        printf("%d\n",ans1);        for(int i=1,la=1;i<=n;i=la+1)        {            la=min(n/(n/i),m/(m/i));      //分块；            ans2+=(ll)(sum[la]-sum[i-1])*(n/i)*(m/i);        }//        printf("%d\n",ans2);        ans=ans2-ans1;        printf("Case %d: %I64d\n",++co,ans);    }    return 0;}