1. 程式人生 > >BZOJ 2820 YY的GCD 莫比烏斯反演

BZOJ 2820 YY的GCD 莫比烏斯反演

cstring std swap cst pac lap Go hint name

2820: YY的GCD

Description

神犇YY虐完數論後給傻×kAc出了一題給定N, M,求1<=x<=N, 1<=y<=M且gcd(x, y)為質數的(x, y)有多少對kAc這種 傻×必然不會了,於是向你來請教……多組輸入

Input

第一行一個整數T 表述數據組數接下來T行,每行兩個正整數,表示N, M

Output

T行,每行一個整數表示第i組數據的結果

Sample Input

2
10 10
100 100

Sample Output

30
2791

HINT

T = 10000

N, M <= 10000000

思路:

題目中描述的式子即為:

$\sum\limits_{p}is[p]\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{M}[gcd(i,j)=p]$

把$p$除到前面可以化簡得 :

$\sum\limits_{p}is[p]\sum\limits_{i=1}^{\lfloor\frac {N}{p}\rfloor}\sum\limits_{j=1}^{\lfloor\frac {M}{p}\rfloor}[gcd(i,j)=1]$

出現$[gcd(i,j)=1]$的形式,考慮莫比烏斯反演,化簡得:

$\sum\limits_{p}is[p]\sum\limits_{i=1}^{\lfloor\frac {N}{p}\rfloor}\sum\limits_{j=1}^{\lfloor\frac {M}{p}\rfloor}\sum\limits_{d|gcd(i,j)}\mu(d)$


把d提到前面得

$\sum\limits_{p}is[p]\sum\limits_{d=1}^{\lfloor\frac {\min(N,M)}{p}\rfloor}\mu(d)\sum\limits_{i=1}^{\lfloor\frac {N}{dp}\rfloor}\sum\limits_{j=1}^{\lfloor\frac {M}{dp}\rfloor}$

設$Q = dp$ 枚舉$Q$化簡為

$\sum\limits_{Q=1}^{\min(N,M)}\lfloor\frac{N}{Q}\rfloor \lfloor\frac{M}{Q}\rfloor\sum\limits_{p|Q}is[p]\mu(\frac{Q}{p})$

設函數$f(n) = \sum\limits_{p|n}is[p]\mu(\frac{n}{p})$

考慮得出$f(n)$

有以下幾種情況 :
1. 若 $f(n)$ 為質數
值即為$\mu(1) = 1$
2. 若 $n % p == 0$ 則$f(n \times p)$可以化成$\sum\limits_{d|n\times p}is[d]\mu(\frac{n\times p}{d})$

考慮當$d!=p$時$\frac{n\times p}{d}$有多個$p$
對$\sum$的貢獻為0,所以此時$f(n\times p)=\mu(n)$
1. 若 $n % p != 0$ , $f(n\times p)$可以化為$f(n)\times \mu(p) + f(p)\times \mu(n)$ 。我們又知道$\mu(p) = -1$,$f(p) = 1$,所以$f(n \times p)=\mu(n)-f(n)$

我們再處理一下前綴和, 老套路分$\sqrt n$塊計算。得到結果
```cpp

#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;
const int N = 11000000;
bool np[N+10];
int mu[N+10], tot, pr[N+10], sum[N+10];
int f[N+10];
int t, n, m;
void init() {
    mu[1]=1;
    for(int i=2;i<=N;i++) {
        if(!np[i]) {
            pr[++tot]=i;
            mu[i]=-1;
            f[i]=1;
        }
        for(int j=1;j<=tot&&i*pr[j]<=N;j++) {
            np[i*pr[j]]=1;
            if(!(i%pr[j])) {
                mu[i*pr[j]]=0;
                f[i*pr[j]]=mu[i];
                break;
            }
            mu[i*pr[j]]=-mu[i];
            f[i*pr[j]]=mu[i]-f[i];
        }
        sum[i]=sum[i-1]+f[i];
    }
}
void solve(int n, int m) {
    if(n>m)swap(n,m);
    int lst;
    long long ans=0;
    for(int i=1;i<=n;i=lst+1) {
        lst = min(n/(n/i), m/(m/i));
        ans+=1ll*(sum[lst]-sum[i-1])*(n/i)*(m/i);
    }
    printf("%lld\n", ans);
}
int main() {
    init();
    scanf("%d",&t);
    while(t--) {
        scanf("%d%d",&n,&m);
        solve(n,m);
    }
}

BZOJ 2818 雙倍經驗

```

BZOJ 2820 YY的GCD 莫比烏斯反演