hdu4746莫比乌斯反演+分块

http://blog.csdn.net/mowayao/article/details/38875021

题意： 5000组样例。 问你[1,n] 和 [1,m]中有多少对数的GCD的素因子个数小于p。

思路：

首先考虑一个相对简单的版本： [1,a] 和 [1,b] 有多少对的数  满足GCD <= d

首先定义两个函数：A(a,b,d) 表示 GCD(a,b) = d的对数，B（a,b,d）表示GCD（a,b） 是d的倍数的对数 易得 B(a,b,d) = (a/d)*(b/d) 根据容斥原理：

B(a,b,i) 前面的系数正是莫比乌斯函数的值。

那么公式可以写成：

A(a,b,1) =  u(1)*B(a,b,1) + u(2)*B(a,b,2) + u(3) *B(a,b,3) + u(4)*B(a,b,4) + u(5)*B(a,b,5) + u(6)*B(a,b,6)...........

A(a,b,2) =  u(1)*B(a,b,2) + u(2)*B(a,b,4) + u(3) *B(a,b,6) + u(4)*B(a,b,8) + u(5)*B(a,b,10) +u(6)*B(a,b,12).........

A(a,b,3) =  u(1)*B(a,b,3) + u(2)*B(a,b,6) + u(3) *B(a,b,9) + u(4)*B(a,b,12) + u(5)*B(a,b,15) +u(6)*B(a,b,18).......

A(a,b,4) =  u(1)*B(a,b,4) + u(2)*B(a,b,8) + u(3) *B(a,b,12) + u(4)*B(a,b,16) + u(5)*B(a,b,20) +u(6)*B(a,b,24).....

答案就是

A(a,b,1)+A(a,b,2)+A(a,b,3)+......A(a,b,d) =   u(1)*B(a,b,1)+(u(1)+u(2))*B(a,b,2) + ....... (u(1)+u(2)+u(3)+u(6))*B(a,b,6)........

可见A(a,b,d) 前的系数为  sigma(u(i)) (i为d的约数) =  C（a,b,d）

 

然后，这一题还有一个限制条件，就是要使素因子的个数小于等于p,那么我们定义这个函数D(a,b,d,p) 表示B(a,b,d) 前的系数，那么我们只要从C(a,b,d)中选出一些满足条件的系数即可。 用一个数组F[d][cnt] (cnt为素因子个数)记录，数组表示的是d的因子的素因子个数为cnt的影响因子大小。先计算完单个，再计算前缀和（接下来有用）。(我们可以知道只有那些素因数 小于等于p的才会用到B(a,b,d),因此只要在B(a,b,d)的位置留下他相应的值就ok了)接着，我们发现对于某个d,会满足B(a,b,d) = (B,a,b,d+x),而且  这个 x = min(a/(a/d),b/(b/d)) ,那么整个式子的计算会呈现块状，因此计算这个区间的时候可以用前缀和。

#include <iostream>
#include <algorithm>
#include <string.h>
#include <cstdio>
#include <vector>
using namespace std;
const int maxn=500099;
int mu[maxn];
int prime[maxn],primenum[maxn];
bool isprime[maxn];
int F[maxn][20];
void getmu()
{
     mu[1]=1;
     memset(isprime,true,sizeof(isprime));
     isprime[0]=isprime[1]=false;
     int cnt=0;
     primenum[1]=0;
     for(int i=2; i<maxn; i++)
        {
             if(isprime[i])
                {
                    prime[cnt++]=i;
                    mu[i]=-1;
                    primenum[i]=1;
                }
                for(int j=0; j<cnt && (prime[j]*i)<maxn; j++)
                {
                    primenum[i*prime[j]]=primenum[i]+1;
                    mu[i*prime[j]]=-mu[i];
                    isprime[i*prime[j]]=false;
                    if( (i%prime[j]) == 0)
                        {
                             mu[i*prime[j]]=0;break;
                        }
                }
        }
}
void getmF()
{
     memset(F,0,sizeof(F));
     for(int i=1; i<maxn; i++){
        for(int j=i; j<maxn; j+=i)
        {
            F[j][primenum[i]]+=mu[j/i];
        }
     }
     for(int i=1; i<maxn; i++)
        for(int j=0; j<=19;j++)
         F[i][j]+=F[i-1][j];
     for(int i=1; i<maxn; i++)
        for(int j=1; j<=19; j++)
         F[i][j]+=F[i][j-1];
}
long long solve(int n,int m, int p)
{
    long long ans=0;
    int ed=0;
    for(int i=1; i<=n; i++)
        {
            ed=min( n/(n/i),m/(m/i));
            ans+=1LL * ( F[ed][p]-F[i-1][p] )*(n/i)*(m/i);
            i=ed;
        }
        return ans;
}
int main()
{
   getmu();
   getmF();
    int cas;
    scanf("%d",&cas);
    for(int cc=1; cc<=cas; cc++)
        {
               int n,m,p;
               scanf("%d%d%d",&n,&m,&p);
               if(p>19)
                {
                    printf("%I64d\n",1LL*n*m); continue;
                }
               if(n>m)swap(n,m);
               long long ans=solve(n,m,p);
              printf("%I64d\n",ans);
        }

    return 0;
}