luoguP4466 [国际集训队]和与积 莫比乌斯反演

Posted 2021-01-27 reverymoon

---

自然想到枚举(gcd(a, b))，不妨设其为(d)，并且(a = di, b = dj(a > b))

那么(frac{ab}{a + b} = frac{dij}{i + j})

由于此时有((i,j) = 1)，因此((i, i + j) = (j, i + j) = 1)

那么，当且仅当(i + j | d)时，((i, j))数对对答案有贡献

对答案有多少的贡献呢？(frac{n}{i(i + j)}) 没有想到这一步

理由是(d = k(i + j))，那么只需满足(ki(i + j) leq n)

当(i > sqrt n)时，((i,j))对答案绝对没有贡献

所以答案为(sum limits_{d = 1}^{sqrt n} sum limits_{i = 1}^{sqrt n / d} sum limits_{j} [(i, j) = 1]frac{n}{i(i + j)})

莫比乌斯反演，得到

(sum limits_{d = 1}^{sqrt n} mu(d) sum limits_{i = 1}^{sqrt n / d} sum limits_{j} frac{n}{d^2i(i + j)})

对内层数论分块统计答案即可一开始把不分块的复杂度算错了，以为能过

---

分析一下复杂度上界

首先考虑对于确定的(d)，枚举(i, j)的复杂度

(sum limits_{i = 1}^{sqrt n / d} sqrt frac{n}{d^2 i} = frac{sqrt n}{d} * sum limits_{i = 1}^{sqrt n / d} frac{1}{sqrt i})

用归纳法可以证明，右边那个东西(leq 2 sqrt {sqrt n / d})

所以对于一个(d)而言，需要(frac{n^{frac{3}{4}}}{d^{frac{3}{2}}})的复杂度

由于(frac{1}{1} + frac{1}{2^{1.5}} + ... + frac{1}{n^{1.5}} leq 3)

所以复杂度就是(O(n^{frac{3}{4}}))

然后绝对跑不到这个上界....

---

#include <cmath>
#include <ctime>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;

#define ri register int
#define ll long long
#define rep(io, st, ed) for(ri io = st; io <= ed; io ++)
#define drep(io, ed, st) for(ri io = ed; io >= st; io --)

const int sid = 1e5 + 5;
    
ll ans;
int n, sq, tot;
int mu[sid], nop[sid], pr[sid];
inline void Sieve(int m) {
    mu[1] = 1;
    rep(i, 2, m) {
        if(!nop[i]) {
            pr[++ tot] = i;
            mu[i] = -1;
        }
        rep(j, 1, tot) {
            int p = i * pr[j];
            if(p > m) break; nop[p] = 1;
            if(i % pr[j] == 0) break;
            mu[p] = -mu[i];
        }
    }
}


int main() {
    cin >> n;
    sq = sqrt(n) + 1;
    Sieve(sq);
    rep(d, 1, sq) {
        if(!mu[d]) continue;
        int p = d * d;
        rep(i, 1, sq / d) {
            int fs = n / d / d / i;
            for(ri ii = i + 1, jj; ii <= (i << 1) - 1 && ii <= fs; ii = jj + 1) {
                jj = min(fs / (fs / ii), (i << 1) - 1);
                ans += 1ll * mu[d] * (jj - ii + 1) * (fs / ii);
            }
        }
    }
    printf("%lld
", ans);
    return 0;
}