题目链接：Project Euler 745

定义 $g(n)$ 表示最大的能够整除 $n$ 的完全平方数。例如 $g(18) = 9, g(19) = 1$。

定义 $S(n) = \sum_{i = 1}^{n} g(i)$，求 $S(10^{14}) \bmod (10^{9} + 7)$

Solution

设 $g(x) = i^{2}$，由于 $\frac{x}{i^{2}}$ 不能含有平方因子，可以枚举 $i$ 并枚举 $\frac{x}{i^{2}}$ 的值，得到：

$$ \sum_{i = 1}^{\lfloor\sqrt{n}\rfloor} i^{2} \sum_{j = 1}^{\left\lfloor\frac{n}{i^{2}}\right\rfloor} \mu^{2}(j) $$

后面的 $\sum$ 本质上是 $1 \sim \left\lfloor\frac{n}{i^{2}}\right\rfloor$ 中没有平方因子的数字个数，设其为 $h\left(\left\lfloor\frac{n}{i^{2}}\right\rfloor\right)$，可以容斥求 $h(n)$：

$$ \begin{align*} h(n) = n - \sum_{i} \left\lfloor\frac{n}{p_{i}^{2}}\right\rfloor + \sum_{i < j} \left\lfloor\frac{n}{(p_{i} p_{j})^{2}}\right\rfloor - \ldots \end{align*} $$

设 $(p_{i_{1}} p_{i_{2}} p_{i_{3}} \ldots p_{i_{k}})^{2} = x^2$，则 $\left\lfloor\frac{n}{x^{2}}\right\rfloor$ 项的系数为 $(-1)^{k}$ 且 $x$ 满足无平方因子，式子改写为：

$$ h(n) = \sum_{i = 1}^{\lfloor\sqrt{n}\rfloor} \left\lfloor\frac{n}{i^{2}}\right\rfloor \mu(i) $$

代入原式：

$$ \begin{align*} & \sum_{i = 1}^{\lfloor\sqrt{n}\rfloor} i^{2} h\left(\left\lfloor\frac{n}{i^{2}}\right\rfloor\right) \\ = & \sum_{i = 1}^{\lfloor\sqrt{n}\rfloor} i^{2} \sum_{j = 1}^{\left\lfloor\sqrt{\frac{n}{i^{2}}}\right\rfloor} \left\lfloor\frac{n}{(ij)^{2}}\right\rfloor \mu(j) \\ = & \sum_{p = 1}^{\lfloor\sqrt{n}\rfloor} \left\lfloor\frac{n}{p^{2}}\right\rfloor \sum_{i \mid p} \mu(i) \left(\frac{p}{i}\right)^{2} \end{align*} $$

注意到后面的 $\sum$ 是一个积性函数，记为 $f(p)$，则有：

$$ \begin{align*} f(n) & = \sum_{i \mid n} \mu(i) \left(\frac{n}{i}\right)^{2} \\ f(p) & = p^{2} - 1 \\ f(p^{k}) & = p^{2k} - p^{2k - 2} \end{align*} $$

时间复杂度：$\mathcal O(\sqrt{n})$，其中 $n = 10^{14}$。

Code

答案：$94586478$。

#include <bits/stdc++.h>
typedef long long int64;

const int N = 1e7 + 5, MOD = 1e9 + 7;
int n, tot, prime[N], f[N];

inline void add(int &x, const int y) {
    (x += y) >= MOD && (x -= MOD);
}
void sieve(int n) {
    std::bitset<N> flg;
    f[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!flg[i]) {
            prime[++tot] = i;
            f[i] = (int64)i * i % MOD - 1;
        }
        for (int j = 1; j <= tot && i * prime[j] <= n; j++) {
            flg[i * prime[j]] = 1;
            if (i % prime[j] == 0) {
                f[i * prime[j]] = (int64)f[i] * prime[j] % MOD * prime[j] % MOD;
                break;
            }
            f[i * prime[j]] = (int64)f[i] * ((int64)prime[j] * prime[j] % MOD - 1) % MOD;
        }
    }
}
int main() {
    int64 n = 1e14;
    sieve((int)sqrt(n));
    int ans = 0;
    for (int i = 1; (int64)i * i <= n; i++) {
        add(ans, (int64)(n / ((int64)i * i)) % MOD * f[i] % MOD);
    }
    printf("%d\n", ans);
    return 0;
}