HDU MophuesMophues

                                                Mophues

题目：

   求1 <=x <= a,1 <= y <= b。求gcd(x,y)的值得素因子个数不超过P的对数。

算法:

   前面的做了这种类型的好几道了，都要吐了。T_T

  如果，此题没有要求素因子不超过P则就是前面博文中给出过的，Mobius + 分块思想。

  而这题多了一个限制P则在求解前缀的时候加上这个条件就好了。

  为什么想到用Mobius呢？有人给出了结论：

[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) 根据容斥原理： A(a,b,d) = (1)*B(a,b,d) + (-1)*B(a,b,2*d)+ (-1) *B(a,b,3*d) +(0)*B(a,b,4*d)+(-1)*B(a,b,5*d)+(1)*B(a,b,6*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的影响因子大小。先计算完单个，再计算前缀和（接下来有用）。 接着，我们发现对于某个d,会满足B(a,b,d) = (B,a,b,d+x),而且 这个 x = min(a/(a/d),b/(b/d)) ,那么整个式子的计算会呈现块状，因此计算这个区间的时候可以用前缀和。 

#include <iostream>
#include <algorithm>
#include <queue>
#include <stack>
#include <vector>
#include <map>
#include <cstdio>
#include <cstring>
using namespace std;

typedef long long LL;
typedef pair<int,int> P;

inline int read(){
    int x = 0,f = 1; char ch = getchar();
    while(ch < '0'||ch > '9'){if(ch == '-')f=-1;ch = getchar();}
    while(ch >= '0'&&ch <= '9'){x = x * 10 + ch -'0';ch = getchar();}
    return x*f;
}

//

const int MAXN = 500000 + 10;
LL primes,mu[MAXN+10],phi[MAXN+10],prime[MAXN+10];
LL pre[MAXN+10][22],pricnt[MAXN+10];
LL n,m,p;
bool com[MAXN+1];

void getPrim(){
    mu[1] = 1; primes = 0;
    
    memset(pricnt,0,sizeof(pricnt));
    
    for(int i = 2;i <= MAXN;++i){
        if(!com[i]){
            mu[i] = -1;
            prime[primes++] = i;
            pricnt[i] = 1;
        }
        for(int j = 0;j < primes&&i*prime[j] <= MAXN;++j){
            com[i*prime[j]] = 1;
            pricnt[i*prime[j]] = pricnt[i] + 1;
            if(i%prime[j]) mu[i*prime[j]] = -mu[i];
            else { mu[i*prime[j]] = 0; break; }
        }
    }
}

void getMobi(){
    for(int i = 1;i <= MAXN;++i){
        for(int j = i;j <= MAXN;j += i){
            pre[j][pricnt[i]] += mu[j/i];
        }
    }

    for(int i = 1;i <= MAXN;++i){
        for(int j = 0;j < 20;++j){
            pre[i][j] += pre[i-1][j];
        }
    }

    for(int i = 0;i <= MAXN;++i){
        for(int j = 1;j < 20;++j){
            pre[i][j] += pre[i][j-1];
        }
    }
}

int main()
{
    getPrim();
    getMobi();
    int T;
    T = read();
    while(T--){
        n = read(); m = read(); p = read();
        if(p >= 19){
            printf("%I64d\n",(LL)n*m);
            continue;
        }
        if(n > m) swap(n,m);

        LL ans = 0;

        for(int i = 1,last;i <= n;i = last + 1){
           last = min(n/(n/i),m/(m/i));
           ans += LL(pre[last][p] - pre[i-1][p]) * (n/i) * (m/i);
        }

        printf("%I64d\n",ans);
    }
    return 0;
}

不知道是不是分块处理的拙了，居然要用那多么的时间！