hdu 4746 Mophues 莫比乌斯反演+前缀和优化

Mophues

题意：给出n, m, p，求有多少对a, b满足gcd(a, b)的素因子个数<=p，(其中1<=a<=n, 1<=b<=m)

有Q组数据；(n, m, P <= 5×10⁵. Q <=5000).

参考:ACdreamers

思路:对于hdu1695 GCD来说，由于只需要求gcd = k的个数，所以我们可以按照不优化莫比乌斯公式直接求，这样求解的时间复杂度为O(n);

若这题我们也不优化，答案累加需要双重循环:

rep1(i,1,n){
     if(num[i] <= p){
           for(int j = i;j <= n;j += i) //优化对象
               ans += mu(c[j/i])*(n/j)*(m/j);
    }
}

时间复杂度直接为O(n^2),其中n = min(n,m);num[i]表示i的素因子个数；

优化:对于hdu 1695我们是按照公式2即枚举d,然后得到d的倍数j,乘以mu[]的sigma和；这时发现后面的F[j]一直在变，这是因为我们求解的是特定的f[d].但是在这一题中求的不是对于特定的gcd == d的个数，而是一个区间的个数(并且还是gcd的素因子的个数的区间)。我们就需要看对于每一个F[n]的总系数Σmu[]为多少?(这样就可以线性处理了)如果没有素因子的限制，我们就可以直接在预处理出mu[i]之后，枚举i的倍数j,累加即可；现在有系数了，其实也就是只多了一个维度，我们将j的因子的素因子个数加到对应的位置([j][num[i]]),这样再前缀求和就可以弄出素因子个数在p以内了；

原本以为预处理出了每个j的Σmu后即可不再优化了，即每组数据线性处理，但是数据组数过多，导致还是TLE了。这里还需要用到分块的思想；即再一次计算前缀和mu[1...j][k],因为对于d ε [i,n/(n/i)],有n/d = n/i;这样我们不需要枚举i了，直接整块处理即可；

时间复杂度:预处理是O(nlog(n))(也是空间复杂度)；之后每组数据时间复杂度为O(sqrt(min(n,m))),分块加速；

#include<iostream>
#include<cstdio>
#include<cstring>
#include<string.h>
#include<algorithm>
#include<vector>
#include<cmath>
#include<stdlib.h>
#include<time.h>
#include<stack>
#include<set>
#include<map>
#include<queue>
using namespace std;
#define rep0(i,l,r) for(int i = (l);i < (r);i++)
#define rep1(i,l,r) for(int i = (l);i <= (r);i++)
#define rep_0(i,r,l) for(int i = (r);i > (l);i--)
#define rep_1(i,r,l) for(int i = (r);i >= (l);i--)
#define MS0(a) memset(a,0,sizeof(a))
#define MS1(a) memset(a,-1,sizeof(a))
#define MSi(a) memset(a,0x3f,sizeof(a))
#define inf 0x3f3f3f3f
#define lson l, m, rt << 1
#define rson m+1, r, rt << 1|1
typedef pair<int,int> PII;
#define A first
#define B second
#define MK make_pair
typedef __int64 ll;
template<typename T>
void read1(T &m)
{
    T 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();}
    m = x*f;
}
template<typename T>
void read2(T &a,T &b){read1(a);read1(b);}
template<typename T>
void read3(T &a,T &b,T &c){read1(a);read1(b);read1(c);}
template<typename T>
void out(T a)
{
    if(a>9) out(a/10);
    putchar(a%10+'0');
}
const int N = 5e5+7;
const int M = 19;
int num[N],c[N],mu[N][M];
int cal(int n,int p)
{
    int cnt = 0;
    while(n%p == 0){
        cnt++;
        n /= p;
    }
    return cnt;
}
void getprime()
{
    for(int i = 2;i < N;i++)if(!num[i]){//素数;
        for(int j = i;j < N;j += i){
            int t = cal(j,i);
            num[j] += t;
            if(t > 1) c[j] = -1;
            else if(c[j] >= 0) c[j]++;
        }
    }
}
int mobius(int n)
{
    if(n == -1) return 0;
    if(n & 1) return -1;
    return 1;
}
void init()
{
    getprime();
    rep0(i,1,N){
        for(int j = i; j < N;j += i){
            mu[j][num[i]] += mobius(c[j/i]);// 对于每一个j的因数i,最后累加时都有一个mu[j/i]*F[j];
        }
    }
    rep0(i,1,N){
        rep0(j,1,M)
            mu[i][j] += mu[i][j-1];// 计算在p范围内有多少有效
    }
    rep0(j,0,M){
        rep0(i,1,N)
            mu[i][j] += mu[i-1][j];
    }
}
int main()
{
    init();
    int n,m,p,i,j,T;
    read1(T);
    while(T--){
        read3(n,m,p);
        if(p >= M){
            printf("%I64d
",1LL*n*m);
            continue;
        }
        ll ans = 0;
        if(n > m) swap(n,m);
        for(i = 1;i <= n;i = j + 1){//分块处理
            j = min(n/(n/i),m/(m/i));//j代表的是值，并不是公约数；　　
            ans += 1LL*(mu[j][p] - mu[i-1][p])*(n/i)*(m/i);
        }
        printf("%I64d
",ans);
    }
    return 0;
}