HDU 1695 莫比乌斯反演

莫比乌斯反演

对于定义在非负整数上的两个函数F(x)， f(x) : 若
$F(n) = \sum_{d | n} f(d)$
则

$$f(n) = \sum_{d | n} u(d)F(\frac{n}{d})　(1)$$
其中： $u(d)$就是莫比乌斯函数, 它的定义如下
$$u(d) = \begin {cases} 1, & d = 1 \\ (-1)^k, & d = p_1p_2p_3···p_k, p_i 是质数，i\epsilon[1, k]\\\ 0, & 其它 \end {cases} (2)$$

常见性质 ：

$$\sum_{d | n} u(d) = \begin {cases} 1, & n = 1\\\ 0, & n >１ \end {cases}(3)$$
证明 : (1)
$$\sum_{d | n} u(d)F(\frac{n}{d})=\sum_{d | n}u(d)\sum_{d' | \frac{n} {d}}f(d')\\\ = \sum_{d' | \frac{n}{d} }u(d')\sum_{d | n}f(d)根据(3)\\\ = \sum_{d'|1}u(d')f(n) = f(n)得证$$

还有一种形式的反演变换:
若$F(n) = \sum_{n| d} f(d)$　

$$f(n) = \sum_{n| d} u(\frac{d}{n})F(d)证明和(1)类似$$

典型的反演:
$n = \sum_{d|n}\phi(d)$
$\phi(n) = \sum_{d|n}u(d) * \frac{n}{d}$

/*莫比乌斯函数线性筛法*/
int mobius[N], prime[N];
bool vis[N];
void getMobius(){
    memset(vis, 0, sizeof vis);
    memset(mobius, 0, sizeof mobius);
    mobius[1] = 1;
    int cnt = 0;
    for(int i = 2; i < N; ++i){
        if(!vis[i]) {
            mobius[i] = -1;/*质数**/
            prime[++cnt] = i;
        }
        for(int j = 1; j <= cnt && i * prime[j] < N; ++j){
            vis[prime[j] * i] = 1;
            if(i % prime[j] == 0) {
                mobius[i * prime[j]] = 0;
/**某个质因子出现两个以上*/
                break;
            }
            mobius[i * prime[j]] =  mobius[i];
          }
    }
}

题解

分别在[1, n]和[1, m]中找出一对x, y满足gcd(x, y) = k
就是求在[1, n / k]和[1,m / k]中找出一对x’, y’满足gcd(x’, y’) = 1
暴力枚举极端情况O(n^2)n最大100000会T
所以定义$f(d)为[1, N]和[1, M]中gcd(x, y) = d的数对$
在定义$F(x)为[1, N]和[1, M]中gcd(x, y)是d的倍数的数对$

明显$F(x) = N / x * M / x$
且有公式$F(x) = \sum_{x | d}f(d)$

根据莫比乌斯反演， 得：
$f(x) = \sum_{x | d}u(\frac{d}{x})F(d)$
因此

$$f(1) = \sum_{i = 1}^{min(N,M)} u(i) F(i)$$
然后在去掉重复的对数：[1, min(N, M)]有一半重复就可解

code

#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>
using namespace std;
typedef long long ll;

const int N = 100005;

int mobius[N], prime[N];
bool vis[N];
int a, b, c, d, k;

void init(){
    memset(vis, 0, sizeof vis);
    memset(mobius, 0, sizeof mobius);
    mobius[1] = 1;
    int cnt = 0;
    for(int i = 2; i < N; ++i){
        if(!vis[i]) {
            mobius[i] = -1;
            prime[++cnt] = i;
        }
        for(int j = 1; j <= cnt && i * prime[j] < N; ++j){
            vis[prime[j] * i] = 1;
            if(i % prime[j] == 0) {
                mobius[i * prime[j]] = 0;
                break;
            }
            mobius[i * prime[j]] = -mobius[i];
        }
    }

}

inline void input(){
    cin >> a >> b >> c >> d >> k;
}

void solve(){
    if(b < d) swap(b, d);
    b /= k;
    d /= k;
    ll ans1 = 0, ans2 = 0;
    for(int i = 1; i <= d; ++i){
        ans1 += mobius[i] * 1ll * (b / i) * (d / i);
        ans2 += mobius[i] * 1ll * (d / i) * (d / i);
    }

    cout << ans1 - (ans2 >> 1) << endl;
}

int main(){

    int cas = 1, t;
    cin >> t;
    init();

    while(t--){
        input();
        cout << "Case " << cas++ << ": ";
        if(k == 0) {
                cout << "0" << endl;continue;
        }
        solve();
    }

    return 0;
}