蓝桥杯2023 C++大学A组 J.翻转硬币 题解

· · 题解

前置知识点

本题是数论题，需要对初等数论以下几个知识有所了解：

以上知识点的学习可参考 OI Wiki。

题解

分析

题意分析

题意即每次选择一个数 i，将 i 的全部倍数翻转。容易看出，一个数翻转两次等价于不翻转，所以每个数至多被选择一次。

同时可以看出，每个硬币只可能被它的因数（包括它本身）翻转，例如硬币 12 只能通过 1, 2, 3, 4, 6, 12 翻转。

初始时硬币 1 朝下，因为 1 的因数只有它自己，所以必须翻转 1 的倍数，此时除 1 之外所有硬币朝下。

然后考虑翻转编号为质数 p 的硬币，p 的因数只有 1 和 p，1 已经选择过，不可能再选择一次，所以必定选择一次质数 p 以翻转。

然后考虑 n = p_1 p_2，也就是可以分解为两个质数乘积的合数。n 已经被 1, p_1, p_2 翻转一次，此时该硬币朝下，剩余的因数只有 n 本身，所以必定选择 n。

继续考虑 n = p_1 p_2 p_3，n 已经被 1, p_1, p_2, p_3, p_1 p_2, p_1 p_3, p_2 p_3 （共 7 个数）翻转，硬币朝下，剩余的因数同样只有 n 本身，所以必定选择 n。

以此类推，所有可以分解为若干个互异质因数之积的整数 n = p_1 p_2 \cdots p_k 都需要选择一次。这样的 n （包括 1）被称为无平方因子数 。

下面证明除此以外的其他数（即有平方因子数）都已经翻转到正面。设有平方因子数 n 的所有互异质因数为 p_1, \dots, p_k，容易看出 n \ne p_1 \cdots p_k（否则和有平方因子矛盾，有平方因子代表存在质数 p 满足 p^2 \mid n），则 n 已经被 1, p_1, \cdots, p_k, p_1 p_2, \cdots, p_{k - 1} p_k, \cdots, p_1\cdots p_k 无平方因子数筛去。n 的因数中，无平方因子数的个数为 2^k（等价于 \{p_1, \dots, p_k\} 的子集数），2^k 是偶数，且最大数 n \ne p_1 \cdots p_k，所以硬币 n 已经被翻转了偶数次，依然是正面。

因此只需要统计 [1, n] 中的无平方因子数个数。无平方因子数的 Möbius 函数值为 \pm 1，其他数为 0，因此问题等价于求 Möbius 函数的绝对值（或平方）的前缀和：

\sum_{i = 1}^n \mu^2(k).

\mu^2 的前缀和性质

首先证明 Möbius 函数平方的一个重要性质：

\mu^2(n) = \sum_{d^2 \mid n} \mu(d).

证明：设 n 的标准分解式为 p_1^{\alpha_1} \cdots p_k^{\alpha_k}，则 d^2 \mid n 说明 d 的标准分解式只能包含次数 \alpha_i \ge 2 的质因数 p_i。

设 I(k) = [k = 1]，n' = p_1^{\lfloor \alpha_1 / 2 \rfloor} \cdots p_k^{\lfloor \alpha_k / 2 \rfloor}，则 d \mid n'，于是

\sum_{d^2 \mid n} \mu(d) = \sum_{d \mid n'} \mu(d) = I(n')

最后一步运用了经典 Dirichlet 卷积结论 \sum_{d \mid n} \mu(d) = I(n)。

