hdu 4746 Mophues 莫比乌斯 分块优化

Description

As we know, any positive integer C ( C >= 2 ) can be written as the multiply of some prime numbers: 
    C = p1×p2× p3× ... × pk 
which p1, p2 ... pk are all prime numbers.For example, if C = 24, then: 
    24 = 2 × 2 × 2 × 3 
    here, p1 = p2 = p3 = 2, p4 = 3, k = 4 

Given two integers P and C. if k<=P( k is the number of C's prime factors), we call C a lucky number of P. 

Now, XXX needs to count the number of pairs (a, b), which 1<=a<=n , 1<=b<=m, and gcd(a,b) is a lucky number of a given P ( "gcd" means "greatest common divisor"). 

Please note that we define 1 as lucky number of any non-negative integers because 1 has no prime factor.

 

Input

The first line of input is an integer Q meaning that there are Q test cases. 
Then Q lines follow, each line is a test case and each test case contains three non-negative numbers: n, m and P (n, m, P <= 5×10 5. Q <=5000).

 

Output

For each test case, print the number of pairs (a, b), which 1<=a<=n , 1<=b<=m, and gcd(a,b) is a lucky number of P.

 

Sample Input

210 10 010 10 1 

 

Sample Output

6393 

然后由莫比乌斯反演，我们知道f(d)=g(d)*u(1)+g(2*d)*u(2)+g(3*d)*u(3)+....

考虑结果ans，ans为所有f(d)且F(d)<=P的和，

我们枚举1<=i<=n，如果i是某个d的倍数且F(d)<=P，那么ans+=g(i)*u(i/d)=[n/i]*[m/i]*u(i/d)。那么这个怎么计算能更快一点？

我们设G(i)为容斥因子：G(i)=sum{u(i/d) | F(d)<=P} 这个值可以nlogn预处理出来，然后我们只需要ans+=G(i)*[n/i]*[m/i]即可

这样的话总的复杂度为O(n*q)还是会T的样子

然后我们注意到[n/i]*[m/i]在一定的范围内是不变的，这个范围是[i,min(n/(n/i),m/(m/i)]，这样我们可以预处理出G(i)的前缀和，然后加快运算（复杂度网上说是sqrt(n)的。。。）

这样总的复杂度是O(q*sqrt(n)+nlog(n))大概这样.

#include<iostream>#include<cstdio>#include<cmath>#include<algorithm>#include<cstring>#include<vector>#include<set>#include<map>#include<queue>#include<bitset>#define fi first#define se second#define pii pair<int,int>#define ll long long#define inf 1<<30#define eps 1e-8using namespace std;const int maxn=500010;int num[maxn];int cnt[19][maxn];int n,m,p;int mu[maxn],prime[maxn],tot;bool check[maxn];void getmu(){    int N=500000;    memset(check,0,sizeof(check));    mu[1]=1,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];        }    }    for(int i=1;i<=N;i++) {        int u=i;        int v=0;        for(int j=0;j<tot;j++) {            if(prime[j]*prime[j]>u)                break;            while(u%prime[j]==0) {                u/=prime[j];                v++;            }        }        if(u>1) v++;        num[i]=v;    }}void init(){    memset(cnt,0,sizeof(num));    int N=500000;    for(int i=1;i<=N;i++) {        for(int j=i;j<=N;j+=i) {            cnt[num[i]][j]+=mu[j/i];        }    }    for(int i=1;i<=N;i++) {        for(int j=1;j<=18;j++) {            cnt[j][i]+=cnt[j-1][i];        }    }    for(int i=0;i<=18;i++) {        for(int j=1;j<=N;j++) {            cnt[i][j]+=cnt[i][j-1];        }    }}int main(){    getmu();    init();    int t;    scanf("%d",&t);    while(t--) {        scanf("%d%d%d",&n,&m,&p);        if(p>18) {            printf("%I64d\n",(ll)n*m);            continue;        }        ll ans=0;        if(n>m) swap(n,m);        for(int i=1;i<=n;) {            int j=min(n/(n/i),m/(m/i));            //cout<<i<<" "<<j<<endl;            ans+=(ll)(cnt[p][j]-cnt[p][i-1])*(n/i)*(m/i);            i=j+1;        }        printf("%I64d\n",ans);    }    return 0;}