[SPOJ pgcd]【容斥原理】【常数】

【地址】http://www.spoj.pl/problems/PGCD/

【题目大意】

10组数据，给定(m,n)，问gcd(i,j)=素数的对数(1<=i<=m,1<=j<=n)   m,n<=10^7

【算法】

不妨设m

先筛出素数和miu函数。

然后暴力容斥。

for (int d=2;d<=m;d++) if (prime[d])

  for (int j=1;j<=m/d;j++)

     ans+=miu(j)*(m/d/j)*(n/d/j)

这样是 m lg m的

miu是容斥系数。

特别地miu[1]=1。

miu[i]=1      i为Square-Free-Number并且i有偶数个素因子

miu[i]=-1     i为Square-Free-Number并且i有奇数个素因子

miu[i]=0      other

这样会TLE。

那么看我们所求。

sigma{    miu[j]*( m/(d*j) )*n/( n/(d*j) )   }    其中d*j<=m && d为素数

尝试枚举d*j的值。令x=d*j，由于d为素数，令x=p1^k1 * p2^k2 * p3^k3 … 其中pi为素数。

那么现在看看会出什么结果。

1、

假如ki全为1。那么从里面拿出一个素数d，会剩下一个Square-Free-Number。

设x有t个素因子。

j含有的素因子个数必然为t-1。所以所有的j的miu值都一样。而取法有t种于是ans+= ( ((t-1)&1)?-1:1 ) * t * (m/x) * (n/x)

2、

假如ki有一个是2，其它全是1，那么只能拿出那个ki=2的，否则剩下的不是SFN，miu值=0，无须计算。

这样设拿走了那个ki=2的一个素因子以后，剩下t个素因子。那么显然 ans+=( (t&1)?-1:1 )*(m/x)*(n/x)

3、

其它的不可能拿走一个素数以后得到SFN了。于是miu值都=0，不必考虑。

 

于是。。。枚举x。然后按照上面的情况分别搞就行了。

实际上都是P[x]*(m/x)*(n/x)的形式。由于上面的分类是针对x这个数值的，所以可以通过筛选法把P[x]预处理出来。

然后for (int x=2;x<=m;x++) ans+=P[x]*(m/x)*(n/x)就行了。

但是这样还不够快。

注意到m/x和n/x都只有sqrt级别个结果。于是随着x的递增，(m/x)*(n/x)只有sqrt(m)+sqrt(n)级别个结果。

通过对P[]建立前缀和数组，可以迅速知道区间的P[]之和。那么对于(m/x)*(n/x)相同的部分，直接跳过一整段，

ans+=(S_P[r]-S_P[l-1])*(m/x)*(n/x)     其中(m/l) * (n/l) == (m/r) * (n/r) 。

 

【时间复杂度】

一开始的筛O(n)的

后面每个询问sqrt(n)的。

 

【其它】

出现了很诡异的事件- -，我用O(n)的筛，45S。

用O(n lg n)的筛，31S。

不明真相中。

 

【CODE】

#include

#include

#include

#include

#include

using namespace std;

const int N=(int)1e7;

int F[N+5];

int cnt_1[N+5];

int cnt_2[N+5];

int S[N+5];

 

void Get_Prime(){

    S[2]=0;

    for (int i=2;i<=N;i++){

        if (cnt_1[i]==0){

            F[i]=i;

            cnt_1[i]=1;

//这个筛是n lg n的- -，但不知道为什么，居然比线性筛快。求指导。。。

/*

            for (int j=i+i,c=2;j<=N;j+=i,c++){

              cnt_1[j]++;

              if (c==i){

                cnt_2[j]++;

                if (cnt_2[j]==1 && j/i/i%i==0) cnt_2[j]++;

                c=0;

              }

            }

*/

        }

 

        for (int j=2,k=i+i;j<=F[i];j++,k+=i){

            if (k>N) break;

            F[k]=j;

            cnt_2[k]=cnt_2[i]+(j==F[i]);  //计算k有多少个i^2形式的因子

            cnt_1[k]=cnt_1[i]+1;    //计算k有多少个因子。

        }

 

        int t;

        int &c1=cnt_1[i];

        int &c2=cnt_2[i];

        if (c2==0)

          if (c1&1)

            t=c1;

          else

            t=-c1;

        else

            if (c2==1) t=( (c1&1)?1:-1 );

            else t=0;

        S[i+1]=S[i]+t;

    }

}

 

int main(){

    Get_Prime();

    int Tc,m,n,d1,d2,Next_1,Next_2,Next;

    scanf("%d",&Tc);

    while (Tc–){

        scanf("%d%d",&m,&n);

        if (m>n) swap(m,n);

        long long ans=0;

        for (int i=2;i<=m;){

            d1=m/i;

            d2=n/i;

            Next_1=m/d1+1;          //Next_1>i && m/Next_1>m/i

            Next_2=n/d2+1;          //they are the same

            Next=Next_1

            ans+=(long long)(S[Next]-S[i])*d1*d2;

            i=Next;

        }

        printf("%lldn",ans);

    }

}