全球百事通！DZY Loves Math

题面

对于正整数 \(n\)，定义 \(f(n)\) 为 \(n\) 所含质因子的最大幂指数。例如 \(f(1960)=f(23×51×72)=3,f(10007)=1,f(1)=0\)。给定正整数 \(a,b\)，求下式的值：

(资料图片)

\[\sum_{i = 1}^a\sum_{j = 1}^bf(\gcd(a,b))\]

题解

在下文的推导中，为避免歧义，设 \(N = \min(a,b), M = \max(a,b)\)，显然 \(N \le M\)。

\[\displaystyle \begin{aligned} & \sum_{i = 1}^N \sum_{j = 1}^M f(\gcd(a,b))\\= & \sum_{d = 1}^N f(d) \sum_{i = 1}^N \sum_{j = 1}^M \left[\gcd(i,j) = d\right] \\= & \sum_{d = 1}^N f(d) \sum_{i = 1}^{\left\lfloor\frac{N}{d}\right\rfloor}\sum_{j = 1}^{\left\lfloor\frac{M}{d}\right\rfloor} \left[\gcd(i,j) = 1\right] \\ = & \sum_{d = 1}^N f(d) \sum_{i = 1}^{\left\lfloor\frac{N}{d}\right\rfloor}\sum_{j = 1}^{\left\lfloor\frac{M}{d}\right\rfloor} \varepsilon(\gcd(i,j)) \\= & \sum_{d = 1}^N f(d) \sum_{i = 1}^{\left\lfloor\frac{N}{d}\right\rfloor}\sum_{j = 1}^{\left\lfloor\frac{M}{d}\right\rfloor} \sum _{t \mid \gcd(i, j)} \mu(t)\\= & \sum_{d = 1}^N f(d) \sum_{i = 1}^{\left\lfloor\frac{N}{d}\right\rfloor}\sum_{j = 1}^{\left\lfloor\frac{M}{d}\right\rfloor} \sum _{t \mid i \land t \mid j} \mu(t)\\\end{aligned}\]

设 \(T = dt\)，那么：

\[\displaystyle \begin{aligned} & \sum_{d = 1}^N f(d) \sum_{i = 1}^{\left\lfloor\frac{N}{d}\right\rfloor}\sum_{j = 1}^{\left\lfloor\frac{M}{d}\right\rfloor} \sum _{t \mid i \land t \mid j} \mu(t)\\= & \sum_{T = 1}^N {\left\lfloor\frac{N}{T}\right\rfloor} {\left\lfloor\frac{M}{T}\right\rfloor} \sum_{d \mid T} f(d) \mu(\frac{T}{d})\end{aligned}\]

设 \(h = f * \mu\)，那么原式可化为：

\[\displaystyle \sum_{T = 1}^N {\left\lfloor\frac{N}{T}\right\rfloor} {\left\lfloor\frac{M}{T}\right\rfloor} h(T)\]

下面考虑求函数 \(h(n)\) 的值，首先对 \(n\) 进行质因数分解：

\[\displaystyle n = \prod_{i = 1}^m p_i^{c_i} \ ( \ p_i \in \mathbb{P}, c_i \ge 1 \ )\]

发现可以把 \(h(n)\) 写成如下形式：

\[\displaystyle h(n) = \sum_{ab = n} f(a) \mu(b)\]

考虑到莫比乌斯函数 \(\mu(n)\) 的性质：如果 \(n\) 中含有平方质因子那么 \(\mu(n) = 0\)。所以可以得出能产生贡献的 \(b\) 即满足 \(\mu(b) \ne 0\) 的 \(b\) 一定满足：

\[\displaystyle b = \prod_{i = 1}^m p_i^{d_i} \ ( \ p_i \in \mathbb{P}, _i \in \{0, 1\} \ )\]

所以可以得出 \(f(a) = l \lor f(a) - l - 1\)。

设 $l = \max \limits_{i = 1}^m c_i,k = \sum \limits_{i = 1}^m \left[ c_i = l\right] $。接下来按 \(k \ne m\) 和 \(k = m\) 两种情况分类讨论 \(h(n) 的值\)。

当 \(k \ne m\) 时，按 \(f(a) = l\) 和 \(f(a) = l - 1\) 两种子情况讨论。

当 \(f(a) = l\) 时，设在 \(k\) 个满足 \(c_i = l\) 的质数中选了 \(t\) 个，在另外 \(m - k\) 个质数中选了 \(s\) 个，那么可以得出贡献为：

\[\displaystyle \begin{aligned} & \sum_{s = 0} ^ {m - k} \sum_{t = 0}^{k - 1} \dbinom{k}{t} \times l \times (-1) ^ {s + t} \times \dbinom{m - k}{s}\\= & \sum_{t = 0}^{k - 1} \dbinom{k}{t} \times l \times (-1) ^ {t} \sum_{s = 0} ^ {m - k} (-1) ^ {s} \times 1^{m - k - s} \dbinom{m - k}{s} \\= & 0\end{aligned}\]

当 \(f(a) = l - 1\) 时，\(k\) 个满足 \(c_i = l\) 的质数中一定全部选上，设在另外 \(m - k\) 个质数中选了 \(s\) 个，那么可以得出贡献为：

