HDU 4746 Mophues（莫比乌斯反演）

这道题是跟着秦总博客学的：附上链接：秦总博客

用$f(d)$表示满足$d=gcd(x,y)$且$1<=x<=n,1<=y<=m$的对数。

用$F(d)$表示满足$d|gcd(x,y)$且$1<=x<=n,1<=y<=m$的对数。

可知：

$F(d)=\lfloor\frac{n}{d}\rfloor*\lfloor\frac{m}{d}\rfloor$

$f(x)=\sum_{x|d}\mu(\frac{d}{x})F(d)$

设$k$是满足质因子个数小于$p$的数。枚举每一个$k$，$\sum_{k}f(k)$就是所求的答案。

$ans=\sum_{k}f(k)=\sum_{k}\sum_{k|d}\mu(\frac{d}{k})F(d)$

展开：

$ans=\sum_{k}f(k)=\sum_{k}\sum_{k|d}\mu(\frac{d}{k})\lfloor\frac{n}{d}\rfloor*\lfloor\frac{m}{d}\rfloor$

算一下现在的复杂度：

$o(n)$线筛出$\mu(x)$，枚举每一个$k$，因为对于每一个$k$对应的满足条件的

$d$有$\frac{n}{k}$个。

一个事实：$o(\sum_{i=1}^n{\frac{1}{i}})$复杂度为$o(logn)$

所以每一个询问复杂度是$o(n*logn)$，总的复杂度就为$o(q*n*logn+n)$，可以看出来跪了。

每次询问能接受的复杂度应该是$o(logn)$或者是$o(\sqrt{n})$级别的。

再来看表达式：

$ans=\sum_{k}f(k)=\sum_{k}\sum_{k|d}\mu(\frac{d}{k})\lfloor\frac{n}{d}\rfloor*\lfloor\frac{m}{d}\rfloor$

可以推出：

$ans=\sum_{d}\lfloor\frac{n}{d}\rfloor*\lfloor\frac{m}{d}\rfloor*\sum_{k|d}\mu(\frac{d}{k})$

预处理出$sum(d)=\sum_{k|d}\mu(\frac{d}{k})$，这里还有一个问题：$k$的质因子树要保证小于$p$个。

$1...5*10^5$质因子数最多有$log(5*10^5)$个，即不超过19个，这样原来的$sum(d)$再加上一维：$sum(p,d)$即可。

此时，$\sum_{k|d}\mu(\frac{d}{k})$在每次查找中的复杂度变成了$o(1)$，现在的问题变成了如何在合适的复杂度内找出$\sum_{d}\lfloor\frac{n}{d}\rfloor*\lfloor\frac{m}{d}\rfloor$的值？

这里用到分块加速的方法：

因为$\lfloor{{n}\over{i}}\rfloor$是向下取整，对于一个$i$我们一定能找到一个区间$[i,i']$使得所有$\lfloor{{n}\over{i}}\rfloor$在该区间内取值都相同。

例：
设$p=\lfloor{{n}\over{i}}\rfloor$，
当$p<={{n}\over{j}}<p+1$时$\lfloor{{n}\over{j}}\rfloor$的取值都是$p$。
$p<={{n}\over{j}}<p+1\implies{\lfloor{{n}\over{p+1}}\rfloor}<j<={\lfloor{{n}\over{p}}}\rfloor$

枚举$i$只需要找到右端点即可，所有我们只需要$j$的右端点。即：

$i'={\lfloor{n\over{\lfloor{n\over{i}}\rfloor}}\rfloor}$

原式中是$\sum_{d}\lfloor\frac{n}{d}\rfloor*\lfloor\frac{m}{d}\rfloor$，每次右端点都取$min(\lfloor{n\over{\lfloor{n\over{i}}\rfloor}}\rfloor,\lfloor{m\over{\lfloor{m\over{i}}\rfloor}}\rfloor)$即可

这样在求和中不需要$o(n)$的遍历，看别人的博客说复杂度是$o(\sqrt{n})$，虽然我不会证。

此时$\sum_{d}\lfloor\frac{n}{d}\rfloor*\lfloor\frac{m}{d}\rfloor$的复杂度降到了$o(\sqrt{n})$。每次$sum(p,d)=\sum_{k|d}\mu(\frac{d}{k})$不再是一个完整区间的和，再维护一个前缀后即可在$o(1)$复杂度内搞定$\sum_{k|d}\mu(\frac{d}{k})$。$sum(p,d)$的值需要$o(nogn)$复杂度搞定。

所以总的复杂度就是$o(q*\sqrt{n}+n+nlogn)$

下面是代码部分：

#include <bits/stdc++.h>
#define LL long long
#define FOR(i,x,y)  for(int i = x;i < y;++ i)
#define IFOR(i,x,y) for(int i = x;i > y;-- i)

using namespace std;

const int maxn = 500050;

int prime[maxn],mu[maxn],cnt[maxn];
int sum[20][maxn];
bool check[maxn];

void Mobius(){
    memset(check,false,sizeof(check));
    mu[1] = 1;
    prime[0] = 0;
    cnt[1] = 0;
    FOR(i,2,maxn){
        if(!check[i]){
            prime[++prime[0]] = i;
            mu[i] = -1;
            cnt[i] = 1;
        }
        FOR(j,1,prime[0]+1){
            if(i*prime[j] >= maxn)  break;
            check[i*prime[j]] = true;
            cnt[i*prime[j]] = cnt[i] + 1;
            if(i % prime[j]){
                mu[i*prime[j]] = -mu[i];
            }
            else{
                mu[i*prime[j]] = 0;
                break;
            }
        }
    }
}

int n,m,p;

void init(){
    Mobius();
    memset(sum,0,sizeof(sum));
    //sum[i][j]表示所有能整除j的质因子个数为i的(x,y)的对数
    FOR(i,1,maxn){
        for(int j = i;j < maxn;j += i){
            sum[cnt[i]][j] += mu[j/i];
        }
    }
    //sum[i][j]表示所有能整除j的质因子个数<=i的(x,y)的对数
    FOR(j,1,maxn){
        FOR(i,1,20){
            sum[i][j] += sum[i-1][j];
        }
    }
    //sum[i][j]表示sum[i][0]+sum[i][1]+...+sum[i][j]的前缀和
    FOR(j,1,maxn){
        FOR(i,0,20){
            sum[i][j] += sum[i][j-1];
        }
    }
}

void work(){
    if(p >= 20) {printf("%I64d\n",(LL)n*(LL)m);return;}
    int mx = min(n,m);
    LL ans = 0;
    for(int i = 1;i <= mx;){
        LL a = n/i,b = m/i;
        LL now = min(n/a,m/b);
        ans += a*b*(sum[p][now]-sum[p][i-1]);
        i = now+1;
    }
    printf("%I64d\n",ans);
}

int main(){
    //freopen("test.in","r",stdin);
    init();
    int T;  scanf("%d",&T);
    while(T--){
        scanf("%d%d%d",&n,&m,&p);
        work();
    }
    return 0;
}