SPOJ PGCD（mobius反演+分块+线性筛）

这道题T了，看了大神代码才知道我把这道题想得太简单了，思路基本上是下面这个博客的。 神牛博客 首先

第一步，

求的是1< =i< =m,1<= j< =n,gcd(i,j)为素数的个数 枚举小于m且小于n的素数p。 然后转化为求1< =i< =m,1<= j< =n，gcd（i，j)为p的个数， （这个就是hdu1695所求的， 解法详见http://blog.csdn.net/abc13068938939/article/details/52198163）

for(int i=0;i<tot&&prime[i]<=a;i++) { int ta=a/prime[i],tb=b/prime[i];//gcd(i,j)=k ==> gcd(i/k,j/k)=1 for(int j=1;j<=ta;j++) res+=mu[j]*(ta/j)*(tb/j); }

当然，T。。。

第二步

回想形式 sigma(F(d) * miu(d / p)) p为素数 =F(d) * sigma(miu(d / p)) p为素数 那么其中的F(d)只和d有关。。。可以枚举d，那么就要求出P[d] = sigma(miu(d / p))了 p是d的某一个素因子。。。 如果d = p1 ^ k1 * p2 ^ k2 * …… pm ^ km。 if (k1,k2,…km都等于1)，那么d/pi就都是m-1个不同的素数相同，并且这样的d/pi总共有m个，所以P[d]=(((m-1)&1)?-1:1)*m; else if(k1,k2,…km只有一个kii等于2，其余都等于1)，那么d/pi(除了i=ii之外)就都有重复的素数，miu[d/pi]=0，所以P[d]=miu[d/pii]= (m&1)?-1:1; else P[d]=0; 那么问题来了怎么知道d的k1,k2,k3,…km，因数分解？太麻烦复杂度太高了。 我们用totnum[i]表示i的全部素数因子个数（含重复），diffnum[i]表示i的全部不同的素数因子个数。 那么（k1,k2,…km都等于1）就相当于totnum[i]==diffnum[i]， (k1,k2,…km只有一个kii等于2，其余都等于1)就相当于totnum[i]==diffnum[i]+1。 初始化i是质数 totnum[i]=diffnum[i]=1; 这里我们想到线性筛在求质数时可以将每个合数的最小质因子求出来，记为minfac[i]， 这样对于每一个合数i totnum[i]=totnum[i/minfac[i]]+1; diffnum[i]=diffnum[i/minfac[i]]+((i/minfac[i])%minfac[i]!=0);

for(int i=2;i<=a&&i<=b;i++) res+=(LL)p[i]*(a/i)*(b/i);

然而还是T T_T||

境界三： 分块 又是经典的分块，考虑F(d) = (n / d) * (m / d)，这个东西相邻的部分值是一样的。 对于这种分块不熟的可以做CQOI2007 余数之和sum。

下面代码的p[i]是上面所讲p[i]的前缀和。 AC代码

#include <iostream> #include <stdio.h> #include <algorithm> #include <stdlib.h> #include <stack> #include <vector> #include <string.h> #include <queue> #define msc(X) memset(X,-1,sizeof(X)) #define ms(X) memset(X,0,sizeof(X)) typedef long long LL; using namespace std; const int MAXN=10000010; bool check[MAXN+10]; int prime[MAXN/3],tot,minumfac[MAXN+10]; int totnum[MAXN+10],diffnum[MAXN+10],p[MAXN+10]; void Moblus(void) { ms(check); tot=0; check[0]=check[1]=true; //mu[1]=1; for(int i=2;i<=MAXN ;i++) { if(!check[i]) { prime[tot++]=i;//mu[i]=-1; } for(int j=0;j<tot;j++) { if(i*prime[j]>MAXN) break; check[i*prime[j]]=true; minumfac[i*prime[j]]=prime[j]; if(i%prime[j]==0){ //mu[i*prime[j]]=0; break; } //else mu[i*prime[j]]=-mu[i]; } } p[1]=0; for(int i=2;i<=MAXN;i++) { if(!check[i]) totnum[i]=diffnum[i]=1; else { totnum[i]=totnum[i/minumfac[i]]+1; diffnum[i]=diffnum[i/minumfac[i]]+ ((i/minumfac[i])%minumfac[i]!=0); } if(totnum[i]==diffnum[i]) p[i]=(((totnum[i]-1)&1)?-1:1)*totnum[i]; else if(totnum[i]-diffnum[i]==1) p[i]=(diffnum[i]&1)?-1:1; else p[i]=0; p[i]+=p[i-1]; } } int main(int argc, char const *argv[]) { int t; cin>>t; Moblus(); while(t--){ int a,b; scanf("%d %d",&a,&b); //if(a>b) {int tmp=a;a=b;b=tmp;}//保证a<=b LL res=0ll; for(int i=2,r;i<=a&&i<=b;i=r+1)//已保证ta <= tb { int d1=a/i,d2=b/i; r=min(a/d1,b/d2); res+=(LL)(p[r]-p[i-1])*d1*d2; } printf("%lld\n",res ); } return 0; }