P3327 莫比乌斯反演+整除分块

Tenshi · 2021-06-01 13:11:15 · 题解

莫比乌斯反演 + 整除分块

分析

首先，我们给出一个结论：

证明：

设 i 的分解式为 i = \prod p_k^{\alpha_k} ，类似地，j = \prod p_k^{\beta_k}

对于质数 p_k ，ij 对应的因数个数为 \alpha_k + \beta_k + 1

而对于右式，p_k 如果想要产生贡献，那么只可能是：x, y 不同时存在因数 p_k，注意到对应的贡献数即为 \alpha_k + \beta_k + 1

因此，由乘法原理，所求证式成立。

为了防止混淆，我们将题面中的 n, m 记为 N, M

由上述结论，题目的式子转化为：

莫比乌斯反演：若 F(n) = \sum_{n|d}f(d)，则有 f(n) = \sum_{n|d} \mu(\frac{d}{n})F(d)

于是我们设 f(n) = \sum_{i=1}^N \sum_{j=1}^M \sum_{x|i} \sum_{y|j} [(x, y) = n] ，对应的 F(n) = \sum_{i=1}^N \sum_{j=1}^M \sum_{x|i} \sum_{y|j} [n | (x, y)]

下面对 F(n) 进行变换：

设 x' = \lfloor \frac{x}{n} \rfloor ~, y'=\lfloor \frac{y}{n} \rfloor ~ , N'=\lfloor \frac{N}{n} \rfloor ~ , M'=\lfloor \frac{M}{n} \rfloor ，进而有

记 h(x) = \sum_{i=1}^x \lfloor \frac{x}{i} \rfloor ，h(x) 可用整除分块处理。

由莫比乌斯反演，f(n) = \sum_{n|d} \mu(\frac{d}{n})F(d) ，我们只需求 f(1)

下面对 f(1) 进行变换：

由整除分块知，h(\lfloor \frac{N}{d} \rfloor)h(\lfloor \frac{M}{d} \rfloor) 取值最多为 \sqrt{N} + \sqrt{M} 块，而对于每一块，可以用前缀和 O(1) 处理出对应的 \mu 值，因此整体的复杂度为 O(T(\sqrt{N} + \sqrt{M}))

细节见代码：

#include<bits/stdc++.h>
using namespace std;

const int S=50050;

int primes[S], cnt;
bool vis[S];
int mu[S], sum[S];
int h[S];

int g(int b, int l){
    return b/(b/l);
}

void init(){
    // init sum[]
    mu[1]=1;
    for(int i=2; i<S; i++){
        if(!vis[i]) primes[cnt++]=i, mu[i]=-1;
        for(int j=0; i*primes[j]<S; j++){
            vis[i*primes[j]]=true;
            if(i%primes[j]==0) break;
            mu[i*primes[j]]=-mu[i];
        }
    }
    for(int i=1; i<S; i++) sum[i]=sum[i-1]+mu[i];

    // init h[]
    for(int x=1; x<S; x++){
        int &v=h[x];
        for(int l=1, r; l<=x; l=r+1){
            r=min(x, g(x, l));
            v+=x/l*(r-l+1);
        }
    }
}

int solve(int N, int M){
    int res=0;

    int n=min(N, M);
    for(int l=1, r; l<=n; l=r+1){
        r=min(n, min(g(N, l), g(M, l)));
        res+=(sum[r]-sum[l-1])*h[N/l]*h[M/l];
    }
    return res;
}

int main(){
    int T; cin>>T;
    init();
    while(T--){
        int N, M; cin>>N>>M;
        cout<<solve(N, M)<<endl;
    }
    return 0;
}