\[\displaystyle \begin{aligned} & \sum_{s = 0} ^ {m - k} (l - 1) \times (-1) ^ {s} \times \dbinom{m - k}{s}\\= & \sum_{s = 0} ^ {m - k} (l - 1) \times (-1) ^ {s}\times 1^{m - k - s} \dbinom{m - k}{s} \\= & (l - 1) \times \sum_{s = 0} ^ {m - k} (-1) ^ {s}\times 1^{m - k - s} \dbinom{m - k}{s}\\= & 0\end{aligned}\]

当 \(k = m\) 时，有 \(f(a) = l\) 和 \(f(a) = l - 1\) 两种子情况，可以得出贡献为：

\[\displaystyle \begin{aligned} & \sum_{s = 0} ^ {m - 1} (l) \times (-1) ^ {s} \times \dbinom{m - k}{s} + (l - 1) \times (-1) ^ {m} \times \dbinom{m - k}{m - k}\\= & \sum_{s = 0} ^ {m} (l) \times (-1) ^ {s} \times \dbinom{m - k}{s} - 1 \times (-1) ^ m\\= & - 1 \times (-1) ^ m \\= & (-1) ^ {m + 1}\end{aligned}\]

综上可以得出 \(h(n)\) 的计算式：

\[h(n)=\begin{cases} 0 & k \ne m\\(-1)^{m + 1} & k = m\end{cases}\]

下面考虑如何高效的预处理出 \(h(n)\) 的值。

考虑一下欧拉筛在筛出合数 \(\displaystyle n = \prod_{i = 1}^m p_i^{c_i} \ ( \ p_i \in \mathbb{P}, c_i \ge 1 , p_i < p_{i + 1})\) 时的路径：

\[p_m \rightarrow p_m^2 \rightarrow p_m^3 \rightarrow \cdots \rightarrow p_m^{c_m} \rightarrow \\p_{m - 1} p_m^{c_m} \rightarrow p_{m - 1}^2 p_m^{c_m} \rightarrow \cdots n\]

也就是说一个合数被筛出的路径是按质因子从大到小的顺序筛出的，所以可以开三个数组 \(preCount,nowCount\) 和 \(factorCount\)，分别记录当前数的最小质因子的幂次，其他所有质因子的幂次和本质不同的质因子的个数。

如果在筛的过程中，设 \(i\) 为当前筛的数， \(j\) 为枚举的质数，\(t = i \cdot j\)，如果 \(j \nmid i\)，也就是说 \(j\) 不是 \(i\) 的质因子，但是其是 \(t\) 的最小质因子，所以 \(nowCount[t] = 1\)，但是 \(preCount[t]\) 的值要根据 \(nowCount[i]\) 和 \(preCount[i]\) 的值来考虑：如果两者相等，直接赋值即可；否则直接赋 \(-1\)，也就是说现在这个欧拉筛上的路径上的数的质因子幂次不可能相等了。如果 \(j \mid i\)，直接累加即可。

Code

#include typedef long long valueType;constexpr valueType maxN = 1e7 + 5;class LineSieve {public:    typedef long long valueType;    typedef std::vector container;private:    valueType size;    container minFactorList;    container primeList;    container preCount, nowCount, factorCount, data, sum;public:    explicit LineSieve(valueType _size_) : size(_size_), minFactorList(_size_ + 1),     preCount(_size_ + 1, 0),nowCount(_size_ + 1, 0),     factorCount(_size_ + 1), data(_size_ + 1),sum(_size_ + 1) {        primeList.reserve((size_t) std::floor(std::log((long double) (_size_ + 1))));        for (valueType i = 2; i <= size; ++i) {            if (minFactorList[i] == 0) {                primeList.push_back(i);                minFactorList[i] = i;                nowCount[i] = 1;                preCount[i] = 0;                factorCount[i] = 1;                data[i] = 1;            }            for (auto const &iter: primeList) {                valueType const t = i * iter;                if (t > size)                    break;                minFactorList[t] = iter;                if (i % iter == 0) {                    nowCount[t] = nowCount[i] + 1;                    preCount[t] = preCount[i];                    factorCount[t] = factorCount[i];                    break;                } else {                    nowCount[t] = 1;                    preCount[t] = (nowCount[i] == preCount[i] || preCount[i] == 0) ? nowCount[i] : -1;                    factorCount[t] = factorCount[i] + 1;                }            }        }        for (int i = 2; i <= size; ++i)            if (nowCount[i] == preCount[i] || preCount[i] == 0)                data[i] = (factorCount[i] & 1) == 1 ? 1 : -1;            else                data[i] = 0;        std::partial_sum(data.begin(), data.end(), sum.begin());    }    valueType ans(valueType x) const {        if (x > size)            throw std::range_error("Larger than Size.");        if (x < 0)            throw std::range_error("Too small.");        return sum[x];    }};int main() {    valueType T;    std::cin >> T;    LineSieve Euler(maxN);    typedef std::function solveFunction;    solveFunction solve = [&Euler](valueType N, valueType M) -> valueType {        if (N > M)            std::swap(N, M);        valueType result = 0;        valueType l = 1, r;        while (l <= N) {            r = std::min(N / (N / l), M / (M / l));            result += (Euler.ans(r) - Euler.ans(l - 1)) * (N / l) * (M / l);            l = r + 1;        }        return result;    };    for (int i = 1; i <= T; ++i) {        valueType N, M;        std::cin >> N >> M;        std::cout << solve(N, M) << "\n";    }    std::cout << std::flush;    return 0;}