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-8
using 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;
}