Luogu 4844 LJJ爱数数

LOJ 6482

设$d = gcd(a, b)$，$xd = a$，$yd = b$，因为$frac{1}{a} + frac{1}{b} = frac{a + b}{ab} = frac{1}{c}$，所以$c(x + y)= xyd$。

因为$d$不整除于$c$，那么$d | (x + y)$，把$d$除过去，

$$frac{x + y}{d} = frac{xy}{c}$$

设这个式子等于$p$，如果$p$不为$1$，那么$p | x$或者$p | y$，$p$不可能同时整除$x, y$（$x, y$互质），但是$p | (x + y)$，所以不成立，得到$p = 1$。

原来的条件就变成了$a + b = d^2$，$c = frac{ab}{d^2}$，$1 leq a,b,c leq n$。

我们可以枚举这个$d$，然后枚举$x$，这样子只要算出和$d$互质的数的数量就可以了。

考虑一下$x$的范围，$x$满足

$$1 leq dx leq n$$

$$1 leq d^2 - dx leq n$$

那么$max(1, d - left lfloor frac{n}{d} ight floor)leq x leq min(d - 1, left lfloor frac{n}{d} ight floor)$。

即为求

$$sum_{d = 1}^{sqrt{2n}}sum_{x = max(1, d - left lfloor frac{n}{d} ight floor)}^{min(d - 1, left lfloor frac{n}{d} ight floor)}[gcd(x, d - x) == 1]$$

有辗转相减法$gcd(x, d - x) = gcd(x, d)$。

求一定范围内的$x$，把它拆成前缀和相减的形式，现在就是要做：

$$sum_{i = 1}^{k}[gcd(i, d) == 1]$$

$$ = sum_{i = 1}^{k}sum_{j | gcd(i, d)}mu (j)$$

$$ = sum_{i = 1}^{k}sum_{j = 1}^{d}mu (j)[j | d][j | i]$$

$$ = sum_{j | d}mu (j)left lfloor frac{k}{j} ight floor$$

枚举$d$的约数即可。

时间复杂度$O(sqrt{n}log(sqrt{n}))$。

Code：

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

const int N = 2e6 + 5;
const int Maxn = 1.5e6;
const int M = 2e7 + 5;

int pCnt = 0, pri[N], mu[N], tot = 0, head[N];
ll n;
bool np[N];

struct Node {
    int to, nxt;
} e[M];

inline void add(int from, int to) {
    e[++tot].to = to;
    e[tot].nxt = head[from];
    head[from] = tot;
}

inline void sieve() {
    mu[1] = 1;
    for (int i = 2; i <= Maxn; i++) {
        if (!np[i]) pri[++pCnt] = i, mu[i] = -1;
        for (int j = 1; j <= pCnt && i * pri[j] <= Maxn; j++) {
            np[i * pri[j]] = 1;
            if (i % pri[j] == 0) {
                mu[i * pri[j]] = 0;
                break;
            }
            mu[i * pri[j]] = -mu[i];
        }
    }
}

inline ll max(ll x, ll y) {
    return x > y ? x : y;
}

inline ll min(ll x, ll y) {
    return x > y ? y : x;
}

inline ll getSum(ll k, int g) {
    ll res = 0;
    for (int i = head[g]; i; i = e[i].nxt)
        res += mu[e[i].to] * (k / e[i].to);
    return res;
}

int main() {
    sieve();
    scanf("%lld", &n);
    
    int lim = sqrt(2LL * n);
    for (int i = 1; i <= lim; i++)
        if (mu[i])
            for (int j = 1; i * j <= lim; j++)
                add(i * j, i);
                    
    ll ans = 0LL;
    for (int i = 1; i <= lim; i++) {
        ll mx = min(n / i, 1LL * i - 1), mn = max(1LL, 1LL * i - (n / i));
        ans += getSum(mx, i) - getSum(mn - 1, i);
    }
    
    printf("%lld
", ans);
    return 0;
}