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

本文同步发表于 知乎 和 洛谷。

题目

（来源：New Online OJ，可以在这个 OJ 上评测该题，也可在 洛谷 P9238 上提交评测）

题目描述

给定 $n$ 个按顺序摆好的硬币，一开始只有第 $1$ 个硬币朝下，其他硬币均朝上。
你每次操作可以选择任何一个整数 $i$ 并将所有满足 $j\mathrm{mod}i=0$ 的位置 $j$ 的硬币翻转。
求最少需要多少次操作可以让所有硬币都朝上。

输入、输出格式

输入一行包含一个整数 $n$。

输出一行包含一个整数表示最少需要的操作次数。

数据范围

• 对于 30% 的评测用例，$n\le 5\times 10^6$；
• 对于 70% 的评测用例，$n\le 10^9$；
• 对于 100% 的评测用例，$1\le n\le 10^{18}$。

样例

样例1

输入

7

输出

6

样例2

输入

1131796

输出

688042

前置知识点

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

• Möbius 函数 $\mu (n)$
• Dirichlet 卷积 $f\ast g$
• 数论分块
• Euler 筛（线性筛）求数论函数的值
• 杜教筛求数论函数前缀和

以上知识点的学习可参考 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$ 翻转，剩余的因数同样只有 $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 , … , p k } \{p_1, \dots, p_k\} {p1​,…,pk​} 的子集数），$2^k$ 是偶数，所以硬币 $n$ 已经被翻转了偶数次，依然是正面。

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

$$\sum _{i=1}^n{\mu }^2(k)$$

前缀和性质

首先证明 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^{\prime }=p_1^{\lfloor {\alpha }_1/2\rfloor }\cdots p_k^{\lfloor {\alpha }_k/2\rfloor }$，则 $d\mid n^{\prime }$，于是
$$\sum _{d^2\mid n}\mu (d)=\sum _{d\mid n^{\prime }}\mu (d)=I(n^{\prime })$$

而 $I(n^{\prime })=1$ 则说明 $n^{\prime }=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^{\prime })={\mu }^2(n)$，命题证毕。

然后据此推导 Möbius 函数平方的一个重要前缀和性质：
$$\sum _{k=1}^n{\mu }^2(k)=\sum _{k=1}^{\lfloor \sqrt{n}\rfloor }\lfloor \frac{n}{k^2}\rfloor \mu (k).$$

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

$$\begin{array}{rl}\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)\lfloor \frac{n}{d^2}\rfloor .\end{array}$$

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

数论分块

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

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

数论分块的关键是，对于给定的 $l$，确定 r = max ⁡ { i ∈ Z : ⌊ n l 2 ⌋ = ⌊ n i 2 ⌋ } r = \max \{i \in \mathbb{Z}: \left \lfloor \frac{n}{l^2} \right \rfloor = \left \lfloor \frac{n}{i^2} \right \rfloor\} r=max{i∈Z:⌊l2n​⌋=⌊i2n​⌋} 的值。设 $y=\lfloor \frac{n}{l^2}\rfloor$，这里需要确定的就是满足 $\lfloor \frac{n}{x^2}\rfloor =y$ 的最大整数 $x$。容易看出 $y\le \frac{n}{x^2}<y+1$，因此 $\frac{n}{y+1}<x^2\le \frac{n}{y}$，$x$ 的最大值即为 $\lfloor \sqrt{\frac{n}{y}}\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$ 的前缀和，否则单次使用杜教筛的时间复杂度将大于预处理复杂度。

数论分块的分块数

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

• 当 $d\le n^{1/3}$ 时，$\lfloor \frac{n}{d^2}\rfloor$ 至多有 $n^{1/3}$ 个值；
• 当 $d>n^{1/3}$ 时，$\frac{n}{d^2}\le n^{1/3}$，由于向下取整，$\lfloor \frac{n}{d^2}\rfloor$ 同样至多有 $n^{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 }$ 的区间个数。此时 $\lfloor \frac{n}{l^2}\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)$，恰好可以通过本题。

代码

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

1. 用 Euler 筛预处理 $k\le 10^7$ 的 Möbius 函数值及其前缀和；
2. 编写函数，使用杜教筛处理 $k>10^7$ 的 Möbius 函数前缀和；
3. 使用数论分块计算 ${\mu }^2$ 的前缀和。

#include <iostream>
#include <cstring>
#include <unordered_map>
#include <unordered_set>
using namespace std;

using ll = long long;

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

ll mu[MAXN + 1], smu[MAXN + 1];
bool is_prime[MAXN + 1];
int primes[MAXN];

void preprocess()
{
    // Euler seive
    constexpr size_t n = MAXN;
    size_t cnt = 0;

    memset(is_prime + 1, 1, n);
    is_prime[0] = is_prime[1] = false;
    mu[1] = 1;
    for (size_t 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];
            }
        }
    }

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

unordered_set<size_t> set; // Store the mu[i] that has been calculated
unordered_map<size_t, ll> fmu;

ll prefix_mu(size_t n)
{
    if (n <= MAXN)
    {
        return smu[n];
    }
    else if (set.find(n) != set.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;
    }
    set.insert(n);
    fmu[n] = ans;
    return ans;
}

// return the largest number x that x ^ 2 <= y
ll lower_sqrt(ll 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
ll prefix_mu2(size_t n)
{
    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));
        l = r + 1;
    }
    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$ 项平均的极限，即

$$\underset{n\to \infty }{\lim }\frac{1}{n}\sum _{k=1}^{n}f(k).$$

在维基百科的 Square-free integer 条目中提到，${\mu }^2$ 的平均值的极限为 ${\pi }^2/6$，即：
$$\underset{n\to \infty }{\lim }\frac{1}{n}\sum _{k=1}^n[\mu (k){]}^2=\frac{6}{{\pi }^2}\approx \mathrm{0.6079.}$$

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

写在最后

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

我只是初步涉猎初等数论和 ACM 数论题目，目前我也只想到这一种解法，如有更好的解法或其他问题，欢迎在评论区仪器交流~