而 I(n') = 1 则说明 n' = 1，\lfloor \alpha_1 / 2 \rfloor = \cdots = \lfloor \alpha_k / 2 \rfloor = 0，也就是 \alpha_1 = \cdots = \alpha_k = 1，对应地 \mu^2(n) = 1。所以 I(n') = \mu^2(n)，命题证毕。

然后据此推导 Möbius 函数平方的一个重要前缀和性质：

\sum_{k = 1}^n \mu^2(k) = \sum_{k = 1}^{\lfloor \sqrt{n} \rfloor} \left \lfloor \frac{n}{k^2} \right \rfloor \mu(k).

证明：将 \mu^2 (k) 用上一定理展开，然后交换求和顺序即可。交换求和顺序的方法参考自 https://zhuanlan.zhihu.com/p/499839696 。

\begin{aligned} \sum_{k = 1}^n \sum_{d^2 \mid k} \mu(d) &= \sum_{d = 1}^{\lfloor \sqrt{n} \rfloor} \sum_{d^2 \mid k, k \le n} \mu(d) \\ &= \sum_{d = 1}^{\lfloor \sqrt{n} \rfloor} \mu(d) \sum_{d^2 \mid k, k \le n} 1 \\ &= \sum_{d = 1}^{\lfloor \sqrt{n} \rfloor} \mu(d) \left \lfloor \frac{n}{d^2} \right \rfloor. \end{aligned}

到这一步，如果使用 Euler 筛预处理出 [1, \sqrt{n}] 的 Möbius 函数值，然后直接计算上述和式，已经可以获得 70% 的分数。但是，上述算法对 n = 10^{18} 的情况依然会超时。

数论分块

我们需要进一步处理，继续降低 n 达到 10^{18} 时的复杂度。观察到和式中包含了 \left \lfloor \frac{n}{d^2} \right \rfloor，可以考虑使用数论分块。

数论分块 \sum f(k) g \left( \left \lfloor \frac{n}{k} \right \rfloor \right) 要求能够预处理出 f(n) 的前缀和，这里 \mu(n) 的前缀和可以用杜教筛来快速计算。但这里 \left \lfloor \frac{n}{d^2} \right \rfloor 的形式和普通的数论分块并不一致，需要推导新的数论分块公式。

数论分块的关键是，对于给定的 l，确定 r = \max \{i \in \mathbb{Z}: \left \lfloor \frac{n}{l^2} \right \rfloor = \left \lfloor \frac{n}{i^2} \right \rfloor\} 的值。设 y = \left \lfloor \frac{n}{l^2} \right \rfloor，这里需要确定的就是满足 \left \lfloor \frac{n}{x^2} \right \rfloor = y 的最大整数 x。容易看出 y \le \frac{n}{x^2} < y + 1，因此 \frac{n}{y + 1} < x^2 \le \frac{n}{y}，x 的最大值即为 \left \lfloor \sqrt{\frac{n}{y}} \right \rfloor。

因此，仿照普通的数论分块，这里的数论分块代码如下：

ll prefix_mu2(size_t n)
{
    // lower_sqrt 表示开方后向下取整
    size_t l = 1, r, bound = lower_sqrt(n);
    ll ans = 0;
    while (l <= bound)
    {
        r = lower_sqrt(n / (n / (l * l)));
        ans += (n / (l * l)) * (prefix_mu(r) - prefix_mu(l - 1)); // prefix_mu 表示 Möbius 函数前缀和
        l = r + 1;
    }
    return ans;
}

其中 lower_sqrt 函数可以用二分法实现。

\mu 前缀和 & 杜教筛

最后只需要计算 k \le 10^9 时 \mu(k) 的前缀和即可。

可以预先用 Euler 筛求出 k \le 10^7 的 \mu 值与前缀和，再用杜教筛处理 k > 10^7 的前缀和。

复杂度分析

杜教筛复杂度

根据 OI-Wiki 杜教筛 的结论，若预处理出前 m^\alpha 的 \mu 前缀和，单次使用杜教筛计算 \mu(i) 的前 m 项和的时间复杂度上界为 O(m^{1 - \alpha})。

这里由于 m \le 10^9，至少应该预处理出前 10^6 的前缀和，否则单次使用杜教筛的时间复杂度将大于预处理复杂度。

数论分块的分块数

由于这里使用了新的数论分块技巧，仿照经典数论分块的方法，这里 \left \lfloor \frac{n}{d^2} \right \rfloor 的取值不会超过 2n^{1 / 3} 个数，从而分块数不会超过 2n^{1 / 3}，证明如下：

每个分块都使用了二分法求开平方值，开平方的总时间复杂度是 O(n ^ {1 / 3} \log n)。

总复杂度

在数论分块中，若预处理出前 n^\alpha 的 \mu 前缀和，则当 r \le n^\alpha 时，计算这个区间的 \mu 前缀和是常数复杂度；当 r > n^\alpha 时，则需要调用杜教筛进行计算。

这里设 \alpha > \frac{1}{3}（在杜教筛复杂度一节，已经至少预处理了 10^6），统计 l > n^\alpha 的区间个数。此时 \left \lfloor \frac{n}{l^2} \right \rfloor 的取值至多有 n^{1 - 2\alpha} 个，也就是至多有 n^{1 - 2\alpha} 个区间。这里每个区间都要调用一次杜教筛，复杂度为 O((\sqrt n)^{1 - 2\alpha}) = O(n^{1 / 2 - \alpha})，所以总的时间复杂度为 O(n^{3 / 2 - 3\alpha})。

由上式看出 \alpha 越接近 1/2 越好，也就是预处理的数越大越好。但由于预处理前 n^{\alpha} 项的的空间复杂度为 O(n^\alpha)，n = 10^{18} 时会 MLE，所以并不能直接令 \alpha = 1 / 2。

如果预处理前 10^7 项，则 \alpha = 7 / 18，时间复杂度为 O(n^{3 / 2 - 3 \times 7 / 18}) = O(n^{1 / 3})。结合前面的开平方复杂度，总的时间复杂度为 O(n^{1 / 3} \log n)，恰好可以通过本题。

内存优化

预处理的数据越多，\alpha 就越大，杜教筛的时间复杂度越低。这里我们将进一步优化内存，使程序的预处理部分达到 5 \times 10^7。

Euler 筛需要一个 is_prime 数组记录被筛去的质数，一个 primes 记录当前得到的质数，还需要一个 mu 数组记录已经筛出的 Möbius 函数值，代码：

constexpr size_t MAXN = 1e7;
int mu[MAXN + 1], smu[MAXN + 1];
bitset<MAXN + 1> is_prime;
unsigned primes[MAXN];

void preprocess()
{
    constexpr size_t n = MAXN;
    size_t cnt = 0;
    is_prime.set();
    is_prime[0] = is_prime[1] = false;
    mu[1] = 1;
    for (ull i = 2; i <= n; i++)
    {
        if  (is_prime[i])
        {
            primes[cnt++] = i;
            mu[i] = -1;
        }
        for (size_t j = 0; j < cnt; j++)
        {
            if ((ll) i * primes[j] > n)
            {
                break;
            }
            is_prime[i * primes[j]] = false;
            if (i % primes[j] == 0)
            {
                mu[i * primes[j]] = 0;
                break;
            }
            else
            {
                mu[i * primes[j]] = -mu[i];
            }
        }
    }
    for (size_t i = 1; i <= n; i++)
    {
        smu[i] = smu[i - 1] + mu[i];
    }
}

但 is_prime、primes 和 mu 数组都只是预处理过程的辅助变量，预处理结束后保留 Möbius 前缀和数组 smu 即可。因此，我们可以使用代码块（一个花括号），将 is_prime 和 primes 变为局部变量，程序跳出代码块时会自动销毁并释放内存。同时，在一开始就将 Möbius 函数值存储在 smu 内，后面直接原地计算前缀和即可。优化后的代码如下：

constexpr size_t MAXN = 5e7;
int smu[MAXN];

void preprocess()
{
    constexpr ull n = MAXN;
    {
        bitset<MAXN + 1> is_prime;
        vector<unsigned int> primes;
        is_prime.set();
        is_prime[0] = is_prime[1] = false;
        smu[1] = 1;
        for (ull i = 2; i <= n; i++)
        {
            if  (is_prime[i])
            {
                primes.push_back(i);
                smu[i] = -1;
            }
            for (unsigned int p: primes)
            {
                if ((ull) i * p > n)
                {
                    break;
                }
                is_prime[i * p] = false;
                if (i % p == 0)
                {
                    smu[i * p] = 0;
                    break;
                }
                else
                {
                    smu[i * p] = -smu[i];
                }
            }
        }
    }
    for (size_t i = 1; i <= n; i++)
    {
        smu[i] += smu[i - 1];
    }
}

这里开 5 \times 10^7 大小的 int 数组，占用内存为 4 \times 5 \times 10^7\ \mathrm{B} = 190.73 \ \mathrm{MB}，内存限制为 256\ \mathrm{MB}，不会 MLE。

代码

汇总以上全部思想，完整的算法步骤如下：

  1. 用 Euler 筛预处理 k \le 10^7 的 Möbius 函数值及其前缀和；
  2. 编写函数，使用杜教筛处理 k > 10^7 的 Möbius 函数前缀和；
  3. 使用数论分块计算 \mu^2 的前缀和。
#include <iostream>
#include <unordered_map>
#include <bitset>
#include <vector>
using namespace std;

using ll = long long;
using ull = unsigned long long;

constexpr size_t MAXN = 5e7;
constexpr size_t INF = 1e9;

int smu[MAXN];

void preprocess()
{
    // Euler seive
    constexpr ull n = MAXN;

    {
        bitset<MAXN + 1> is_prime;
        vector<unsigned int> primes;
        is_prime.set();
        is_prime[0] = is_prime[1] = false;
        smu[1] = 1;
        for (ull i = 2; i <= n; i++)
        {
            if  (is_prime[i])
            {
                primes.push_back(i);
                smu[i] = -1;
            }
            for (unsigned int p: primes)
            {
                if ((ull) i * p > n)
                {
                    break;
                }
                is_prime[i * p] = false;
                if (i % p == 0)
                {
                    smu[i * p] = 0;
                    break;
                }
                else
                {
                    smu[i * p] = -smu[i];
                }
            }
        }
    }

    // Calculate the prefix sum
    for (size_t i = 1; i <= n; i++)
    {
        smu[i] += smu[i - 1];
    }
}

unordered_map<size_t, ll> fmu; // Store calculated prefix sums of mu

ll prefix_mu(size_t n)
{
    if (n <= MAXN)
    {
        return smu[n];
    }
    else if (fmu.find(n) != fmu.end()) // This value has been calculated
    {
        return fmu[n];
    }
    // Dujiao seive
    size_t l = 2, r;
    ll ans = 1;
    while (l <= n)
    {
        r = n / (n / l);
        ans -= (r - l + 1) * prefix_mu(n / l);
        l = r + 1;
    }
    fmu[n] = ans;
    return ans;
}

// return the largest number x that x ^ 2 <= y
ull lower_sqrt(ull y)
{
    size_t l = 1, r = INF;
    ll ans;
    while (l <= r)
    {
        size_t mid = (l + r) / 2;
        if (mid * mid <= y)
        {
            ans = mid;
            l = mid + 1;
        }
        else
        {
            r = mid - 1;
        }
    }
    return ans;
}

// Calculate the result
ull prefix_mu2(size_t n)
{
    size_t l = 1, r, bound = lower_sqrt(n);
    ll ans = 0;
    ll last_predix = 0, this_predix;
    while (l <= bound)
    {
        r = lower_sqrt(n / (n / (l * l)));
        this_predix = prefix_mu(r);
        ans += (n / (l * l)) * (this_predix - last_predix);
        l = r + 1;
        last_predix = this_predix;
    }
    return ans;
}

int main()
{
    ios::sync_with_stdio(false);
    size_t n;
    cin >> n;
    preprocess();
    cout << prefix_mu2(n) << '\n';
    return 0;
}

扩展

在这个程序中，输入比较大的 n，可以发现 \sum \mu^2(i) 总是约等于 n 的 0.6 倍，这背后是否有什么规律呢？

解析数论的研究内容之一就是数论函数的前 n 项平均的极限（或增长速率），即

\lim_{n \to \infty} \frac{1}{n}\sum_{k = 1}^n f(k).

在维基百科的 Square-free integer 条目中提到，\mu^2 的平均值的极限为 6 / \pi^2 ，即：

\lim_{n \to \infty} \frac{1}{n}\sum_{k = 1}^n [\mu(k)]^2 = \frac{6}{\pi^2} \approx 0.6079.

通过这个结论也可以验证算法在 n 较大时输出的结果。

写在最后

本题解也同步发表于知乎。

本题解的主要思想来自知乎文章 https://zhuanlan.zhihu.com/p/521413862 和超理论坛讨论 https://chaoli.club/index.php/8624 ，经过自己的思考、测试和细节补充最终得到这份题解。

本 ACM 蒟蒻只是初步涉猎初等数论，目前也只想到这一种解法，欢迎各位大佬在评论区一起交流